tauola_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/forZ-MEc.h"
21#include "Photos/PhotosHepMCEvent.h"
22#include "Photos/Log.h"
23
24//TAUOLA header files
25#include "Tauola/Tauola.h"
26#include "Tauola/TauolaHepMCEvent.h"
27
28using namespace std;
29using namespace Pythia8;
30using namespace Photospp;
31using namespace Tauolapp;
32
33unsigned long NumberOfEvents = 10000;
34unsigned int EventsToCheck=20;
35
36// elementary test of HepMC typically executed before
37// detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
38// similar test was performed in Fortran
39// we perform it before and after Photos (for the first several events)
40void checkMomentumConservationInEvent(HepMC::GenEvent *evt)
41{
42 //cout<<"List of stable particles: "<<endl;
43
44 double px=0.0,py=0.0,pz=0.0,e=0.0;
45
46 for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
47 p != evt->particles_end(); ++p )
48 {
49 if( (*p)->status() == 1 )
50 {
51 HepMC::FourVector m = (*p)->momentum();
52 px+=m.px();
53 py+=m.py();
54 pz+=m.pz();
55 e +=m.e();
56 //(*p)->print();
57 }
58 }
59 cout.precision(6);
60 cout.setf(ios_base::scientific, ios_base::floatfield);
61 cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
62}
63
64// Example of user function that can be passed to PHOTOS
65// for calculation of anomalous couplings in Z NLO
66double exampleAnmalousCouplingsZNLO(double qp[4],double qm[4],double ph[4],double pp[4],double pm[4],int IDE,int IDF)
67{
68 return 1.0;
69}
70
71int main(int argc,char **argv)
72{
73 HepMC::Pythia8ToHepMC ToHepMC;
74
75 // Initialization of pythia
76 Pythia pythia;
77 Event& event = pythia.event;
78
79 pythia.readFile("tauola_photos_pythia_example.conf");
80 pythia.init();
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
91 cout.setf(ios_base::fixed, ios_base::floatfield);
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 // Example of use of anomalous couplings in Z NLO with user function passed by pointer
122 //Photos::setMeCorrectionWtForZ(true);
123 //PhotosMEforZ::set_VakPol(exampleAnmalousCouplingsZNLO);
124
125 MC_Initialize();
126
127 // Begin event loop
128 for (unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
129 {
130 if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
131 if(!pythia.next()) continue;
132
133 // Convert event record to HepMC
134 HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
135 ToHepMC.fill_next_event(event, HepMCEvt);
136
137 // Run TAUOLA on the event
138 TauolaHepMCEvent * t_event = new TauolaHepMCEvent(HepMCEvt);
139
140 // We may want to undecay previously decayed taus.
141 //t_event->undecayTaus();
142 t_event->decayTaus();
143 delete t_event;
144
145 //Log::LogPhlupa(2,4);
146
147 if(iEvent<EventsToCheck)
148 {
149 cout<<" "<<endl;
150 cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl;
151 checkMomentumConservationInEvent(HepMCEvt);
152 }
153
154 // Run PHOTOS on the event
155 PhotosHepMCEvent evt(HepMCEvt);
156 evt.process();
157
158 if(iEvent<EventsToCheck)
159 {
160 checkMomentumConservationInEvent(HepMCEvt);
161 }
162
163 // Run MC-TESTER on the event
164 HepMCEvent temp_event(*HepMCEvt,false);
165 MC_Analyze(&temp_event);
166
167 // Print out last 5 events
168 if(iEvent>=NumberOfEvents-5)
169 {
170 Log::RedirectOutput(Log::Info());
171 //pythia.event.list();
172 HepMCEvt->print();
174 }
175
176 // Clean up
177 delete HepMCEvt;
178 }
179
180 Log::RedirectOutput(Log::Info());
181 pythia.stat();
183
184 MC_Finalize();
185}
static void RevertOutput()
static void RedirectOutput(void(*func)(), ostream &where=*out)
Definition Log.cxx:93
static void SummaryAtExit()