]>
Commit | Line | Data |
---|---|---|
f97ec6ec | 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 "Pythia.h" | |
11 | #include "HepMCInterface.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 | ||
23 | //TAUOLA header files | |
24 | #include "Tauola/Tauola.h" | |
25 | #include "Tauola/TauolaHepMCEvent.h" | |
26 | ||
27 | using namespace std; | |
28 | using namespace Pythia8; | |
29 | using namespace Photospp; | |
30 | using namespace Tauolapp; | |
31 | ||
32 | unsigned long NumberOfEvents = 10000; | |
33 | unsigned int EventsToCheck=20; | |
34 | ||
35 | // elementary test of HepMC typically executed before | |
36 | // detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html | |
37 | // similar test was performed in Fortran | |
38 | // we perform it before and after Photos (for the first several events) | |
39 | void checkMomentumConservationInEvent(HepMC::GenEvent *evt) | |
40 | { | |
41 | //cout<<"List of stable particles: "<<endl; | |
42 | ||
43 | double px=0.0,py=0.0,pz=0.0,e=0.0; | |
44 | ||
45 | for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); | |
46 | p != evt->particles_end(); ++p ) | |
47 | { | |
48 | if( (*p)->status() == 1 ) | |
49 | { | |
50 | HepMC::FourVector m = (*p)->momentum(); | |
51 | px+=m.px(); | |
52 | py+=m.py(); | |
53 | pz+=m.pz(); | |
54 | e +=m.e(); | |
55 | //(*p)->print(); | |
56 | } | |
57 | } | |
58 | cout.precision(6); | |
59 | cout.setf(ios_base::floatfield); | |
60 | cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl; | |
61 | } | |
62 | ||
63 | int main(int argc,char **argv) | |
64 | { | |
65 | HepMC::I_Pythia8 ToHepMC; | |
66 | ||
67 | // Initialization of pythia | |
68 | Pythia pythia; | |
69 | Event& event = pythia.event; | |
70 | ||
71 | pythia.readString("PartonLevel:ISR = off"); | |
72 | pythia.readString("PartonLevel:FSR = off"); | |
73 | ||
74 | pythia.readString("WeakSingleBoson:ffbar2gmZ = on"); | |
75 | pythia.readString("23:onMode = off"); | |
76 | pythia.readString("23:onIfAny = 15"); | |
77 | pythia.particleData.readString("15:mayDecay = off"); //<- uncomment for pythia+tauola | |
78 | ||
79 | //pythia.init( -2212, -2212, 14000.0); //proton proton collisions | |
80 | pythia.init( 11, -11, 91.187); //electron positron collisions | |
81 | ||
82 | // TAUOLA and PHOTOS initialization | |
83 | Tauola::initialize(); | |
84 | Photos::initialize(); | |
85 | ||
86 | Photos::setInfraredCutOff(0.01/91.187); // 10MeV for scale to M_Z=91.187 | |
87 | //Photos::setDoubleBrem(false); | |
88 | //Photos::setExponentiation(false); | |
89 | ||
90 | Log::SummaryAtExit(); | |
91 | cout.setf(ios::fixed); | |
92 | //Log::LogInfo(false) //To turn printing of last five events and pythia statistics off | |
93 | ||
94 | // Example setup - suppress processing of whole Z0 decay, | |
95 | // leaving only the Z0 -> tau+ tau- decay and whole branch starting | |
96 | // from tau- to be processed | |
97 | //Photos::suppressBremForBranch(0,23); | |
98 | //Photos::forceBremForDecay (2,23,15,-15); | |
99 | //Photos::forceBremForBranch(0,15); | |
100 | ||
101 | // Force mass of electron and positron to be 0.000511 | |
102 | //Photos::forceMass(11,0.000511); | |
103 | ||
104 | // Force mass of electron and positron to be taken | |
105 | // from event record instead of being calculated from 4-vector | |
106 | //Photos::forceMassFromEventRecord(11); | |
107 | ||
108 | // Exclude particles with given status code from being processed | |
109 | // or taken into account during momentum conservation calculation | |
110 | //Photos::ignoreParticlesWithStatus(3); | |
111 | ||
112 | // Remove status code from the list of ignored status codes | |
113 | //Photos::DeIgnoreParticlesWithStatus(3); | |
114 | ||
115 | // Force writing history decay products for vertices | |
116 | // modified i.e. with added photons. These particles will | |
117 | // have the provided status code. Photos will ignore | |
118 | // all particles with this status code. | |
119 | //Photos::createHistoryEntries(true,3); | |
120 | ||
121 | MC_Initialize(); | |
122 | ||
123 | // Begin event loop | |
124 | for (unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent) | |
125 | { | |
126 | if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl; | |
127 | if(!pythia.next()) continue; | |
128 | ||
129 | // Convert event record to HepMC | |
130 | HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent(); | |
131 | ToHepMC.fill_next_event(event, HepMCEvt); | |
132 | ||
133 | // Run TAUOLA on the event | |
134 | TauolaHepMCEvent * t_event = new TauolaHepMCEvent(HepMCEvt); | |
135 | ||
136 | // We may want to undecay previously decayed taus. | |
137 | //t_event->undecayTaus(); | |
138 | t_event->decayTaus(); | |
139 | delete t_event; | |
140 | ||
141 | //Log::LogPhlupa(2,4); | |
142 | ||
143 | if(iEvent<EventsToCheck) | |
144 | { | |
145 | cout<<" "<<endl; | |
146 | cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl; | |
147 | checkMomentumConservationInEvent(HepMCEvt); | |
148 | } | |
149 | ||
150 | // Run PHOTOS on the event | |
151 | PhotosHepMCEvent evt(HepMCEvt); | |
152 | evt.process(); | |
153 | ||
154 | if(iEvent<EventsToCheck) | |
155 | { | |
156 | checkMomentumConservationInEvent(HepMCEvt); | |
157 | } | |
158 | ||
159 | // Run MC-TESTER on the event | |
160 | HepMCEvent temp_event(*HepMCEvt,false); | |
161 | MC_Analyze(&temp_event); | |
162 | ||
163 | // Print out last 5 events | |
164 | if(iEvent>=NumberOfEvents-5) | |
165 | { | |
166 | Log::RedirectOutput(Log::Info()); | |
167 | //pythia.event.list(); | |
168 | HepMCEvt->print(); | |
169 | Log::RevertOutput(); | |
170 | } | |
171 | ||
172 | // Clean up | |
173 | delete HepMCEvt; | |
174 | } | |
175 | ||
176 | Log::RedirectOutput(Log::Info()); | |
177 | pythia.statistics(); | |
178 | Log::RevertOutput(); | |
179 | ||
180 | MC_Finalize(); | |
181 | } |