StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StarPythia8.cxx
1 #include "StarPythia8.h"
2 ClassImp(StarPythia8);
3 
4 #include "TDatabasePDG.h"
5 #include "TParticlePDG.h"
6 
7 #include "StarGenerator/UTIL/StarRandom.h"
8 #include "StarGenerator/EVENT/StarGenPPEvent.h"
9 #include "StarGenerator/EVENT/StarGenEPEvent.h"
10 #include "TString.h"
11 #include "TSystem.h"
12 
13 #ifndef Pythia8_version
14 #error "Pythia8_version is not defined"
15 #endif
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 // ----------------------------------------------------------------------------
25 // ----------------------------------------------------------------------------
26 StarPythia8::StarPythia8(const char *name) : StarGenerator(name)
27 {
28 
29  // Create the new instance of pythia, specifying the path to the
30  // auxilliary files required by pythia.
31  // NOTE: When adding new versions of Pythia8, we need to specify
32  // the version of the code in the source code
33  TString path = "StRoot/StarGenerator/"; path+= Pythia8_version; path+="/share/Pythia8/xmldoc/";
34  {
35  ifstream in(path);
36  if (!in.good()) { path = "$(STAR)/"+path; }
37  path = gSystem->ExpandPathName(path.Data());
38  }
39 
40 
41  Info(GetName(),Form("MC version is %s data at %s",Pythia8_version,path.Data()));
42  Info(GetName(),Form("Configuration files at %s",path.Data()));
43 
44  mPythia = new Pythia8::Pythia( path.Data() );
45  mPythia -> setRndmEnginePtr( new PyRand() );
46 
47 }
48 // ----------------------------------------------------------------------------
49 // ----------------------------------------------------------------------------
50 // ----------------------------------------------------------------------------
52 {
53 
54  //
55  // Remap pythia8 to pdg particles
56  //
57  map<TString,TString> particles;
58  particles["electron"]="e-";
59  particles["proton"]="proton";
60 
61  TString myBlue = particles[mBlue]; if ( myBlue == "" ) myBlue=mBlue;
62  TString myYell = particles[mYell]; if ( myYell == "" ) myYell=mYell;
63 
64  //
65  // Initialize pythia based on the frame and registered beam momenta
66  // TODO: Switch to StarParticleDB
67  //
68  TDatabasePDG &pdg = (*TDatabasePDG::Instance());
69  TParticlePDG *blue = pdg.GetParticle(myBlue); assert(blue);
70  TParticlePDG *yell = pdg.GetParticle(myYell); assert(yell);
71  //
72  if ( mFrame == "CMS" ) InitCMS ( blue->PdgCode(), yell->PdgCode() );
73  if ( mFrame == "FIXT" ) InitFIXT( blue->PdgCode(), yell->PdgCode() );
74  if ( mFrame == "3MOM" ) Init3MOM( blue->PdgCode(), yell->PdgCode() );
75  if ( mFrame == "4MOM" ) Init4MOM( blue->PdgCode(), yell->PdgCode() );
76  if ( mFrame == "5MOM" ) Init5MOM( blue->PdgCode(), yell->PdgCode() );
77  //
78  // Setup event record based upon the beam species
79  //
80  if ( mBlue == "proton" ) mEvent = new StarGenPPEvent();
81  else assert(0); // Pythia 8 does not (yet) support e+p collisions
82  //
83  // Make several particles stable which may cross the beam pipe,
84  // and so the simulation package must be allowed to decide to
85  // decay them or not.
86  //
87  Set("111:onMode=0"); // pi0 stable to permit mother/daughter in star record
88  Set("211:onMode=0"); // pi+/- stable
89  Set("221:onMode=0"); // eta stable
90  Set("321:onMode=0"); // K+/- stable
91  Set("310:onMode=0"); // K short
92  Set("130:onMode=0"); // K long
93  Set("3122:onMode=0"); // Lambda 0
94  Set("3112:onMode=0"); // Sigma -
95  Set("3222:onMode=0"); // Sigma +
96  Set("3212:onMode=0"); // Sigma 0
97  Set("3312:onMode=0"); // Xi -
98  Set("3322:onMode=0"); // Xi 0
99  Set("3334:onMode=0"); // Omega -
100  //
101  return StMaker::Init();
102  //
103 }
104 // ----------------------------------------------------------------------------
106 {
107 
108  //
109  // Generate the event. This is where the action happens. The rest of this
110  // method is just copying data from pythia into the event record.
111  //
112  mPythia -> next();
113 
114  //
115  // Get the pythis event record
116  //
117  Pythia8::Event &event = mPythia->event;
118 
119  //
120  // Get the number of particles in the event. Pythia 8 include a "0th" entry,
121  // which represents the system in the event record. We will skip over this.
122  //
123  mNumberOfParticles = event.size() - 1;
124 
125 
126  if ( mBlue == "proton" ) FillPP( mEvent );
127  else FillEP( mEvent );
128 
129  // Loop over all particles in the event
130  for ( int idx=1; idx <= mNumberOfParticles; idx++ )
131  {
132 
133  int id = event[idx].id();
134 // int stat = event.statusHepMC(idx);
135  int stat = event[idx].statusHepMC();
136  int mother1 = event[idx].mother1();
137  int mother2 = event[idx].mother2();
138  int daughter1 = event[idx].daughter1();
139  int daughter2 = event[idx].daughter2();
140  double px = event[idx].px();
141  double py = event[idx].py();
142  double pz = event[idx].pz();
143  double energy = event[idx].e();
144  double mass = event[idx].m();
145  double vx = event[idx].xProd(); // mm
146  double vy = event[idx].yProd(); // mm
147  double vz = event[idx].zProd(); // mm
148  double vt = event[idx].tProd(); // mm/c
149 
150  mEvent -> AddParticle( stat, id, mother1, mother2, daughter1, daughter2, px, py, pz, energy, mass, vx, vy, vz, vt );
151 
152  }
153 
154  return kStOK;
155 }
156 // ----------------------------------------------------------------------------
157 // ----------------------------------------------------------------------------
158 // ----------------------------------------------------------------------------
159 void StarPythia8::FillEP( StarGenEvent *myevent )
160 {
161  //
162  // Fill event-wise information
163  //
164  StarGenEPEvent *event = (StarGenEPEvent *)myevent;
165  Pythia8::Info &info = mPythia->info;
166 
167  event -> idBlue = info.idA(); // Beam particle
168  event -> idYell = info.idB(); // Beam particle
169  event -> process = info.code();
170  event -> subprocess = (info.hasSub())? 0 : info.codeSub();
171 
172  int id = info.id1();
173  double x = info.x1();
174  double xPdf = info.pdf1();
175  int valence = info.isValence1();
176  if ( TMath::Abs(id)>6 )
177  {
178  id = info.id2();
179  x = info.x2();
180  xPdf = info.pdf2();
181  valence = info.isValence2();
182  }
183 
184  event -> idParton = id;
185  event -> xParton = x;
186  event -> xPdf = xPdf;
187 
188  event -> Q2 = info.Q2Fac(); // factorization scale
189  event -> valence = valence;
190 
191 // event -> sHat = info.sHat();
192 // event -> tHat = info.tHat();
193 // event -> uHat = info.uHat();
194 // event -> ptHat = info.pTHat();
195 // event -> thetaHat = info.thetaHat();
196 // event -> phiHat = info.phiHat();
197 
198  event -> weight = info.weight();
199 
200 }
201 // ----------------------------------------------------------------------------
202 // ----------------------------------------------------------------------------
203 // ----------------------------------------------------------------------------
204 void StarPythia8::FillPP( StarGenEvent *myevent )
205 {
206  //
207  // Fill event-wise information
208  //
209  StarGenPPEvent *event = (StarGenPPEvent *)myevent;
210  Pythia8::Info &info = mPythia->info;
211 
212  event -> idBlue = info.idA(); // Beam particle
213  event -> idYell = info.idB(); // Beam particle
214  event -> process = info.code();
215  event -> subprocess = (info.hasSub())? 0 : info.codeSub();
216 
217  event -> idParton1 = info.id1();
218  event -> idParton2 = info.id2();
219  event -> xParton1 = info.x1();
220  event -> xParton2 = info.x2();
221  event -> xPdf1 = info.pdf1();
222  event -> xPdf2 = info.pdf2();
223  event -> Q2fac = info.Q2Fac(); // factorization scale
224  event -> Q2ren = info.Q2Ren(); // renormalization scale
225  event -> valence1 = info.isValence1();
226  event -> valence2 = info.isValence2();
227 
228  event -> sHat = info.sHat();
229  event -> tHat = info.tHat();
230  event -> uHat = info.uHat();
231  event -> ptHat = info.pTHat();
232  event -> thetaHat = info.thetaHat();
233  event -> phiHat = info.phiHat();
234 
235  event -> weight = info.weight();
236 
237 }
238 // ----------------------------------------------------------------------------
239 // ----------------------------------------------------------------------------
240 // ----------------------------------------------------------------------------
242 {
243 
244  StarGenStats stats( GetName(), "Pythia 8 Run Statistics" );
245  Pythia8::Info &info = mPythia->info;
246 
247  // Print pythia statistics
248  mPythia->stat();
249 
250  stats.nTried = info.nTried();
251  stats.nSelected = info.nSelected();
252  stats.nAccepted = info.nAccepted();
253  stats.sigmaGen = info.sigmaGen();
254  stats.sigmaErr = info.sigmaErr();
255  stats.sumWeightGen = info.weightSum();
256 
257  stats.nFilterSeen = stats.nAccepted;
258  stats.nFilterAccept = stats.nAccepted;
259 
260  stats.Dump();
261 
262  // Return a copy of the class we just created
263  return stats;
264 
265 }
266 
267 
268 // ----------------------------------------------------------------------------
269 int StarPythia8::InitCMS( int blue, int yell )
270 {
271  Set( Form("Beams:idA=%i",blue) ); // +z --> towards EEMC
272  Set( Form("Beams:idB=%i",yell) ); // -z --> towards ETOF
273  Set( "Beams:frameType=1"); // 1=CMS
274  Set( Form("Beams:eCM=%f",mRootS) );
275 //mPythia->init( blue, yell, mRootS );
276  mPythia->init();
277  return 0;
278 }
279 // ----------------------------------------------------------------------------
280 int StarPythia8::InitFIXT( int blue, int yell )
281 {
282  Set( Form("Beams:idA=%i",blue) ); // +z --> towards EEMC
283  Set( Form("Beams:idB=%i",yell) ); // -z --> towards ETOF
284  Set( "Beams:frameType=2"); // 2=fixed target mode
285  double eA = (mRootS>0)? mRootS : 0.0;
286  double eB = (mRootS<0)? mRootS : 0.0;
287  Set( Form("Beams:eA=%f", eA) );
288  Set( Form("Beams:eB=%f", eB) );
289  mPythia->init();
290  return 0;
291 }
292 // ----------------------------------------------------------------------------
293 int StarPythia8::Init3MOM( int blue, int yell )
294 {
295  LOG_INFO << "Init3/4/5MOM not implemented yet" << endm;
296  assert(0);
297  // cout << "Initializing 3MOM: " << endl;
298  // mBlueMomentum.Print();
299  // mYellMomentum.Print();
300 
301 
302  // mPythia -> init( blue, yell,
303  // mBlueMomentum.X(), mBlueMomentum.Y(), mBlueMomentum.Z(),
304  // mYellMomentum.X(), mYellMomentum.Y(), mYellMomentum.Z() );
305 
306  return 0;
307 }
308 // ----------------------------------------------------------------------------
309 int StarPythia8::Init4MOM( int blue, int yell )
310 {
311  return Init3MOM( blue, yell );
312 }
313 // ----------------------------------------------------------------------------
314 int StarPythia8::Init5MOM( int blue, int yell )
315 {
316  return Init3MOM( blue, yell );
317 }
318 // ----------------------------------------------------------------------------
Double_t mRootS
CMS energy or incident beam momentum for fixed target collisions.
int Generate()
Generate one event.
static StarRandom & Instance()
Obtain the single instance of the random number generator.
Definition: StarRandom.cxx:87
TString mYell
Name of the yellow beam particle (-z)
Event record class tailored to PP kinematics.
ABC for defining event generator interfaces.
Definition: StarGenerator.h:34
StarGenEvent * mEvent
Generated event.
void FillEP(StarGenEvent *event)
(Optional) Method to fill a DIS event
Double_t flat() const
Return a random number uniformly distributed between 0 and 1.
Definition: StarRandom.cxx:68
End of run statistics for event generators.
Definition: StarGenStats.h:21
Base class for event records.
Definition: StarGenEvent.h:81
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
Event record class tailored to DIS kinemaics.
Definition: Stypes.h:40
TString mFrame
Frame of the collision, i.e. CMS, FIXT, 3MOM, 4MOM, 5MOM.
TString mBlue
Name of the blue beam particle (+z)
void FillPP(StarGenEvent *event)
(Optional) Method to fill a PP event
int Init()
Initialize the event generator.
Definition: StarPythia8.cxx:51
StarGenStats Stats()
Return end-of-run statistics.
void Set(const char *s)
Pass a string to Pythia8::Pythia::readString(), for user configuration.
Definition: StarPythia8.h:91