photosLCG_pythia_example.cxx
1/**
2 * Example of use of photos C++ interface. Pythia events are
3 * generated first and photos used for FSR.
4 *
5 * @author LCG & Tomasz Przedzinski
6 * @date 6 May 2011
7 */
8
9//pythia header files
10#include "Pythia8/Pythia.h"
11#include "Pythia8Plugins/HepMC2.h"
12
13//PHOTOS header files
14#include "Photos/Photos.h"
15#include "Photos/PhotosHepMCEvent.h"
16#include "Photos/Log.h"
17
18using namespace std;
19using namespace Pythia8;
20using namespace Photospp;
21
22bool ShowersOn=true;
23unsigned long NumberOfEvents = 100000;
24
25// Calculates energy ratio between (l + bar_l) and (l + bar_l + X)
26double calculate_ratio(HepMC::GenEvent *evt, double *ratio_2)
27{
28 double ratio = 0.0;
29 for(HepMC::GenEvent::particle_const_iterator p=evt->particles_begin();p!=evt->particles_end(); p++)
30 {
31 // Search for Z
32 if( (*p)->pdg_id() != 23 ) continue;
33
34 // Ignore this Z if it does not have at least two daughters
35 if( !(*p)->end_vertex() || (*p)->end_vertex()->particles_out_size() < 2 ) continue;
36
37 // Sum all daughters other than photons
38 double sum = 0.0;
39 for(HepMC::GenVertex::particle_iterator pp = (*p)->end_vertex()->particles_begin(HepMC::children);
40 pp != (*p)->end_vertex()->particles_end(HepMC::children);
41 ++pp)
42 {
43 // (*pp)->print();
44 if( (*pp)->pdg_id() != 22 ) sum += (*pp)->momentum().e();
45 }
46
47 // Calculate ratio and ratio^2
48 ratio = sum / (*p)->momentum().e();
49 *ratio_2 = sum*sum / ( (*p)->momentum().e() * (*p)->momentum().e() );
50
51 // Assuming only one Z decay per event
52 return ratio;
53 }
54 return 0.0;
55}
56
57int main(int argc,char **argv)
58{
59 // Initialization of pythia
60 HepMC::Pythia8ToHepMC ToHepMC;
61 Pythia pythia;
62 Event& event = pythia.event;
63 //pythia.settings.listAll();
64
65 pythia.readFile("photosLCG_pythia_example.conf");
66 pythia.init();
67
68 Photos::initialize();
69
70 // Turn on NLO corrections - only for PHOTOS 3.2 or higher
71
72 //Photos::setMeCorrectionWtForZ(zNLO);
73 //Photos::maxWtInterference(4.0);
74 //Photos::iniInfo();
75
77 cout.setf(ios_base::fixed, ios_base::floatfield);
78
79 double ratio_exp = 0.0, ratio_alpha = 0.0;
80 double ratio_exp_2 = 0.0, ratio_alpha_2 = 0.0;
81 double buf = 0.0;
82
83 Photos::setDoubleBrem(true);
84 Photos::setExponentiation(true);
85 Photos::setInfraredCutOff(0.000001);
86
87 Log::Info() << "---------------- First run - EXP ----------------" <<endl;
88
89 // Begin event loop
90 for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
91 {
92 if(iEvent%10000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
93 if (!pythia.next()) continue;
94
95 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
96 ToHepMC.fill_next_event(event, HepMCEvt);
97 //HepMCEvt->print();
98
99 //Log::LogPhlupa(1,3);
100
101 // Call photos
102 PhotosHepMCEvent evt(HepMCEvt);
103 evt.process();
104
105 ratio_exp += calculate_ratio(HepMCEvt,&buf);
106 ratio_exp_2 += buf;
107
108 //HepMCEvt->print();
109
110 // Clean up
111 delete HepMCEvt;
112 }
113
114 Photos::setDoubleBrem(false);
115 Photos::setExponentiation(false);
116 Photos::setInfraredCutOff(1./91.187); // that is 1 GeV for
117 // pythia.init( 11, -11, 91.187);
118
119 Log::Info() << "---------------- Second run - ALPHA ORDER ----------------" <<endl;
120
121 // Begin event loop
122 for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
123 {
124 if(iEvent%10000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
125 if (!pythia.next()) continue;
126
127 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
128 ToHepMC.fill_next_event(event, HepMCEvt);
129 //HepMCEvt->print();
130
131 //Log::LogPhlupa(1,3);
132
133 // Call photos
134 PhotosHepMCEvent evt(HepMCEvt);
135 evt.process();
136
137 ratio_alpha += calculate_ratio(HepMCEvt,&buf);
138 ratio_alpha_2 += buf;
139
140 //HepMCEvt->print();
141
142 // Clean up
143 delete HepMCEvt;
144 }
145
146 pythia.stat();
147
148 ratio_exp = ratio_exp / (double)NumberOfEvents;
149 ratio_exp_2 = ratio_exp_2 / (double)NumberOfEvents;
150
151 ratio_alpha = ratio_alpha / (double)NumberOfEvents;
152 ratio_alpha_2 = ratio_alpha_2 / (double)NumberOfEvents;
153
154 double err_exp = sqrt( (ratio_exp_2 - ratio_exp * ratio_exp ) / (double)NumberOfEvents );
155 double err_alpha = sqrt( (ratio_alpha_2 - ratio_alpha * ratio_alpha) / (double)NumberOfEvents );
156
157 cout.precision(6);
158 cout.setf(ios_base::fixed, ios_base::floatfield);
159 cout << "********************************************************" << endl;
160 cout << "* Z -> l + bar_l + gammas *" << endl;
161 cout << "* E(l + bar_l) / E(l + bar_l + gammas) ratio *" << endl;
162 cout << "* *" << endl;
163 cout << "* PHOTOS - EXP: " << ratio_exp <<" +/- "<<err_exp <<" *" << endl;
164 cout << "* PHOTOS - ALPHA ORDER: " << ratio_alpha <<" +/- "<<err_alpha<<" *" << endl;
165 cout << "********************************************************" << endl;
166
167}
static void SummaryAtExit()