TEvtGen/Photos/examples/photos_standalone_example.cxx

   1 /**
   2  * Example of photos usage.
   3  * Events are loaded from pre-generated set featuring Z0 -> tau+ tau- decays
   4  * and processed by photos.
   5  *
   6  * @author Tomasz Przedzinski
   7  * @date 17 July 2010
   8  */
   9
  10 //HepMC header files
  11 #include "HepMC/IO_GenEvent.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 Photospp;
  20
  21 int EventsToCheck=20;
  22
  23 // elementary test of HepMC typically executed before
  24 // detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
  25 // similar test was performed in Fortran
  26 // we perform it before and after Photos (for the first several events)
  27 void checkMomentumConservationInEvent(HepMC::GenEvent *evt)
  28 {
  29         //cout<<"List of stable particles: "<<endl;
  30
  31         double px=0.0,py=0.0,pz=0.0,e=0.0;
  32
  33         for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
  34               p != evt->particles_end(); ++p )
  35         {
  36                 if( (*p)->status() == 1 )
  37                 {
  38                         HepMC::FourVector m = (*p)->momentum();
  39                         px+=m.px();
  40                         py+=m.py();
  41                         pz+=m.pz();
  42                         e +=m.e();
  43                         //(*p)->print();
  44                 }
  45         }
  46   cout.precision(6);
  47   cout.setf(ios_base::floatfield);
  48         cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
  49 }
  50
  51 int main()
  52 {
  53         HepMC::IO_GenEvent file("photos_standalone_example.dat",std::ios::in);
  54
  55         Photos::initialize();
  56         Photos::setInfraredCutOff(0.001/200);
  57
  58         int photonAdded=0,twoAdded=0,moreAdded=0,evtCount=0;
  59         // Begin event loop. Generate event.
  60         while(true)
  61         {
  62                 // Create event
  63                 HepMC::GenEvent *HepMCEvt = new HepMC::GenEvent();
  64                 file.fill_next_event(HepMCEvt);
  65                 if(file.rdstate()) break;
  66                 evtCount++;
  67                 int buf = -HepMCEvt->particles_size();
  68
  69                 //cout << "BEFORE:"<<endl;
  70                 //HepMCEvt->print();
  71
  72                 if(evtCount<EventsToCheck)
  73                 {
  74                         cout<<"                                          "<<endl;
  75                         cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl;
  76                         checkMomentumConservationInEvent(HepMCEvt);
  77                 }
  78
  79                 // Process by photos
  80                 PhotosHepMCEvent evt(HepMCEvt);
  81                 evt.process();
  82
  83                 if(evtCount<EventsToCheck)
  84                 {
  85                         checkMomentumConservationInEvent(HepMCEvt);
  86                 }
  87
  88                 buf+=HepMCEvt->particles_size();
  89                 if(buf==1)      photonAdded++;
  90                 else if(buf==2) twoAdded++;
  91                 else if(buf>2)  moreAdded++;
  92
  93                 //cout << "AFTER:"<<endl;
  94                 //HepMCEvt->print();
  95
  96                 //clean up
  97                 delete HepMCEvt;
  98         }
  99
 100         // Print results
 101         cout.precision(2);
 102         cout.setf(ios::fixed);
 103         cout<<endl;
 104         if(evtCount==0)
 105         {
 106                 cout<<"Something went wrong with the input file: photos_standalone_example.dat"<<endl;
 107                 cout<<"No events were processed."<<endl<<endl;
 108                 return 0;
 109         }
 110         cout<<"Summary (whole event processing):"<<endl;
 111         cout<<evtCount   <<"\tevents processed"<<endl;
 112         cout<<photonAdded<<"\ttimes one photon added to the event           \t("<<(photonAdded*100./evtCount)<<"%)"<<endl;
 113         cout<<twoAdded   <<"\ttimes two photons added to the event          \t("<<(twoAdded*100./evtCount)<<"%)"<<endl;
 114         cout<<moreAdded  <<"\ttimes more than two photons added to the event\t("<<(moreAdded*100./evtCount)<<"%)"<<endl<<endl;
 115         cout<<"(Contrary to results from MC-Tester, these values are technical and infrared unstable)"<<endl<<endl;
 116 }