← Back to branch summary

~maddevelopers/mg5amcnlo/new_clustering

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/srcPythia8/Pythia82_enhanceUH.cc

  • Committer: Rikkert Frederix
  • Date: 2021-09-09 15:51:40 UTC
  • mfrom: (78.75.502 3.2.1)
  • Revision ID: frederix@physik.uzh.ch-20210909155140-rg6umfq68h6h47cf
merge with 3.2.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/srcPythia8/Pythia82_enhanceUH.cc
 
1
// Driver for Pythia 8. Reads an input file dynamically created on
 
2
// the basis of the inputs specified in MCatNLO_MadFKS_PY8.Script 
 
3
#include "Pythia8/Pythia.h"
 
4
#include "Pythia8Plugins/HepMC2.h"
 
5
#include "Pythia8Plugins/aMCatNLOHooks.h"
 
6
#include "Pythia8Plugins/CombineMatchingInput.h"
 
7
#include "HepMC/GenEvent.h"
 
8
#include "HepMC/IO_GenEvent.h"
 
9
#include "HepMC/IO_BaseClass.h"
 
10
#include "HepMC/IO_HEPEVT.h"
 
11
#include "HepMC/HEPEVT_Wrapper.h"
 
12
#include "fstream"
 
13
#include "LHEFRead.h"
 
14
 
 
15
using namespace Pythia8;
 
16
 
 
17
//MZ add UserHooks
 
18
class EnhanceHooks : public UserHooks {
 
19
 
 
20
public:
 
21
 
 
22
  // Constructor and destructor do nothing.
 
23
  EnhanceHooks() {}
 
24
  ~EnhanceHooks() {}
 
25
 
 
26
  // Enhance real-emission rate. Thus no trial-emission enhancement.
 
27
  bool canEnhanceEmission() { return true;}
 
28
  bool canEnhanceTrial()    { return false;}
 
29
 
 
30
  // Function to return the weight enhance factor.
 
31
  double enhanceFactor(string name) {
 
32
    if (name == "isr:Q2QA") return 10.;
 
33
    return 1.0;
 
34
  }
 
35
 
 
36
  // Function to return the vetoing probability.
 
37
  double vetoProbability(string name) {
 
38
    if (name == "isr:Q2QA") return 0.1;
 
39
    return 0.0;
 
40
  }
 
41
 
 
42
};
 
43
//MZ
 
44
 
 
45
 
 
46
extern "C" {
 
47
  extern struct {
 
48
    double EVWGT;
 
49
  } cevwgt_;
 
50
}
 
51
#define cevwgt cevwgt_
 
52
 
 
53
extern "C" { 
 
54
  void pyabeg_(int&,char(*)[50]);
 
55
  void pyaend_(double&);
 
56
  void pyanal_(int&,double(*));
 
57
}
 
