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