FairRoot/PandaRoot
Functions | Variables
photosLCG_pythia_example.cxx File Reference
#include "Pythia.h"
#include "HepMCInterface.h"
#include "Photos/Photos.h"
#include "Photos/PhotosHepMCEvent.h"
#include "Photos/Log.h"

Go to the source code of this file.

Functions

double calculate_ratio (HepMC::GenEvent *evt, double *ratio_2)
 
int main (int argc, char **argv)
 

Variables

bool ShowersOn =true
 
unsigned long NumberOfEvents = 100000
 

Function Documentation

double calculate_ratio ( HepMC::GenEvent *  evt,
double *  ratio_2 
)

Definition at line 26 of file photosLCG_pythia_example.cxx.

References p.

Referenced by main().

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 }
p
Double_t p
Definition: anasim.C:58
evt
int evt
Definition: checkhelixhit.C:36
int main ( int  argc,
char **  argv 
)

Definition at line 57 of file photosLCG_pythia_example.cxx.

References calculate_ratio(), evt, chigen::initialize(), NumberOfEvents, Photospp::PhotosEvent::process(), chigen::pythia::pythia, ShowersOn, and sqrt().

58 {
59  // Initialization of pythia
60  HepMC::I_Pythia8 ToHepMC;
61  Pythia pythia;
62  Event& event = pythia.event;
63  //pythia.settings.listAll();
64 
65  if(!ShowersOn)
66  {
67  //pythia.readString("HadronLevel:all = off");
68  pythia.readString("HadronLevel:Hadronize = off");
69  pythia.readString("SpaceShower:QEDshower = off");
70  pythia.readString("SpaceShower:QEDshowerByL = off");
71  pythia.readString("SpaceShower:QEDshowerByQ = off");
72  }
73  pythia.readString("PartonLevel:ISR = on");
74  pythia.readString("PartonLevel:FSR = off");
75 
76  pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
77  pythia.readString("23:onMode = off");
78  pythia.readString("23:onIfAny = 13");
79  //pythia.readString("23:onIfAny = 11");
80  pythia.init( 11, -11, 91.187); //e+ e- collisions
81 
82  Photos::initialize();
83 
84  // Turn on NLO corrections - only for PHOTOS 3.2 or higher
85 
86  //Photos::setMeCorrectionWtForZ(zNLO);
87  //Photos::maxWtInterference(4.0);
88  //Photos::iniInfo();
89 
90  Log::SummaryAtExit();
91  cout.setf(ios::fixed);
92 
93  double ratio_exp = 0.0, ratio_alpha = 0.0;
94  double ratio_exp_2 = 0.0, ratio_alpha_2 = 0.0;
95  double buf = 0.0;
96 
97  Photos::setDoubleBrem(true);
98  Photos::setExponentiation(true);
99  Photos::setInfraredCutOff(0.000001);
100 
101  Log::Info() << "---------------- First run - EXP ----------------" <<endl;
102 
103  // Begin event loop
104  for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
105  {
106  if(iEvent%10000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
107  if (!pythia.next()) continue;
108 
109  HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
110  ToHepMC.fill_next_event(event, HepMCEvt);
111  //HepMCEvt->print();
112 
113  //Log::LogPhlupa(1,3);
114 
115  // Call photos
116  PhotosHepMCEvent evt(HepMCEvt);
117  evt.process();
118 
119  ratio_exp += calculate_ratio(HepMCEvt,&buf);
120  ratio_exp_2 += buf;
121 
122  //HepMCEvt->print();
123 
124  // Clean up
125  delete HepMCEvt;
126  }
127 
128  Photos::setDoubleBrem(false);
129  Photos::setExponentiation(false);
130  Photos::setInfraredCutOff(1./91.187); // that is 1 GeV for
131  // pythia.init( 11, -11, 91.187);
132 
133  Log::Info() << "---------------- Second run - ALPHA ORDER ----------------" <<endl;
134 
135  // Begin event loop
136  for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
137  {
138  if(iEvent%10000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
139  if (!pythia.next()) continue;
140 
141  HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
142  ToHepMC.fill_next_event(event, HepMCEvt);
143  //HepMCEvt->print();
144 
145  //Log::LogPhlupa(1,3);
146 
147  // Call photos
148  PhotosHepMCEvent evt(HepMCEvt);
149  evt.process();
150 
151  ratio_alpha += calculate_ratio(HepMCEvt,&buf);
152  ratio_alpha_2 += buf;
153 
154  //HepMCEvt->print();
155 
156  // Clean up
157  delete HepMCEvt;
158  }
159 
160  pythia.statistics();
161 
162  ratio_exp = ratio_exp / (double)NumberOfEvents;
163  ratio_exp_2 = ratio_exp_2 / (double)NumberOfEvents;
164 
165  ratio_alpha = ratio_alpha / (double)NumberOfEvents;
166  ratio_alpha_2 = ratio_alpha_2 / (double)NumberOfEvents;
167 
168  double err_exp = sqrt( (ratio_exp_2 - ratio_exp * ratio_exp ) / (double)NumberOfEvents );
169  double err_alpha = sqrt( (ratio_alpha_2 - ratio_alpha * ratio_alpha) / (double)NumberOfEvents );
170 
171  cout.precision(6);
172  cout.setf(ios::fixed);
173  cout << "********************************************************" << endl;
174  cout << "* Z -> l + bar_l + gammas *" << endl;
175  cout << "* E(l + bar_l) / E(l + bar_l + gammas) ratio *" << endl;
176  cout << "* *" << endl;
177  cout << "* PHOTOS - EXP: " << ratio_exp <<" +/- "<<err_exp <<" *" << endl;
178  cout << "* PHOTOS - ALPHA ORDER: " << ratio_alpha <<" +/- "<<err_alpha<<" *" << endl;
179  cout << "********************************************************" << endl;
180 
181 }
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
ShowersOn
bool ShowersOn
Definition: photos_pythia_example.cxx:27
NumberOfEvents
int NumberOfEvents
Definition: photos_hepevt_example.cxx:16
evt
int evt
Definition: checkhelixhit.C:36
chigen::initialize
void initialize()
Definition: ChiGenContext.cxx:216
calculate_ratio
double calculate_ratio(HepMC::GenEvent *evt, double *ratio_2)
Definition: photosLCG_pythia_example.cxx:26
chigen::pythia::pythia
Pythia8::Pythia * pythia
Definition: ChiGenContext.cxx:174

Variable Documentation

unsigned long NumberOfEvents = 100000

Definition at line 23 of file photosLCG_pythia_example.cxx.

bool ShowersOn =true

Definition at line 22 of file photosLCG_pythia_example.cxx.

Generated on Wed Apr 3 2019 08:44:19 for FairRoot/PandaRoot by   doxygen 1.8.5