StRoot  1
photosLCG_pythia_example.cxx
1
10 #include "Pythia8/Pythia.h"
11 #include "Pythia8Plugins/HepMC2.h"
12
14 #include "Photos/Photos.h"
15 #include "Photos/PhotosHepMCEvent.h"
16 #include "Photos/Log.h"
17
18 using namespace std;
19 using namespace Pythia8;
20 using namespace Photospp;
21
22 bool ShowersOn=true;
23 unsigned long NumberOfEvents = 100000;
24
25 // Calculates energy ratio between (l + bar_l) and (l + bar_l + X)
26 double calculate_ratio(HepMC::GenEvent *evt, double *ratio_2)
27 {
28  double ratio = 0.0;
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
57 int 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
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
76  Log::SummaryAtExit();
77  cout.setf(ios::fixed);
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::fixed);
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 }
const particle iterator
Definition: GenEvent.h:464
particle_const_iterator particles_begin() const
begin particle iteration
Definition: GenEvent.h:507
The GenEvent class is the core of HepMC.
Definition: GenEvent.h:155
particle_const_iterator particles_end() const
end particle iteration
Definition: GenEvent.h:511