photos_test.cxx
1/**
2 * Main program for testing photos C++ interface.
3 * Pythia events are generated first and photos used for FSR.
4 *
5 * @author Nadia Davidson and Tomasz Przedzinski
6 * @date 10 May 2011
7 */
8
9//Pythia header files
10#include "Pythia8/Pythia.h"
11#include "Pythia8Plugins/HepMC2.h"
12
13//MC-TESTER header files
14#include "Generate.h"
15#include "HepMCEvent.H"
16#include "Setup.H"
17
18//PHOTOS header files
19#include "Photos/Photos.h"
20#include "Photos/PhotosHepMCEvent.h"
21#include "Photos/Log.h"
22
23using namespace std;
24using namespace Pythia8;
25using namespace Photospp;
26
27unsigned long NumberOfEvents = 10000;
28unsigned int EventsToCheck=20;
29
30// elementary test of HepMC typically executed before
31// detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
32// similar test was performed in Fortran
33// we perform it before and after Photos (for the first several events)
34void checkMomentumConservationInEvent(HepMC::GenEvent *evt)
35{
36 //cout<<"List of stable particles: "<<endl;
37
38 double px=0.0,py=0.0,pz=0.0,e=0.0;
39
40 for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
41 p != evt->particles_end(); ++p )
42 {
43 if( (*p)->status() == 1 )
44 {
45 HepMC::FourVector m = (*p)->momentum();
46 px+=m.px();
47 py+=m.py();
48 pz+=m.pz();
49 e +=m.e();
50 //(*p)->print();
51 }
52 }
53 cout.precision(6);
54 cout.setf(ios_base::scientific, ios_base::floatfield);
55 cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
56}
57
58// Finds X Y -> 6 -6 decay and converts it to 100 -> 6 -6, where 100 = X + Y
59// Method used only in test for t bar t pair production
60void fixForMctester(HepMC::GenEvent *evt)
61{
62 for(HepMC::GenEvent::particle_const_iterator p=evt->particles_begin();p!=evt->particles_end(); p++)
63 if((*p)->pdg_id()==6)
64 {
65 HepMC::GenParticle *pt = *p;
66 int id=(* pt->production_vertex()->particles_in_const_begin() )->pdg_id();
67 if(id!=21 && id!=11 && id>5) continue;
68
69 // Get first mother and add 2x second mother to it
70 HepMC::GenParticle *X = (* pt->production_vertex()->particles_in_const_begin());
71 HepMC::GenParticle *Y = (* ++(pt->production_vertex()->particles_in_const_begin()) );
72 HepMC::FourVector fX = X->momentum();
73 HepMC::FourVector fY = Y->momentum();
74 HepMC::FourVector fXY(fX.px()+fY.px(),fX.py()+fY.py(),fX.pz()+fY.pz(),fX.e()+fY.e());
75 X->set_momentum(fXY);
76 // Unique ID for MC-Tester to analyze
77 X->set_pdg_id(100);
78
79 // Set 2nd mother as decayed and delete it from production vertex
80 Y->set_status(1);
81 (* Y->production_vertex()->particles_in_const_begin())->set_status(1);
82 pt->production_vertex()->remove_particle(Y);
83 return;
84 }
85}
86
87int main(int argc,char **argv)
88{
89
90 // Program needs at least 2 parameters
91 if(argc<3)
92 {
93 cout<<endl<<"Usage: "<<argv[0]<<" <pythia_conf> <no_events> [ <special_mode> ]"<<endl;
94 cout<<endl<<" eg. "<<argv[0]<<" pythia_W.conf 0 10000 4 0"<<endl;
95 cout<<endl;
96 return -1;
97 }
98
99 HepMC::Pythia8ToHepMC ToHepMC;
100
101 // Initialization of pythia
102 Pythia pythia;
103 Event& event = pythia.event;
104
105 /********************************************************
106 Read input parameters from console. List of parameters:
107 1. Pythia configuration filename
108 2. Number of events
109 3. Special mode - default(off), ttbar, NLO
110 4. Photos - use 1-photon mode on/off
111
112 Example where all input parameters are used:
113
114 ./photos_test.exe pythia_W.conf 100000 0 0
115 - use pythia_W.conf
116 - generate 100 000 events
117 - default configuration (not using any special mode)
118 - Photos is not using 1-photon mode (default option, except for WmunuNLO and ZmumuNLO)
119 *********************************************************/
120
121 // 1. Load pythia configuration file (argv[1], from console)
122 pythia.readFile(argv[1]);
123
124 // 2. Get number of events (argv[2], from console)
125 NumberOfEvents = atoi(argv[2]);
126
127 pythia.init();
129
132
133 bool topDecays = false;
134 bool pairEmission = false;
135
136 // 3. Check if we're using any special mode
137 if(argc>3)
138 {
139 // Top decays
140 if(atoi(argv[3])==1) topDecays=true;
141 // NLO mode
142 else if(atoi(argv[3])==2)
143 {
146 //Photos::meCorrectionWtForScalar(true);
147 }
148 // Pairs emission
149 else if(atoi(argv[3])==4)
150 {
151 pairEmission = true;
152 }
153 }
154
155 // 4. Check if we're using 1-photon mode
156 if(argc>4 && atoi(argv[4])==1)
157 {
162 }
163
164 // Photon emission is turned on by default
165 // Use this flag to turn it off if needed
166 //Photos::setPhotonEmission(false);
167 Photos::setPairEmission(pairEmission);
168
169 MC_Initialize();
170
173
174 // Begin event loop
175 for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
176 {
177 cout.setf(ios_base::fixed, ios_base::floatfield);
178 if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
179
180 // --- Event no transmitted inside Photos::
181 Photos::setEventNo(iEvent);
182 // --- may be useful for temporary event specific test prints.
183
184 if (!pythia.next()) continue;
185
186 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
187 ToHepMC.fill_next_event(event, HepMCEvt);
188
189 if(iEvent<EventsToCheck)
190 {
191 cout<<" "<<endl;
192 cout<<"Momentum conservation check BEFORE/AFTER Photos"<<endl;
193 checkMomentumConservationInEvent(HepMCEvt);
194 }
195
196 // Run PHOTOS on the event
197 PhotosHepMCEvent evt(HepMCEvt);
198 evt.process();
199
200 if(iEvent<EventsToCheck)
201 {
202 checkMomentumConservationInEvent(HepMCEvt);
203 }
204
205 // Top decays - we mess with the event so MC-TESTER can work on it as in LC analysis case
206 if(topDecays) fixForMctester(HepMCEvt);
207
208 // Run MC-TESTER on the event
209 HepMCEvent temp_event(*HepMCEvt,false);
210 MC_Analyze(&temp_event);
211
212 // Clean up
213 delete HepMCEvt;
214 }
215 pythia.stat();
216 MC_Finalize();
217
218 // Additional printout for pairs
219 if( pairEmission ) {
220 // PAIR emission
221 // Test with formula 11 from UTHEP-93-0301 M. Skrzypek ...
222 const double PI=3.141592653589793238462643;
223 const double ALFINV= 137.01;
224
225 double deno=log(91/2/0.000511);
226 deno=deno*deno*(deno+log(2.))*4;
227 double delta=5;//.125;//0.125; //0.25;
228 double L=log(2*delta/0.000511)-5.0/6.0;
229 double num=0.99* 4.0/3.0*(L*L*L/3.+ (31./36.- PI*PI/6.)*L+0.5940);
230 printf (" >>> Soft pairs emission --- probability tests <<< \n");
231 printf (" Delta = %15.8f GeV (set the same in pairs.cxx): \n", delta);
232 printf (" Z->ee: \n");
233 printf (" Abslolute= %15.8f Relative to crude= %15.8f \n",num/PI/PI/ALFINV/ALFINV/0.99, num/deno);
234
235
236 deno=log(91/2/0.1056);
237 deno=deno*deno*(deno+log(2.))*4;
238 L=log(2*delta/0.1056)-5.0/6.0;
239 num=0.99* 4.0/3.0*(L*L*L/3.+ (31./36.- PI*PI/6.)*L+0.5940);
240 printf (" Z->mumu: \n");
241 printf (" Abslolute= %15.8f Relative to crude= %15.8f \n",num/PI/PI/ALFINV/ALFINV/0.99, num/deno);
242 }
243
244}
static void SummaryAtExit()
static void initialize()
Definition: Photos.cxx:53
static void setMeCorrectionWtForW(bool corr)
Definition: Photos.cxx:426
static void setExponentiation(bool expo)
Definition: Photos.cxx:403
static void setInfraredCutOff(double cut_off)
static void maxWtInterference(double interference)
static void setMeCorrectionWtForZ(bool corr)
Definition: Photos.cxx:431
static void setPairEmission(bool ifpair)
Definition: Photos.cxx:416
static void setEventNo(int iEvt)
static void setDoubleBrem(bool doub)
static void iniInfo()
Definition: Photos.cxx:181