← Back to branch summary

~maddevelopers/mg5amcnlo/simple_unlops

« back to all changes in this revision

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

  • Committer: olivier-mattelaer
  • Date: 2021-11-12 09:29:31 UTC
  • mfrom: (967.1.15 3.3.0)
  • Revision ID: olivier-mattelaer-20211112092931-4ec9qfzgxkeovqog
version 3.3.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/srcPythia8/Pythia83.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
extern "C" {
 
18
  extern struct {
 
19
    double EVWGT;
 
20
  } cevwgt_;
 
21
}
 
22
#define cevwgt cevwgt_
 
23
 
 
24
extern "C" { 
 
25
  void pyabeg_(int&,char(*)[50]);
 
26
  void pyaend_(double&);
 
27
  void pyanal_(int&,double(*));
 
28
}
 
29
 
 
30
int main() {
 
31
  Pythia pythia;
 
32
 
 
33
  int cwgtinfo_nn;
 
34
  char cwgtinfo_weights_info[1024][50];
 
35
  double cwgt_ww[1024];
 
36
 
 
37
  string inputname="Pythia8.cmd",outputname="Pythia8.hep";
 
38
 
 
39
  pythia.readFile(inputname.c_str());
 
40
 
 
41
  //Create UserHooks pointer for the FxFX matching. Stop if it failed. Pass pointer to Pythia.
 
42
  CombineMatchingInput combined;
 
43
  //UserHooks* matching            = NULL;
 
44
 
 
45
  string filename = pythia.word("Beams:LHEF");
 
46
 
 
47
  MyReader read(filename);
 
48
  read.lhef_read_wgtsinfo_(cwgtinfo_nn,cwgtinfo_weights_info);
 
49
  pyabeg_(cwgtinfo_nn,cwgtinfo_weights_info);
 
50
 
 
51
  int nAbort=10;
 
52
  int nPrintLHA=1;
 
53
  int iAbort=0;
 
54
  int iPrintLHA=0;
 
55
  int nstep=5000;
 
56
  int iEventtot=pythia.mode("Main:numberOfEvents");
 
57
  int iEventshower=pythia.mode("Main:spareMode1");
 
58
  string evt_norm=pythia.word("Main:spareWord1");
 
59
  int iEventtot_norm=iEventtot;
 
60
  if (evt_norm != "sum"){
 
61
    iEventtot_norm=1;
 
62
  }
 
63
 
 
64
  //FxFx merging
 
65
  bool isFxFx=pythia.flag("JetMatching:doFxFx");
 
66
  if (isFxFx) {
 
67
    combined.setHook(pythia);
 
68
    //matching = combined->getHook(pythia);
 
69
    //if (!matching) {
 
70
    //  std::cout << " Failed to initialise jet matching structures.\n"
 
71
    //            << " Program stopped.";
 
72
    //  return 1;
 
73
    //}
 
74
    //pythia.setUserHooksPtr(matching);
 
75
    int nJmax=pythia.mode("JetMatching:nJetMax");
 
76
    double Qcut=pythia.parm("JetMatching:qCut");
 
77
    double PTcut=pythia.parm("JetMatching:qCutME");
 
78
    if (Qcut <= PTcut || Qcut <= 0.) {
 
79
      std::cout << " \n";
 
80
      std::cout << "Merging scale (shower_card.dat) smaller than pTcut (run_card.dat)"
 
81
                << Qcut << " " << PTcut << "\n";
 
82
      return 0;
 
83
    }
 
84
  }
 
85
 
 
86
  // Initialise Pythia.
 
87
  if (!pythia.init()) {
 
88
    cout << "Error: could not initialise Pythia" << endl;
 
89
    return 0;
 
90
  };
 
91
 
 
92
  HepMC::IO_BaseClass *_hepevtio;
 
93
  HepMC::Pythia8ToHepMC ToHepMC;
 
94
  HepMC::IO_GenEvent ascii_io(outputname.c_str(), std::ios::out);
 
95
  double nSelected;
 
96
  int nTry;
 
97
  double norm;
 
98
 
 
99
  // Cross section
 
100
  double sigmaTotal  = 0.;
 
101
  int iLHEFread=0;
 
102
 
 
103
  for (int iEvent = 0; ; ++iEvent) {
 
104
    if (!pythia.next()) {
 
105
      if (++iAbort < nAbort) continue;
 
106
      break;
 
107
    }
 
108
    // the number of events read by Pythia so far
 
109
    nSelected=double(pythia.info.nSelected());
 
110
    // normalisation factor for the default analyses defined in pyanal_
 
111
    norm=iEventtot_norm*iEvent/nSelected;
 
112
 
 
113
    if (nSelected >= iEventshower) break;
 
114
    if (pythia.info.isLHA() && iPrintLHA < nPrintLHA) {
 
115
      pythia.LHAeventList();
 
116
      pythia.info.list();
 
117
      pythia.process.list();
 
118
      pythia.event.list();
 
119
      ++iPrintLHA;
 
120
    }
 
121
 
 
122
    double evtweight = pythia.info.weight();
 
123
    double normhepmc;
 
124
    // Add the weight of the current event to the cross section.
 
125
    normhepmc = 1. / double(iEventshower);
 
126
    if (evt_norm != "sum") {
 
127
      sigmaTotal  += evtweight*normhepmc;
 
128
    } else {
 
129
      sigmaTotal  += evtweight*normhepmc*iEventtot;
 
130
    }
 
131
 
 
132
    HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
 
133
    ToHepMC.fill_next_event( pythia, hepmcevt );
 
134
 
 
135
    //define the IO_HEPEVT
 
136
    _hepevtio = new HepMC::IO_HEPEVT;
 
137
    _hepevtio->write_event(hepmcevt);
 
138
    
 
139
    //event weight
 
140
    cevwgt.EVWGT=hepmcevt->weights()[0];
 
141
 
 
142
    //call the FORTRAN analysis for this event. First, make sure to
 
143
    //re-synchronize the reading of the weights with the reading of
 
144
    //the event. (They get desynchronised if an event was rejected).
 
145
    nTry=pythia.info.nTried();
 
146
    for (; iLHEFread<nTry ; ++iLHEFread) {
 
147
      read.lhef_read_wgts_(cwgt_ww);
 
148
    }
 
149
    cwgt_ww[0]=cevwgt.EVWGT;
 
150
    pyanal_(cwgtinfo_nn,cwgt_ww);
 
151
 
 
152
    if (iEvent % nstep == 0 && iEvent >= 100){
 
153
      pyaend_(norm);
 
154
    }
 
155
    delete hepmcevt;
 
156
  }
 
157
  pyaend_(norm);
 
158
 
 
159
  pythia.stat();
 
160
  if (isFxFx){
 
161
    std::cout << " \n";
 
162
    std::cout << "*********************************************************************** \n";
 
163
    std::cout << "*********************************************************************** \n";
 
164
    std::cout << "Cross section, including FxFx merging is: "
 
165
              << sigmaTotal << "\n";
 
166
    std::cout << "*********************************************************************** \n";
 
167
    std::cout << "*********************************************************************** \n";
 
168
  }
 
169
 
 
170
  return 0;
 
171
}