58
 
 
59
int main() {
 
60
  Pythia pythia;
 
61
  //MZ pt0 min
 
62
  pythia.settings.forceParm("SpaceShower:pT0Ref", 0.);
 
63
  pythia.settings.forceParm("TimeShower:pT0Ref", 0.);
 
64
  //
 
65
 
 
66
  int cwgtinfo_nn;
 
67
  char cwgtinfo_weights_info[1024][50];
 
68
  double cwgt_ww[1024];
 
69
 
 
70
 
 
71
  //MZ Set up a user hook and send it in.
 
72
  UserHooks* enhanceHooks = 0;
 
73
  enhanceHooks = new EnhanceHooks();
 
74
  pythia.setUserHooksPtr( enhanceHooks);
 
75
  //MZ
 
76
 
 
77
  string inputname="Pythia8.cmd",outputname="Pythia8.hep";
 
78
 
 
79
  pythia.readFile(inputname.c_str());
 
80
 
 
81
  //Create UserHooks pointer for the FxFX matching. Stop if it failed. Pass pointer to Pythia.
 
82
  CombineMatchingInput* combined = NULL;
 
83
  UserHooks* matching            = NULL;
 
84
 
 
85
  string filename = pythia.word("Beams:LHEF");
 
86
 
 
87
  MyReader read(filename);
 
88
  read.lhef_read_wgtsinfo_(cwgtinfo_nn,cwgtinfo_weights_info);
 
89
  pyabeg_(cwgtinfo_nn,cwgtinfo_weights_info);
 
90
 
 
91
  int nAbort=10;
 
92
  int nPrintLHA=1;
 
93
  int iAbort=0;
 
94
  int iPrintLHA=0;
 
95
  int nstep=5000;
 
96
  int iEventtot=pythia.mode("Main:numberOfEvents");
 
97
  int iEventshower=pythia.mode("Main:spareMode1");
 
98
  string evt_norm=pythia.word("Main:spareWord1");
 
99
  int iEventtot_norm=iEventtot;
 
100
  if (evt_norm != "sum"){
 
101
    iEventtot_norm=1;
 
102
  }
 
103
 
 
104
  //FxFx merging
 
105
  bool isFxFx=pythia.flag("JetMatching:doFxFx");
 
106
  if (isFxFx) {
 
107
    matching = combined->getHook(pythia);
 
108
    if (!matching) {
 
109
      std::cout << " Failed to initialise jet matching structures.\n"
 
110
                << " Program stopped.";
 
111
      return 1;
 
112
    }
 
113
    pythia.setUserHooksPtr(matching);
 
114
    int nJmax=pythia.mode("JetMatching:nJetMax");
 
115
    double Qcut=pythia.parm("JetMatching:qCut");
 
116
    double PTcut=pythia.parm("JetMatching:qCutME");
 
117
    if (Qcut <= PTcut || Qcut <= 0.) {
 
118
      std::cout << " \n";
 
119
      std::cout << "Merging scale (shower_card.dat) smaller than pTcut (run_card.dat)"
 
120
                << Qcut << " " << PTcut << "\n";
 
121
      return 0;
 
122
    }
 
123
  }
 
124
 
 
125
  pythia.init();
 
126
 
 
127
  HepMC::IO_BaseClass *_hepevtio;
 
128
  HepMC::Pythia8ToHepMC ToHepMC;
 
129
  HepMC::IO_GenEvent ascii_io(outputname.c_str(), std::ios::out);
 
130
  double nSelected;
 
131
  int nTry;
 
132
  double norm;
 
133
 
 
134
  // Cross section
 
135
  double sigmaTotal  = 0.;
 
136
  int iLHEFread=0;
 
137
 
 
138
  for (int iEvent = 0; ; ++iEvent) {
 
139
    if (!pythia.next()) {
 
140
      if (++iAbort < nAbort) continue;
 
141
      break;
 
142
    }
 
143
    //MZ compute the enhanced weight
 
144
 
 
145
    double enhance = enhanceHooks->getEnhancedEventWeight();
 
146
    
 
147
    //MZ
 
148
    // the number of events read by Pythia so far
 
149
    nSelected=double(pythia.info.nSelected());
 
150
    // normalisation factor for the default analyses defined in pyanal_
 
151
    norm=iEventtot_norm*iEvent/nSelected;
 
152
 
 
153
    if (nSelected >= iEventshower) break;
 
154
    if (pythia.info.isLHA() && iPrintLHA < nPrintLHA) {
 
155
      pythia.LHAeventList();
 
156
      pythia.info.list();
 
157
      pythia.process.list();
 
158
      pythia.event.list();
 
159
      ++iPrintLHA;
 
160
    }
 
161
 
 
162
    double evtweight = pythia.info.weight() * enhance;//MZ
 
163
    double normhepmc;
 
164
    // Add the weight of the current event to the cross section.
 
165
    normhepmc = 1. / double(iEventshower);
 
166
    if (evt_norm != "sum") {
 
167
      sigmaTotal  += evtweight*normhepmc;
 
168
    } else {
 
169
      sigmaTotal  += evtweight*normhepmc*iEventtot;
 
170
    }
 
171
 
 
172
    HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
 
173
    ToHepMC.fill_next_event( pythia, hepmcevt );
 
174
 
 
175
    //define the IO_HEPEVT
 
176
    _hepevtio = new HepMC::IO_HEPEVT;
 
177
    _hepevtio->write_event(hepmcevt);
 
178
    
 
179
    //event weight
 
180
    cevwgt.EVWGT=hepmcevt->weights()[0]*enhance;//MZ
 
181
 
 
182
    //call the FORTRAN analysis for this event. First, make sure to
 
183
    //re-synchronize the reading of the weights with the reading of
 
184
    //the event. (They get desynchronised if an event was rejected).
 
185
    nTry=pythia.info.nTried();
 
186
    for (; iLHEFread<nTry ; ++iLHEFread) {
 
187
      read.lhef_read_wgts_(cwgt_ww);
 
188
    }
 
189
    cwgt_ww[0]=cevwgt.EVWGT;
 
190
    pyanal_(cwgtinfo_nn,cwgt_ww);
 
191
 
 
192
    if (iEvent % nstep == 0 && iEvent >= 100){
 
193
      pyaend_(norm);
 
194
    }
 
195
    delete hepmcevt;
 
196
  }
 
197
  pyaend_(norm);
 
198
 
 
199
  pythia.stat();
 
200
  if (isFxFx){
 
201
    std::cout << " \n";
 
202
    std::cout << "*********************************************************************** \n";
 
203
    std::cout << "*********************************************************************** \n";
 
204
    std::cout << "Cross section, including FxFx merging is: "
 
205
              << sigmaTotal << "\n";
 
206
    std::cout << "*********************************************************************** \n";
 
207
    std::cout << "*********************************************************************** \n";
 
208
  }
 
209
 
 
210
  return 0;
 
211
}