LOEventGenerationBias: py8_mg5amc_bias_interface.cc

File py8_mg5amc_bias_interface.cc, 4.2 KB (added by Valentin, 4 years ago)
Line 
1#include "Pythia8/Pythia.h"
2#include <stdlib.h>
3#include <math.h>
4#include <vector>
5#include <sstream>
6#include "param_card.h"
7#include "run_card.h"
8#include "proc_card.h"
9
10bool initialization_done = false;
11
12Pythia8::Pythia pythia;
13// Stringstream of the events in which MG5aMC will write and PY8 read.
14std::stringstream EvtSS;
15
16void generate_header_SS(std::stringstream & headerSS,
17                                                const int &      Pythia8BeamA,
18                                                const int &      Pythia8BeamB,
19                                                const double &   PythiaBeamEnergyA,
20                                                const double &   PythiaBeamEnergyB) {
21        headerSS<<"<LesHouchesEvents version=\"3.0\">\n";
22        headerSS<<"<header>\n";
23        // Writing proc_card.dat in the banner
24        headerSS<<"<MG5ProcCard>\n";
25        headerSS<<ProcCard<<"\n";       
26        headerSS<<"</MG5ProcCard>\n";
27        // Writing run_card.dat in the banner   
28        headerSS<<"<MGRunCard>\n";
29        headerSS<<RunCard<<"\n";       
30        headerSS<<"</MGRunCard>\n";
31        // Writing param_card.dat in the banner
32        headerSS<<"<slha>\n";
33        headerSS<<ParamCard<<"\n";     
34        headerSS<<"</slha>\n";
35        headerSS<<"</header>\n";
36        // Writing the init block (only beam related information is relevant here)     
37        headerSS<<"<init>\n";
38        headerSS<<Pythia8BeamA<<" "<<Pythia8BeamB<<" ";
39        headerSS<<std::scientific<<PythiaBeamEnergyA<<" "<<std::scientific<<PythiaBeamEnergyB<<" ";
40        // We fill in the rest of the init information with blanks
41        headerSS<<"0 0 -1 -1 -3 1\n";
42        headerSS<<"-1.0 0.0 1.0 1\n";
43        headerSS<<"<generator name='MadGraph5_aMC@NLO'>please cite 1405.0301</generator>\n";
44        headerSS<<"</init>\n";
45//      headerSS<<"</LesHouchesEvents>\n";
46        return;
47}
48
49extern "C" {
50//Fortran interface
51  void py8_bias_weight_(const double &   eCM,
52                                                const int &      Pythia8BeamA,
53                                                const double &   PythiaBeamEnergyA,
54                                                const int &      Pythia8BeamB,
55                                                const double &   PythiaBeamEnergyB,
56                                const char *     Pythia8EvtRecord,
57                                const double *   p,
58                                                const int &      nParticles,
59                                                const double &   MurScale,
60                                                const double &   AlphaQCD,
61                                                const double &   AlphaQED,
62                                                const int *      Pythia8ID,
63                                                const int *      Pythia8MotherOne,
64                                                const int *      Pythia8MotherTwo,
65                                                const int *      Pythia8ColorOne,
66                                                const int *      Pythia8ColorTwo,
67                                                const int *      Pythia8Status,
68                                                const int *      Pythia8Helicities,
69                                                double &         OutputBiasWeight    )
70  {
71        if (!initialization_done) {
72//     /!\ WARNING this example is tailored for p p > z > e+ e- j for now.
73
74                pythia.readString("Merging:doKtMerging=on");
75                pythia.readString("Merging:Process=pp>e+e-");
76                pythia.readString("Merging:TMS=30.0");
77                pythia.readString("Merging:nJetMax=2");
78                pythia.readString("Merging:applyVeto=off");
79                pythia.readString("Beams:readLHEFheaders=on");
80                pythia.readString("Merging:includeWeightInXsection=off");
81
82                pythia.readString("hadronlevel:all=off");
83                pythia.readString("partonlevel:all=off");
84                pythia.readString("SpaceShower:QEDshowerByL=off");
85                pythia.readString("SpaceShower:QEDshowerByQ=off");
86                pythia.readString("ProcessLevel:resonanceDecays=off");
87                pythia.readString("BeamRemnants:primordialKT=off");
88                pythia.readString("TimeShower:QEDshowerByQ=off");
89                pythia.readString("TimeShower:QEDshowerByL=off");
90                pythia.readString("partonlevel:mpi=off");
91                pythia.readString("PartonLevel:FSRinResonances=off");
92                pythia.readString("PartonLevel:Remnants=off");
93                pythia.readString("Check:event=off");
94               
95                pythia.settings.mode("Beams:frameType",4);             
96                std::stringstream headerSS;
97                generate_header_SS(headerSS, Pythia8BeamA,Pythia8BeamB,PythiaBeamEnergyA,PythiaBeamEnergyB);
98                EvtSS.str(headerSS.str());
99                pythia.init(&EvtSS, &headerSS);
100                initialization_done = true;
101        }
102
103        std::string EvtStr(Pythia8EvtRecord);
104        EvtStr = EvtStr.substr(0,EvtStr.find("</event>")+8);
105        // Restart the buffer, clearing EOF state
106        EvtSS.clear();
107        // Set the stringstream buffer to contain the event only
108        EvtSS.str(EvtStr);
109
110//      std::cout<<"---> Event transmitted to Pythiaa8 <---\n"<<EvtSS.str()<<std::endl;
111
112        // Compute the merging weight. All other shower processsing has been disabled in the initialisation above.
113        pythia.next(); 
114//      pythia.event.list();
115    // This is the CKKW Sudakov merging weight from the trial shower.
116        std::cout<<pythia.info.mergingWeight()<<std::endl;
117//      exit(0);
118
119        OutputBiasWeight = pythia.info.mergingWeight();
120  }
121}