StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
taumain_pythia_example.c
1 
10 #include "Tauola/Log.h"
11 #include "Tauola/Tauola.h"
12 #include "Tauola/TauolaHepMCEvent.h"
13 
14 //pythia header files
15 #ifdef PYTHIA8180_OR_LATER
16 #include "Pythia8/Pythia.h"
17 #include "Pythia8/Pythia8ToHepMC.h"
18 #else
19 #include "Pythia.h"
20 #include "HepMCInterface.h"
21 #endif
22 
23 //MC-TESTER header files
24 #include "Generate.h"
25 #include "HepMCEvent.H"
26 #include "Setup.H"
27 
28 #include "tauola_print_parameters.h"
29 using namespace std;
30 using namespace Pythia8;
31 using namespace Tauolapp;
32 
33 int NumberOfEvents = 10000;
34 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 Tauola (for the first several events)
40 void 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 
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::floatfield);
61  cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
62 }
63 
64 void redMinus(TauolaParticle *minus)
65 {
66  //
67  // this method can be used to redefine branching ratios in decay of tau-
68  // either generally, or specific to tau- with pointer *minus.
69  //
70  // Pointer *minus can be used to define properties of decays for taus
71  // at specific point(s) in the event tree. Example:
72  // vector<TauolaParticle*> x=minus->getMothers();
73  // and define special versions depending on x.
74  //
75  // Any combination of methods
76  // Tauola::setTauBr(int mode, double br);
77  // Tauola::setTaukle(double bra1, double brk0, double brk0b,double brks);
78  // can be called here
79 }
80 
81 void redPlus(TauolaParticle *plus)
82 {
83  //
84  // this method can be used to redefine branching ratios in decay of tau+
85  // either generally, or specific to tau+ with pointer *plus.
86  //
87  // Pointer *plus can be used to define properties of decays for tau
88  // at specific point(s) in the event tree. Example:
89  // vector<TauolaParticle*> x=plus->getMothers();
90  // and define special versions depending on x.
91  //
92  // Any combination of methods
93  // Tauola::setTauBr(int mode, double br);
94  // Tauola::setTaukle(double bra1, double brk0, double brk0b,double brks);
95  // can be called here
96 }
97 
98 int main(int argc,char **argv){
99 
100  Log::SummaryAtExit();
101 
102  // Initialization of pythia
103  Pythia pythia;
104  Event& event = pythia.event;
105 
106  // Pythia8 HepMC interface depends on Pythia8 version
107 #ifdef PYTHIA8180_OR_LATER
108  HepMC::Pythia8ToHepMC ToHepMC;
109 #else
110  HepMC::I_Pythia8 ToHepMC;
111 #endif
112 
113  //pythia.readString("HadronLevel:all = off");
114  pythia.readString("HadronLevel:Hadronize = off");
115  pythia.readString("SpaceShower:QEDshowerByL = off");
116  pythia.readString("SpaceShower:QEDshowerByQ = off");
117  pythia.readString("PartonLevel:ISR = off");
118  pythia.readString("PartonLevel:FSR = off");
119 
120  // Tauola is currently set to undecay taus. Otherwise, uncomment this line.
121  // Uncommenting it will speed up the test significantly
122  // pythia.particleData.readString("15:mayDecay = off");
123 
124  pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
125  pythia.readString("23:onMode = off");
126  pythia.readString("23:onIfAny = 15");
127  // pythia.readString("23:onIfMatch = 15 -15");
128 
129  pythia.init( 11, -11, 92.); //electron positron collisions
130 
131  // Set up Tauola
132 
133  // Set Tauola decay mode (if needed)
134  // Tauola::setSameParticleDecayMode(19); //19 and 22 contains K0
135  // Tauola::setOppositeParticleDecayMode(19); // 20 contains eta
136 
137  // Set Higgs scalar-pseudoscalar mixing angle
138  // Tauola::setHiggsScalarPseudoscalarMixingAngle(0.7853);
139  // Tauola::setHiggsScalarPseudoscalarPDG(25);
140 
141  Tauola::initialize();
142 
143  tauola_print_parameters(); // Prints TAUOLA parameters (residing inside its library): e.g. to test user interface
144 
145  // Our default units are GEV and MM, that will be outcome units after TAUOLA
146  // if HepMC unit variables are correctly set.
147  // with the following coice you can fix the units for final outcome:
148  // Tauola::setUnits(Tauola::GEV,Tauola::MM);
149  // Tauola::setUnits(Tauola::MEV,Tauola::CM);
150 
151  // Other usefull settings:
152  // Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off.
153  // Tauola::setTauLifetime(0.0); //new tau lifetime in mm
154  // Tauola::spin_correlation.setAll(false);
155 
156  // Log::LogDebug(true);
157 
158  // Tauola::setRedefineTauMinus(redMinus); // activates execution of routine redMinus in TAUOLA interface
159  // Tauola::setRedefineTauPlus(redPlus); // activates execution of routine redPlus in TAUOLA interface
160 
161  MC_Initialize();
162 
163  // Begin event loop. Generate event.
164  for (int iEvent = 0; iEvent < NumberOfEvents; ++iEvent){
165 
166  if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<endl;
167  if (!pythia.next()) continue;
168 
169  // Convert event record to HepMC
170  HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
171 
172  // Conversion needed if HepMC use different momentum units
173  // than Pythia. However, requires HepMC 2.04 or higher.
174  HepMCEvt->use_units(HepMC::Units::GEV,HepMC::Units::MM);
175 
176  ToHepMC.fill_next_event(event, HepMCEvt);
177 
178  if(iEvent<EventsToCheck)
179  {
180  cout<<" "<<endl;
181  cout<<"Momentum conservation chceck BEFORE/AFTER Tauola"<<endl;
182  checkMomentumConservationInEvent(HepMCEvt);
183  }
184 
185  // Run TAUOLA on the event
186  TauolaHepMCEvent * t_event = new TauolaHepMCEvent(HepMCEvt);
187 
188  // Since we let Pythia decay taus, we have to undecay them first.
189  t_event->undecayTaus();
190  t_event->decayTaus();
191  delete t_event;
192 
193  if(iEvent<EventsToCheck)
194  {
195  checkMomentumConservationInEvent(HepMCEvt);
196  }
197 
198  Log::Debug(5)<<"helicites = "<<Tauola::getHelPlus()<<" "<<Tauola::getHelMinus()
199  <<" electroweak wt= "<<Tauola::getEWwt()<<endl;
200 
201  // Run MC-TESTER on the event
202  HepMCEvent temp_event(*HepMCEvt,false);
203  MC_Analyze(&temp_event);
204 
205  // Print some events at the end of the run
206  if(iEvent>=NumberOfEvents-5){
207  pythia.event.list();
208  HepMCEvt->print();
209  }
210 
211  // Clean up HepMC event
212  delete HepMCEvt;
213  }
214  pythia.statistics();
215  MC_Finalize();
216 
217  // This is an access to old FORTRAN info on generated tau sample.
218  // That is why it refers to old version number (eg. 2.7) for TAUOLA.
219  //Tauola::summary();
220 }
221 
const particle iterator
Definition: GenEvent.h:464
particle_const_iterator particles_begin() const
begin particle iteration
Definition: GenEvent.h:507
double e() const
return E
Definition: SimpleVector.h:73
double px() const
return px
Definition: SimpleVector.h:70
The GenEvent class is the core of HepMC.
Definition: GenEvent.h:155
double pz() const
return pz
Definition: SimpleVector.h:72
void print(std::ostream &ostr=std::cout) const
dumps to ostr
Definition: GenEvent.cc:277
particle_const_iterator particles_end() const
end particle iteration
Definition: GenEvent.h:511
FourVector is a simple representation of a physics 4 vector.
Definition: SimpleVector.h:42
double py() const
return py
Definition: SimpleVector.h:71
void use_units(Units::MomentumUnit, Units::LengthUnit)
Definition: GenEvent.h:856