single_photos_gun_example.cxx
1/**
2 * Example of use of processParticle routine.
3 * Pythia events are generated and photos used on first tau+ found.
4 *
5 * @author Tomasz Przedzinski
6 * @date 17 July 2010
7 */
8
9//pythia header files
10#include "Pythia8/Pythia.h"
11#include "Pythia8Plugins/HepMC2.h"
12
13//PHOTOS header files
14#include "Photos/Photos.h"
15#include "Photos/PhotosHepMCParticle.h"
16#include "Photos/Log.h"
17
18using namespace std;
19using namespace Pythia8;
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(int argc,char **argv)
53{
54 // Initialization of pythia
55 HepMC::Pythia8ToHepMC ToHepMC;
56 Pythia pythia;
57 Event& event = pythia.event;
58
59 pythia.readFile("single_photos_gun_example.conf");
60 pythia.init();
61
64
65 int NumberOfEvents = 10000;
66 if(argc>1) NumberOfEvents=atoi(argv[1]);
67
68 int photonAdded=0,twoAdded=0,moreAdded=0,tauCount=0;
69
70 // Begin event loop. Generate event.
71 for (int iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
72 {
73 if(iEvent%(NumberOfEvents/10)==0) Log::Info()<<iEvent<<endl;
74 if(!pythia.next()) continue;
75
76 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
77 ToHepMC.fill_next_event(event, HepMCEvt);
78
79 if(iEvent<EventsToCheck)
80 {
81 cout<<" "<<endl;
82 cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl;
83 checkMomentumConservationInEvent(HepMCEvt);
84 }
85
86 // Find tau
87 HepMC::GenParticle *tau=0;
88 for(HepMC::GenEvent::vertex_const_iterator i = HepMCEvt->vertices_begin();i!=HepMCEvt->vertices_end();i++)
89 {
90 for(HepMC::GenVertex::particles_in_const_iterator p=(*i)->particles_in_const_begin();p!=(*i)->particles_in_const_end(); p++)
91 {
92 if((*p)->pdg_id()==15) tau=*p;
93 break;
94 }
95 if(tau) break;
96 }
97 if(tau)
98 {
99 tauCount++;
100 int buf = -HepMCEvt->particles_size();
101
102 // Call photos
104
105 buf+=HepMCEvt->particles_size();
106 if(buf==1) photonAdded++;
107 else if(buf==2) twoAdded++;
108 else if(buf>2) moreAdded++;
109 }
110
111 if(iEvent<EventsToCheck)
112 {
113 checkMomentumConservationInEvent(HepMCEvt);
114 }
115
116 //clean up
117 delete HepMCEvt;
118 }
119 pythia.stat();
120
121 // Print results
122 cout.precision(2);
123 cout.setf(ios_base::fixed, ios_base::floatfield);
124 cout<<endl;
125 if(tauCount==0)
126 {
127 cout<<"Something went wrong with pythia generation."<<endl;
128 cout<<"No taus were processed."<<endl<<endl;
129 return 0;
130 }
131 cout<<"Summary (single tau decay processing):"<<endl;
132 cout<<tauCount <<"\ttaus processed"<<endl;
133 cout<<photonAdded<<"\ttimes one photon added to the decay \t("<<(photonAdded*100./tauCount)<<"%)"<<endl;
134 cout<<twoAdded <<"\ttimes two photons added to the decay \t("<<(twoAdded*100./tauCount)<<"%)"<<endl;
135 cout<<moreAdded <<"\ttimes more than two photons added to the decay\t("<<(moreAdded*100./tauCount)<<"%)"<<endl<<endl;
136 cout<<"(Contrary to results from MC-Tester, these values are technical and infrared unstable)"<<endl<<endl;
137 cout<<"To proccess different number of events use:"<<endl<<" ./single_photos_gun_example <number_of_events>"<<endl<<endl;
138}
139
static void initialize()
Definition: Photos.cxx:53
static void setInfraredCutOff(double cut_off)
static void processParticle(PhotosParticle *p)
Definition: Photos.cxx:225