photos_hepmc3_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 26 January 2020
8 */
9
10// HepMC3 header files
11#include "HepMC3/ReaderAscii.h"
12#include "HepMC3/GenEvent.h"
13#include "HepMC3/Print.h"
14
15// PHOTOS header files
16#include "Photos/Photos.h"
17#include "Photos/PhotosHepMC3Event.h"
18#include "Photos/Log.h"
19
20using namespace std;
21using namespace Photospp;
22
23int EventsToCheck=20;
24
25// elementary test of HepMC typically executed before
26// detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
27// similar test was performed in Fortran
28// we perform it before and after Photos (for the first several events)
29void checkMomentumConservationInEvent(HepMC3::GenEvent &evt)
30{
31 double px = 0.0;
32 double py = 0.0;
33 double pz = 0.0;
34 double e = 0.0;
35
36 for (auto p : evt.particles())
37 {
38 if (p->status() == 1)
39 {
40 HepMC3::FourVector m = p->momentum();
41 px += m.px();
42 py += m.py();
43 pz += m.pz();
44 e += m.e();
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 HepMC3::ReaderAscii input_file("photos_hepmc3_standalone_example.dat");
55
58
59 int photonAdded = 0;
60 int twoAdded = 0;
61 int moreAdded = 0;
62 int evtCount = 0;
63
64 // Begin event loop. Generate event.
65 while (!input_file.failed()) {
66
67 HepMC3::GenEvent evt(Units::GEV, Units::MM);
68
69 // Read event from input file
70 input_file.read_event(evt);
71
72 // If reading failed - exit loop
73 if (input_file.failed()) {
74 break;
75 }
76
77 evtCount++;
78
79 int buf = -evt.particles().size();
80
81 //cout << "BEFORE:"<<endl;
82 //Print::listing(evt);
83
84 if (evtCount < EventsToCheck)
85 {
86 cout << endl;
87 cout << "Momentum conservation chceck BEFORE/AFTER Photos" << endl;
88 checkMomentumConservationInEvent(evt);
89 }
90
91 // Process by photos
92 PhotosHepMC3Event photosEvent(&evt);
93 photosEvent.process();
94
95 if (evtCount < EventsToCheck)
96 {
97 checkMomentumConservationInEvent(evt);
98 }
99
100 buf += evt.particles().size();
101
102 if (buf == 1) {
103 photonAdded++;
104 }
105 else if (buf == 2) {
106 twoAdded++;
107 }
108 else if (buf > 2) {
109 moreAdded++;
110 }
111
112 //cout << "AFTER:" << endl;
113 //Print::listing(evt);
114 }
115
116 input_file.close();
117
118 // Print results
119 cout.precision(2);
120 cout.setf(ios_base::fixed, ios_base::floatfield);
121 cout << endl;
122
123 if (evtCount == 0)
124 {
125 cout<<"Something went wrong with the input file: photos_standalone_example.dat"<<endl;
126 cout<<"No events were processed."<<endl<<endl;
127 return 0;
128 }
129
130 cout << "Summary (whole event processing):" << endl;
131 cout << evtCount << "\tevents processed" << endl;
132 cout << photonAdded << "\ttimes one photon added to the event \t(" << (photonAdded*100./evtCount) << "%)" << endl;
133 cout << twoAdded << "\ttimes two photons added to the event \t(" << (twoAdded*100./evtCount) << "%)" << endl;
134 cout << moreAdded << "\ttimes more than two photons added to the event\t(" << (moreAdded*100./evtCount) << "%)" << endl << endl;
135 cout << "(Contrary to results from MC-Tester, these values are technical and infrared unstable)" << endl << endl;
136}
static void initialize()
Definition Photos.cxx:53
static void setInfraredCutOff(double cut_off)