| 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 |
|
|---|
| 10 | bool initialization_done = false;
|
|---|
| 11 |
|
|---|
| 12 | Pythia8::Pythia pythia;
|
|---|
| 13 | // Stringstream of the events in which MG5aMC will write and PY8 read.
|
|---|
| 14 | std::stringstream EvtSS;
|
|---|
| 15 |
|
|---|
| 16 | void 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 |
|
|---|
| 49 | extern "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 | }
|
|---|
