TEvtGen/Photos/examples/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 "Pythia.h"
  11 #include "HepMCInterface.h"
  12
  13 //PHOTOS header files
  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;
  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
  57 int main(int argc,char **argv)
  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 }