photos_tauola_test.cxx
1/**
2 * Main program for testing photos C++ interface.
3 * Pythia events are generated first, Tauola++ used for tau decays
4 * and photos used for FSR.
5 *
6 * @author Nadia Davidson and Tomasz Przedzinski
7 * @date 10 May 2011
8 */
9
10//Pythia header files
11#include "Pythia8/Pythia.h"
12#include "Pythia8Plugins/HepMC2.h"
13
14//MC-TESTER header files
15#include "Generate.h"
16#include "HepMCEvent.H"
17#include "Setup.H"
18
19//TAUOLA header files
20#include "Tauola/Tauola.h"
21#include "Tauola/TauolaHepMCEvent.h"
22
23//PHOTOS header files
24#include "Photos/Photos.h"
25#include "Photos/PhotosHepMCParticle.h"
26#include "Photos/PhotosHepMCEvent.h"
27#include "Photos/Log.h"
28
29using namespace std;
30using namespace Pythia8;
31using namespace Photospp;
32using namespace Tauolapp;
33
34unsigned long NumberOfEvents = 10000;
35unsigned int EventsToCheck=20;
36
37// elementary test of HepMC typically executed before
38// detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
39// similar test was performed in Fortran
40// we perform it before and after Photos (for the first several events)
41void checkMomentumConservationInEvent(HepMC::GenEvent *evt)
42{
43 //cout<<"List of stable particles: "<<endl;
44
45 double px=0.0,py=0.0,pz=0.0,e=0.0;
46
47 for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
48 p != evt->particles_end(); ++p )
49 {
50 if( (*p)->status() == 1 )
51 {
52 HepMC::FourVector m = (*p)->momentum();
53 px+=m.px();
54 py+=m.py();
55 pz+=m.pz();
56 e +=m.e();
57 //(*p)->print();
58 }
59 }
60 cout.precision(6);
61 cout.setf(ios_base::scientific, ios_base::floatfield);
62 cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
63}
64
65int main(int argc,char **argv)
66{
67
68 // Program needs at least 3 parameters
69 if(argc<4)
70 {
71 cout<<endl<<"Usage: "<<argv[0]<<" <pythia_conf> <no_events> <tauola_mode> [ <alpha_order> <ScalarNLO_mode> ]"<<endl;
72 cout<<endl<<" eg. "<<argv[0]<<" pythia_H.conf 10000 4 0 0"<<endl;
73 cout<<endl;
74 return -1;
75 }
76
77 HepMC::Pythia8ToHepMC ToHepMC;
78
79 // Initialization of pythia
80 Pythia pythia;
81 Event& event = pythia.event;
82
83 /********************************************************
84 Read input parameters from console. List of parameters:
85 1. Pythia configuration filename
86 2. Number of events
87 3. Tauola decay mode (refer to Tauola documentation)
88 4. Photos - use 1-photon mode on/off
89 5. Photos - use ScalarNLO mode on/off
90
91 Example where all input parameters are used:
92
93 ./photos_tauola_test.exe pythia_H.conf 100000 4 0 0
94 - use pythia_H.conf
95 - generate 100 000 events
96 - fix TAUOLA decay to channel 4 (RHORHO_MODE)
97 - Photos is not using 1-photon mode (default option)
98 - Photos is not in ScalarNLO mode (default option)
99 *********************************************************/
100
101 // 1. Load pythia configuration file (argv[1], from console)
102 pythia.readFile(argv[1]);
103
104 // 2. Get number of events (argv[2], from console)
105 NumberOfEvents = atoi(argv[2]);
106
107 // 3. Set Tauola decay mode (argv[3], from console)
108 // argv[3]=3 (tau => pi nu_tau) for Ztautau
109 // argv[3]=4 (tau => pi pi nu_tau) for Htautau
110 Tauola::setSameParticleDecayMode(atoi(argv[3]));
111 Tauola::setOppositeParticleDecayMode(atoi(argv[3]));
112
113 pythia.init();
114 Tauola::initialize();
116
120
121 // 4. Check if we're using 1-photon mode
122 if( argc>4 && atoi(argv[4]) )
123 {
126
127 // Set infrared cutoff to 10MeV for scale M_Z=91.187GeV or 500 GeV
128 if(atoi(argv[2])==1) Photos::setInfraredCutOff(0.01/91.187);
129 else Photos::setInfraredCutOff(0.01/500.);
130
132 }
133
134 // 5. Check if we're in ScalarNLO mode
135 if( argc>5 )
136 {
137 Tauola::setEtaK0sPi(1,1,0);
138
139 // Check if we are using NLO
140 if(atoi(argv[5])) Photos::setMeCorrectionWtForScalar(true);
141 }
142
144
145 MC_Initialize();
146
147 // Begin event loop
148 for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
149 {
150 cout.setf(ios_base::fixed, ios_base::floatfield);
151 if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
152 if(!pythia.next()) continue;
153
154 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
155 ToHepMC.fill_next_event(event, HepMCEvt);
156
157 if(iEvent<EventsToCheck)
158 {
159 cout<<" "<<endl;
160 cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl;
161 checkMomentumConservationInEvent(HepMCEvt);
162 }
163
164 // Run TAUOLA on the event
165 TauolaHepMCEvent * t_event = new TauolaHepMCEvent(HepMCEvt);
166
167 // Since we let Pythia decay taus, we have to undecay them first.
168 t_event->undecayTaus();
169 t_event->decayTaus();
170 delete t_event;
171
172 // Run PHOTOS on the event
173 PhotosHepMCEvent evt(HepMCEvt);
174 evt.process();
175
176 if(iEvent<EventsToCheck)
177 {
178 checkMomentumConservationInEvent(HepMCEvt);
179 }
180
181 // Run MC-TESTER on the event
182 HepMCEvent temp_event(*HepMCEvt,false);
183 MC_Analyze(&temp_event);
184
185 //clean up
186 delete HepMCEvt;
187 }
188 pythia.stat();
189 MC_Finalize();
190}
static void SummaryAtExit()
static void initialize()
Definition: Photos.cxx:53
static void setExponentiation(bool expo)
Definition: Photos.cxx:403
static void setInfraredCutOff(double cut_off)
static void maxWtInterference(double interference)
static void setMeCorrectionWtForScalar(bool corr)
Definition: Photos.cxx:435
static void setDoubleBrem(bool doub)