photos_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 Nadia Davidson
6 * @date 6 July 2009
7 */
8
9//pythia header files
10#include "Pythia8/Pythia.h"
11#include "Pythia8Plugins/HepMC2.h"
12
13//MC-TESTER header files
14#include "Generate.h"
15#include "HepMCEvent.H"
16#include "Setup.H"
17
18//PHOTOS header files
19#include "Photos/Photos.h"
20#include "Photos/PhotosHepMCEvent.h"
21#include "Photos/Log.h"
22
23using namespace std;
24using namespace Pythia8;
25using namespace Photospp;
26
27bool ShowersOn=true;
28unsigned long NumberOfEvents = 10000;
29unsigned int EventsToCheck=20;
30
31// elementary test of HepMC typically executed before
32// detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
33// similar test was performed in Fortran
34// we perform it before and after Photos (for the first several events)
35void checkMomentumConservationInEvent(HepMC::GenEvent *evt)
36{
37 //cout<<"List of stable particles: "<<endl;
38
39 double px=0.0,py=0.0,pz=0.0,e=0.0;
40
41 for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
42 p != evt->particles_end(); ++p )
43 {
44 if( (*p)->status() == 1 )
45 {
46 HepMC::FourVector m = (*p)->momentum();
47 px+=m.px();
48 py+=m.py();
49 pz+=m.pz();
50 e +=m.e();
51 //(*p)->print();
52 }
53 }
54 cout.precision(6);
55 cout.setf(ios_base::scientific, ios_base::floatfield);
56 cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
57}
58
59/* Switch Status of History Entries
60
61 If Photos::createHistoryEntries(true,3) was called, this function changes the
62 status code of photons added by Photos and particles modified by Photos
63 to 3, switching the status of history entries to 1.
64
65 This results leaves all modifications performed by Photos as history entries,
66 while the regular entries represent original, unmodified event.
67
68 This is an example of how such operation can be performed in user analysis.
69 By default, this function is not used. The example of its use is commented
70 out in main event loop.
71
72 NOTE: The algorithm works only on stable particles and assumes that
73 there were no modifications to the order of the particles in
74 which they were written to HepMC by Photos. */
75void switch_history_entries_status(HepMC::GenEvent *evt)
76{
77 for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
78 p != evt->particles_end(); ++p )
79 {
80 if((*p)->status()==3)
81 {
82 if((*p)->pdg_id()==22) continue;
83
84 int barcode = (*p)->barcode();
85
86 HepMC::GenVertex *v = (*p)->production_vertex();
87
88 // History entries are added after photons, so we check what is the
89 // position of current particle relative to photons.
90 int position = 0;
91 int last_photon_position = -1;
92
93 for(HepMC::GenVertex::particles_out_const_iterator p2 = v->particles_out_const_begin();
94 p2 != v->particles_out_const_end(); ++p2)
95 {
96 position++;
97
98 if((*p2)->barcode()==barcode) break;
99
100 if((*p2)->pdg_id()==22) { last_photon_position=position; }
101 }
102
103 // If particle is found prior to photons - it was already processed, so skip it
104 if(last_photon_position<0) continue;
105
106 position -= last_photon_position;
107 HepMC::GenParticle *part = NULL;
108
109 // Now, find the particle that corresponds to this history entry
110 for(HepMC::GenVertex::particles_out_const_iterator p2 = v->particles_out_const_begin();
111 p2 != v->particles_out_const_end(); ++p2)
112 {
113 --position;
114
115 if (position > 0) continue;
116 else if(position == 0) part = *p2;
117 else
118 {
119 // Give all remaining photons status 3
120 if((*p2)->pdg_id()==22 ) (*p2)->set_status(3);
121
122 // Finish if there are no more photons
123 else break;
124 }
125 }
126
127 // Check if this is the particle we are looking for
128 if( part->pdg_id() != (*p)->pdg_id())
129 {
130 cout<<"switch_history_entries_status: mismatch in pdg_id of history entry"<<endl;
131 cout<<"and its corresponding particle. The algorithm does not work correctly."<<endl;
132 exit(-1);
133 }
134
135 // Skip this particle if its status is not 1
136 if(part->status()!=1) continue;
137
138 // Switch status codes of these particles
139 part->set_status(3);
140 (*p)->set_status(1);
141 }
142 }
143}
144
145int main(int argc,char **argv)
146{
147 // Initialization of pythia
148 HepMC::Pythia8ToHepMC ToHepMC;
149 Pythia pythia;
150 Event& event = pythia.event;
151 //pythia.settings.listAll();
152
153 pythia.readFile("photos_pythia_example.conf");
154 pythia.init();
155
156 MC_Initialize();
157
158 Photos::initialize();
159 //Photos::setDoubleBrem(false);
160 //Photos::setExponentiation(false);
161
162 Photos::setInfraredCutOff(0.01/91.187); // 10MeV for scale to M_Z=91.187
163 Photos::maxWtInterference(3.0);
164 //Photos::createHistoryEntries(true,3);
165
166 Photos::iniInfo();
168 cout.setf(ios_base::fixed, ios_base::floatfield);
169
170 // Begin event loop
171 for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
172 {
173 if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
174 if (!pythia.next()) continue;
175
176 // Convert event record to HepMC
177 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
178 ToHepMC.fill_next_event(event, HepMCEvt);
179 //HepMCEvt->print();
180
181 if(iEvent<EventsToCheck)
182 {
183 cout<<" "<<endl;
184 cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl;
185 checkMomentumConservationInEvent(HepMCEvt);
186 }
187
188 //Log::LogPhlupa(1,3);
189
190 // Run PHOTOS on the event
191 PhotosHepMCEvent evt(HepMCEvt);
192 evt.process();
193
194 // Uncomment to turn on switching of the status code of history entries
195 // with the regular entries for stable particles
196 //switch_history_entries_status(HepMCEvt);
197
198 if(iEvent<EventsToCheck)
199 {
200 checkMomentumConservationInEvent(HepMCEvt);
201 }
202
203 //HepMCEvt->print();
204
205 // Run MC-TESTER on the event
206 HepMCEvent temp_event(*HepMCEvt,false);
207 MC_Analyze(&temp_event);
208
209 // Print out last 5 events
210 if(iEvent>=NumberOfEvents-5) HepMCEvt->print();
211
212 // Clean up
213 delete HepMCEvt;
214 }
215 pythia.stat();
216 MC_Finalize();
217}
static void SummaryAtExit()