StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StarPythia8Decayer.cxx
1 #include "StarPythia8Decayer.h"
2 #include "TLorentzVector.h"
3 #include <assert.h>
4 #include "TParticle.h"
5 #include "StarGenerator/UTIL/StarRandom.h"
6 #include "TSystem.h"
7 #include "StMessMgr.h"
8 
9 #include "ParticleData.h"
10 
11 #ifndef Pythia8_version
12 #error "Pythia8_version is not defined"
13 #endif
14 
15 using namespace Pythia8;
16 
17 // ----------------------------------------------------------------------------
18 // Remap pythia's random number generator to StarRandom
19 class PyRand : public Pythia8::RndmEngine {
20 public:
21  double flat() { return StarRandom::Instance().flat(); }
22 };
23 
24 //..................................................................................................
26  mPythia(0),
27  mOwner(false),
28  mDebug(0),
29  mRootS(510.0)
30 {
31 
32  TString path = "StRoot/StarGenerator/"; path+= Pythia8_version; path+="/share/Pythia8/xmldoc/";
33  {
34  ifstream in(path);
35  if (!in.good()) { path = "$(STAR)/"+path; }
36  path = gSystem->ExpandPathName(path.Data());
37  }
38 
39  Info(GetName(),Form("MC version is %s data at %s",Pythia8_version,path.Data()));
40  Info(GetName(),Form("Configuration files at %s",path.Data()));
41 
42  mPythia = new Pythia8::Pythia( path.Data() );
43  mPythia -> setRndmEnginePtr( new PyRand() );
44 
45  mPythia -> readString("ProcessLevel:all = off"); // switch off initial state machinery
46  mPythia -> readString("Check:event = off"); // no consistency checks against IS
47  mPythia -> init();
48 
49 }
50 //..................................................................................................
52 {
53 #define Set(x) mPythia -> readString(x)
54  Set("Beams:idA=2112");
55  Set("Beams:idB=2112");
56  Set("Beams:frameType=1");
57  Set(Form("Beams:eCM=%f",mRootS));
58  mPythia->init();
59 #undef Set
60 }
61 //..................................................................................................
62 void StarPythia8Decayer::Decay( int pdgid, TLorentzVector *_p )
63 {
64 
65  LOG_INFO << "Decay pdgid=" << pdgid << endm;
66 
67  // Clear the event from the last run
68  ClearEvent();
69  // Add the particle to the pythia stack
70  AppendParticle( pdgid, _p );
71  // Get the particle ID and allow it to decay
72  // int id = mPythia -> event[0].id();
73  // mPythia->particleData.mayDecay( id, true );
74  mPythia->next();
75 
76  // Print list of pythia lines on debug
77  if ( mDebug ) mPythia->event.list();
78 }
79 //..................................................................................................
80 int StarPythia8Decayer::ImportParticles( TClonesArray *_array )
81 {
82  // Save the decay products
83  assert(_array);
84  TClonesArray &array = *_array;
85  array.Clear();
86 
87  int nparts = 0;
88 
89  for ( int i=0; i<mPythia->event.size();i++ )
90  {
91  //if ( mPythia->event[i].id() == 90 ) continue; //??
92  new(array[nparts++]) TParticle (
93  mPythia->event[i].id(),
94  mPythia->event[i].status(),
95  mPythia->event[i].mother1(),
96  mPythia->event[i].mother2(),
97  mPythia->event[i].daughter1(),
98  mPythia->event[i].daughter2(),
99  mPythia->event[i].px(), // [GeV/c]
100  mPythia->event[i].py(), // [GeV/c]
101  mPythia->event[i].pz(), // [GeV/c]
102  mPythia->event[i].e(), // [GeV]
103  mPythia->event[i].xProd(), // [mm]
104  mPythia->event[i].yProd(), // [mm]
105  mPythia->event[i].zProd(), // [mm]
106  mPythia->event[i].tProd()); // [mm/c]
107  };
108  return nparts;
109 }
110 //..................................................................................................
111 // not implemented. complain about it. loudly.
112 void StarPythia8Decayer::SetForceDecay( int type ){ assert(0); }
113 void StarPythia8Decayer::ForceDecay(){ assert(0); }
114 float StarPythia8Decayer::GetPartialBranchingRatio( int ipdg ){ assert(0); }
115 void StarPythia8Decayer::ReadDecayTable(){ assert(0); }
116 //..................................................................................................
117 float StarPythia8Decayer::GetLifetime(int pdg)
118 {
119  // return lifetime in seconds of teh particle with PDG number pdg
120  return (mPythia->particleData.tau0(pdg) * 3.3333e-12) ;
121 }
122 //..................................................................................................
123 void StarPythia8Decayer::AppendParticle(int pdg, TLorentzVector* p)
124 {
125  // Append a particle to the stack to be decayed
126  const int status = 23;
127 
128  double pt2 = p->Perp2(); // pt squared
129  double px = p->Px();
130  double py = p->Py();
131  double pz = p->Pz(); // z-component
132  double M2 = p->M2(); // mass squared
133 
134  // Lookup particle entry in pythia8
135  const auto* pde = mPythia->particleData.particleDataEntryPtr( pdg );
136  if ( pde ) {
137  M2 = pde->m0() * pde->m0();
138  }
139 
140  double E = TMath::Sqrt( pt2 + M2 + pz*pz );
141 
142  mPythia->event.append(pdg, status, 0, 0, px, py, pz, E, TMath::Sqrt(M2) );
143 }
144 //..................................................................................................
146 {
147  mPythia->event.reset(); // Reset the event
148 }
149 //..................................................................................................
151 {
152  if (mPythia) delete mPythia; // delete
153 }
static StarRandom & Instance()
Obtain the single instance of the random number generator.
Definition: StarRandom.cxx:87
float GetPartialBranchingRatio(int pdgid)
Return the branching ratio for the spdcified PDG ID.
~StarPythia8Decayer()
Class destructor.
int ImportParticles(TClonesArray *array=0)
Returns the decay products in a TClonesArray of TParticle.
Double_t flat() const
Return a random number uniformly distributed between 0 and 1.
Definition: StarRandom.cxx:68
void Init()
Initializes the decayer.
StarPythia8Decayer()
Class constructor.
void AppendParticle(int pdgid, TLorentzVector *p=0)
Add a particle with specified PDG ID to the stack.
void ClearEvent()
Clear the event.
float GetLifetime(int pdgid)
Return teh lifetime in seconds for the specified particle.
void Set(const char *cmd)
Modify pythia8 behavior.
void Decay(int pdg, TLorentzVector *p=0)
Decays the particle specified by PDG id and lorentz vector.