StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
runEmbeddingSimulation2014.C
1 // macro to instantiate the Geant3 from within
2 // STAR C++ framework and get the starsim prompt
3 // To use it do
4 // root4star starsim.C
5 
10 Int_t hid( Int_t z, Int_t a, Int_t l=0 )
11 {
12  // 10LZZZAAAI
13  return ( 1000000000
14  + 10000000*l
15  + 10000*z
16  + 10*a );
17 }
18 
19 class St_geant_Maker;
20 St_geant_Maker *geant_maker = 0;
21 
22 class StarGenEvent;
23 StarGenEvent *event = 0;
24 
25 class StarPrimaryMaker;
26 StarPrimaryMaker *_primary = 0;
27 
28 class StarKinematics;
30 
31 //TF1* ptDist = 0;
32 //TF1* etaDist = 0;
33 
34 TFile* _tagfile = 0; // file containing the tags
35 TTree* _tags = 0; // tree containing the tags
36 double _vxyz[3] = { 0, 0, 0 }; // event vertex
37 double _sxyz[3] = { 0, 0, 0 }; // additional smearing
38 int _runnumber = 0;
39 int _evtnumber = 0;
40 unsigned int _seed = 0;
41 int _unprimaries = 0;
42 double _eventtime = 0;
43 double _prodtime = 0;
44 double _magfield = -5.005;
45 
46 int _npart = 10; // floor number of tracks per event
47 float _fpart = 0.05; // fraction of track multiplicity to embed
48 int _pid[50]={1,2,3,8,9,11,12,14,15,10017,10007,45,50045,49,50049,37,61053,61054,62053,62054,60008,60009,95,98,97}; //geant3 particle ID
49 TString _pnm[50]={"gamma","positron","electron","pi+","pi-","K+","K-","proton","antiproton","eta","pi0","deuteron","antideuteron","he3","antihe3","D0","HyperT_2body","HyperT_bar_2body","HyperT_3body","HyperT_bar_3body","PQ1730_k0lam","PQ1730_kpluslam","K0s_piplus_piminus","lambda_p_piminus","lambdabar_pbar_piplus"}; // particle names
50 TString _part = "pi+"; // particle to simulate, default pi+
51 TString _part_save;
52 float _ptmn = 0.100; // min pT to simulate [GeV]
53 float _ptmx = 7.000; // max pT to simulate [GeV]
54 float _etamn = -2.0; // min eta to simulate
55 float _etamx = +2.0; // max eta to simulate
56 bool _rapiditymode = true; //default is flat in pseudorapidity (eta)
57 //float _mnvtx = -5.0; // min z vertex
58 //float _mxvtx = +5.0; // max z vertex
59 
60 TString _geometry = "y2014x";
61 TString DBV; // If unset, will fill from tag file "DbV20161018";
62 TString SDT = "sdt20140610";
63 
64 const int maxEvents = 1000;
65 
66 
67 TF1 *ptDist = 0;
68 TF1 *etaDist = 0;
69 
70 //______________________________________________________________________________________
71 void setRngSeed()
72 {
73 
74  TString sid = gSystem->Getenv("JOBID"); sid.Resize(8);
75  TString six = gSystem->Getenv("JOBINDEX");
76 
77  unsigned long long id = gROOT->ProcessLine(Form( "0x%s", sid.Data() ) );
78  unsigned int id1 = (0xffff&id);
79  unsigned int id2 = gROOT->ProcessLine( six.Data() );
80 
81  StarRandom::seed( _seed = id1 + id2 );
83 
84 }
85 
86 //______________________________________________________________________________________
87 
88 void geometry( TString tag, Bool_t agml=true )
89 {
90  TString cmd = "DETP GEOM "; cmd += tag;
91  if ( !geant_maker ) geant_maker = (St_geant_Maker *)chain->GetMaker("geant");
92  geant_maker -> LoadGeometry(cmd);
93  // if ( agml ) command("gexec $STAR_LIB/libxgeometry.so");
94 }
95 
96 //______________________________________________________________________________________
97 
98 void command( TString cmd )
99 {
100  if ( !geant_maker ) geant_maker = (St_geant_Maker *)chain->GetMaker("geant");
101  geant_maker -> Do( cmd );
102 }
103 
104 //______________________________________________________________________________________
105 
106 void trig( int n=1 )
107 {
108 
109  for ( int i=0; i<n; i++ ) {
110 
111  if ( i+1 > maxEvents ) break;
112 
113  // Clear the chain from the previous event
114  chain->Clear();
115 
116  // Read the ttree
117  _tags->GetEntry(i);
118 
119  cout << "_unprimaries = " << _unprimaries << endl;
120 
121  _primary->SetVertex( _vxyz[0], _vxyz[1], _vxyz[2] );
122  _primary->SetSigma(0,0,0);
123 
124  //
125  // Generate 5% of the uncorrected number of primaries recorded in the
126  // tags file. If less than minimum specified, generate the minimum
127  //
128  int npart ;
129  if ( _fpart < 1.0 ) {
130  npart = int(_unprimaries * _fpart) ;
131  if ( npart < _npart ) npart = _npart;
132  }
133  else npart = int(_fpart);
134 
135  if(i==0) _part_save=_part;
136  if(_part_save == "antideuteron" || _part_save == "he3" || _part_save == "antihe3"){
137  if(i==0)
138  _part = "deuteron"; //prime the first event with deuteron for nucleus embedding
139  else
140  _part = _part_save;
141  }
142 
143  // Print the run and vertex info
144  cout << Form("run=%i event=%i seed=%i part=%s npart=%i %i vxyz=%f %f %f",
145  _runnumber,
146  _evtnumber,
147  _seed,
148  _part.Data(),
149  npart,
150  _unprimaries,
151  _vxyz[0],_vxyz[1],_vxyz[2]) << endl;
152 
153  command( Form("RUNG %i %i", _runnumber, _evtnumber-1 ) );
154  chain->SetDateTime( int(_eventtime), int( 100000*(_eventtime-int(_eventtime)) ) );
155 
156 
157  // Standard throw things flat in pt and eta
158  if ( 0==ptDist ) kinematics->Kine( npart, _part, _ptmn, _ptmx, _etamn, _etamx );
159 
160  // Sample spectrum
161  if ( ptDist ) kinematics->Dist( npart, _part, ptDist, etaDist );
162 
163 #if 0
164  //
165  // Make sure BField can be updated on event by event
166  //
167  if (StarMagField::Instance() && StarMagField::Instance()->IsLocked()) {
168  float scale = StarMagField::Instance()->GetFactor();
169  delete StarMagField::Instance();
170  new StarMagField ( StarMagField::kMapped, scale);
171  cout << "New mag field. Don't lock, let it be updated" << endl;
172  }
173 #endif
174 
175 
176  // Generate the event
177  chain->Make();
178 
179  // // Print the event
180  _primary->event()->Print();
181  // command("gprint kine");
182  }
183 }
184 //______________________________________________________________________________________
185 
186 
187 void Kinematics()
188 {
189 
190  // gSystem->Load( "libStarGeneratorPoolPythia6_4_23.so" );
191  gSystem->Load( "libKinematics.so");
192  kinematics = new StarKinematics();
193 
194  //default is flat in eta, switch on flat in 'y'
195  if (_rapiditymode){
196  kinematics->SetAttr( "rapidity", 1 );
197  }
198 
199  _primary->AddGenerator(kinematics);
200 }
201 //______________________________________________________________________________________
202 
203 void runEmbeddingSimulation2014(
204  const char* tagfile="/star/data100/GRID/daq/2014/st_physics_adc_15161060_raw_5000043.tags.root",
205  const char* fzdfile,
206  unsigned int rngSeed=1234,
207  int nevents=-1)
208 {
209 
210  //________________________________________________________
211  //
212  // Open tagfile from where we will obtain the event vertex
213  //
214  _tagfile = TFile::Open(tagfile);
215  _tags = (TTree*) _tagfile -> Get("MoreTags");
216  //
217  _tags->SetBranchAddress( "RunId", &_runnumber );
218  _tags->SetBranchAddress( "EvtId", &_evtnumber );
219  _tags->SetBranchAddress( "VX", &_vxyz[0] );
220  _tags->SetBranchAddress( "VY", &_vxyz[1] );
221  _tags->SetBranchAddress( "VZ", &_vxyz[2] );
222  _tags->SetBranchAddress( "nRefMult", &_unprimaries );
223 
224  _tags->SetBranchAddress( "EvtTime", &_eventtime );
225  _tags->SetBranchAddress( "ProdTime", &_prodtime );
226  _tags->SetBranchAddress( "magField", &_magfield );
227 
228 
229  _tags->GetEntry(0); // read in first event
230 
231  SDT = Form("sdt%i",int(_eventtime));
232  if ( DBV == "" ) DBV = Form("dbv%i",int(_prodtime));
233 
234  // Determine maximum number of events to process
235  if ( nevents < 0 ) nevents = _tags->GetEntries();
236  cout << endl;
237  cout << "####################################################################" << endl;
238  cout << "## Processing nevents = " << nevents << endl;
239  cout << "## geometry : " << _geometry.Data() << endl;
240  cout << "## sdt timestamp: " << SDT.Data() << endl;
241  cout << "## dbv timestamp: " << DBV.Data() << endl;
242  cout << "####################################################################" << endl;
243  cout << endl;
244 
245  //
246  //________________________________________________________
247  //
248 
249  //
250  //________________________________________________________
251  //
252  // Setup the big full chain
253  //
254  //________________________________________________________
255  //
256  gROOT->ProcessLine(".L bfc.C");
257  {
258  TString simple = _geometry; simple += " ";
259  simple += SDT; simple += " ";
260  simple += DBV; simple += " ";
261  simple += " geant gstar usexgeom agml misalign newtpcalignment bigbig ";
262 
263  bfc(0, simple );
264  }
265 
266  //
267  //________________________________________________________
268  //
269  // Load in supporting libraries
270  //________________________________________________________
271  //
272  gSystem->Load( "libVMC.so");
273 
274  gSystem->Load( "StarGeneratorUtil.so" );
275  gSystem->Load( "StarGeneratorEvent.so" );
276  gSystem->Load( "StarGeneratorBase.so" );
277  gSystem->Load( "StarGeneratorDecay.so" );
278  gSystem->Load( "libPythia8_1_62.so");
279  gSystem->Load( "libMathMore.so" );
280 
281  // force gstar load/call...
282  gSystem->Load( "gstar.so" );
283  command("call gstar");
284 
285  //________________________________________________________
286  //
287  // Setup RNG seed and map all ROOT TRandom here
288  //________________________________________________________
289  //
290  //setRngSeed();
291  StarRandom::seed( _seed = rngSeed ); // but will reset based on run number and event number
293 
294  //
295  // Create the primary event generator and insert it
296  // before the geant maker
297  //
298  // StarPrimaryMaker *
299  _primary = new StarPrimaryMaker();
300  {
301  _primary -> SetFileName( "kinematics.starsim.root");
302  chain -> AddBefore( "geant", _primary );
303  }
304 
305 
306  Kinematics();
307 
308  // Particle data
310 
311 
312 #if 0
313  // Only for nucleus listed below
314  // Load STAR Particle DataBase and add the hypertriton definitions. Map to the
315  // decay modes as defined in gstar_part.g
317  pdb.AddParticle("HyperT_2body", new TParticlePDG( "HyperpT_2body", "HyperTriton --> He3 pi-", 2.99131, false, 0.0, +3.0, "hypernucleus", +hid(1,3,1), 0, 61053 ));
318  pdb.AddParticle("HyperT_bar_2body", new TParticlePDG( "HyperT_bar_2body", "AntiHyperTriton --> He3bar pi+", 2.99131, false, 0.0, -3.0, "hypernucleus", -hid(1,3,1), 0, 61054 ));
319  pdb.AddParticle("HyperT_3body", new TParticlePDG( "HyperT_3body", "HyperTriton --> d p pi-", 2.99131, false, 0.0, +3.0, "hypernucleus", +hid(1,3,1), 0, 62053 ));
320  pdb.AddParticle("HyperT_bar_3body", new TParticlePDG( "HyperT_bar_3body", "AntiHyperTriton --> dbar pbar pi+", 2.99131, false, 0.0, -3.0, "hypernucleus", -hid(1,3,1), 0, 62054 ));
321  // Hypertriton will be phase-space decayed by geant
322  pdb.AddParticle("deuteron", new TParticlePDG( "deuteron", "Deuteron", 1.876, true, 0.0, 3.0, "hion", hid(1,2,0), 0, 45 ));
323  pdb.AddParticle("antideuteron", new TParticlePDG( "antideuteron", "anti Deuteron", 1.876, true, 0.0, -3.0, "hion", -hid(1,2,0), 0, 50045 ));
324  pdb.AddParticle("he3", new TParticlePDG( "he3", "Helium-3", 2.809, true, 0.0, 6.0, "hion", hid(2,3,0), 0, 49 ));
325  pdb.AddParticle("antihe3", new TParticlePDG( "antihe3", "anti Helium-3", 2.809, true, 0.0, -6.0, "hion", -hid(2,3,0), 0, 50049 ));
326 
327 #endif
328 
329 #if 1
330  //
331  // Setup decay manager
332  //
333  StarDecayManager *decayMgr = AgUDecay::Manager();
334  StarPythia8Decayer *decayPy8 = new StarPythia8Decayer();
335  decayMgr->AddDecayer( 0, decayPy8 ); // Handle any decay requested
336  decayPy8->SetDebug(1);
337 
338  if(_part == "D0"){
339  // One can add a particle to G3 using...
340  data.AddParticleToG3( "D0", 0.186483E+01, 0.41010E-12, 0., 3, 421, 37, 0, 0 );
341  TParticlePDG* D0 = data.GetParticle("D0");
342  D0->Print();
343  //
344  // Set D0 decay to K+ pi- mode (or cc).
345  //
346  decayPy8->Set("421:onMode = off");
347  decayPy8->Set("421:onIfMatch = 211 321");
348  }
349  if(_part == "eta"){
350  data.AddParticleToG3( "eta",5.47853e-01, 0.50244E-18, 0., 3, 221, 17, 0, 0 );
351  TParticlePDG* eta = data.GetParticle("eta");
352  eta->Print();
353  //
354  // Set eta Dalitz decay to gamma e- e+ mode.
355  //
356  decayPy8->Set("221:onMode = off");
357  decayPy8->Set("221:onIfMatch = 22 11 -11");
358  }
359  if(_part == "pi0"){
360  data.AddParticleToG3( "pi0",1.34977e-01, 0.85200E-16, 0., 3, 111, 7, 0, 0 );
361  TParticlePDG* pi0 = data.GetParticle("pi0");
362  pi0->Print();
363  //
364  // Set pi0 Dalitz decay to gamma e- e+ mode.
365  //
366  decayPy8->Set("111:onMode = off");
367  decayPy8->Set("111:onIfMatch = 22 11 -11");
368  }
369  if(_part == "deuteron" || _part == "antideuteron"){
370  data.AddParticleToG3( "deuteron", 0.1876E+01, 0.10000E+16, 1., 8, 1000010020, 45, 0, 0 );
371  TParticlePDG* deuteron = data.GetParticle("deuteron");
372  deuteron->Print();
373  }
374  if(_part == "antideuteron"){
375  data.AddParticleToG3( "antideuteron", 0.1876E+01, 0.10000E+16, -1., 8, -1000010020, 50045, 0, 0 );
376  TParticlePDG* antideuteron = data.GetParticle("antideuteron");
377  antideuteron->Print();
378  }
379  if(_part == "PQ1730_k0lam") {
380  double bratio[] = { 0.5, 0.5, 0.0, 0., 0., 0.0 };
381  int mode[] = {9598,9597, 0, 0, 0, 0 };
382  data.AddParticleToG3( "PQ1730_k0lam", 0.1730E+01, 1E-16, 0., 3, 912323, 60008, bratio, mode );
383  TParticlePDG* mypart = data.GetParticle("PQ1730_k0lam");
384  mypart->Print();
385  }
386  if(_part == "PQ1730_kpluslam") {
387  double bratio[] = { 0.5, 0.5, 0.0, 0., 0., 0.0 };
388  int mode[] = {9698,9697, 0, 0, 0, 0 };
389  data.AddParticleToG3( "PQ1730_kpluslam", 0.1730E+01, 1E-16, 0., 3, 912323, 60009, bratio, mode );
390  TParticlePDG* mypart = data.GetParticle("PQ1730_kpluslam");
391  mypart->Print();
392  }
393 
394  if(_part == "K0s_piplus_piminus") {
395  double bratio[] = { 0.5, 0.5, 0.0, 0., 0., 0.0 };
396  int mode[] = {908,809, 0, 0, 0, 0 };
397  data.AddParticleToG3( "K0s_piplus_piminus", 0.4977, 0.1E-23, 0., 4, 310, 95, bratio, mode );
398  TParticlePDG* mypart = data.GetParticle("K0s_piplus_piminus");
399  mypart->Print();
400  }
401  // lambda --> p pi- g3id = 98
402  // lambdabar -- pbar pi+ g3id = 97
403  if(_part =="lambdabar_pbar_piplus") {
404  double bratio[] = { 0.5, 0.5, 0.0, 0., 0., 0.0 };
405  int mode[] = {914,1409, 0, 0, 0, 0 };
406  data.AddParticleToG3( "lambdabar_pbar_piplus", 0.1116E+1, 0.2632E-9, 0., 3, -3122, 97, bratio, mode );
407  TParticlePDG* mypart = data.GetParticle("lambdabar_pbar_piplus");
408  mypart->Print();
409  }
410  if(_part =="lambda_p_piminus") {
411  double bratio[] = { 0.5, 0.5, 0.0, 0., 0., 0.0 };
412  int mode[] = {815,1508, 0, 0, 0, 0 };
413  data.AddParticleToG3( "lambda_p_piminus", 0.1116E+1, 0.2632E-9, 0., 3, 3122, 98, bratio, mode );
414  TParticlePDG* mypart = data.GetParticle("lambda_p_piminus");
415  mypart->Print();
416  }
417 
418 
419 #endif
420 
421  //
422  // Initialize primary event generator and all sub makers
423  //
424  _primary -> Init();
425 
426  //
427  // Setup geometry and set starsim to use agusread for input
428  //
429  geometry(Form( " field=%f %s", _magfield, _geometry.Data()));
430  command("gkine -4 0");
431 
432  /*
433  TString outputFile = gSystem->BaseName( _tagfile->GetName() );
434  outputFile.ReplaceAll(".tags.root",".fzd");
435 
436  command(Form( "gfile o %s", outputFile.Data()));
437  */
438  command(Form( "gfile o %s", fzdfile));
439 
440  //
441  // Trigger on nevents
442  //
443  trig( nevents );
444 
445  command("call agexit"); // Make sure that STARSIM exits properly
446 
447 }
448 
449 //______________________________________________________________________________________
450 
451 void Validate()
452 {
453  StarGenEvent* event = _primary->event();
454  StarGenParticle* part = (*event)[1];
455  part->Print();
456  double Etotal = part->GetEnergy();
457 
458  // loop over MC particles
459  TTable* gtable = (TTable* ) chain->GetDataSet("g2t_track");
460  assert(gtable);
461  int ntable = gtable->GetNRows();
462 
463  // gtable->Print(0,ntable);
464  return;
465 
466  double Etest = 0;
467  for ( int itable=1; itable<ntable; itable++ )
468  {
469  g2t_track_st* track = (g2t_track_st*)gtable->At(itable);
470  if ( 0 == track->eg_label ) { Etest += track->e; }
471  }
472 
473  cout << "Egener = " << Etotal << endl;
474  cout << "Egeant = " << Etest << endl;
475  cout << "Violation = " << 100*(Etest / Etotal - 1) << "%" << endl;
476 
477 }
478 //______________________________________________________________________________________
479 
480 void runEmbeddingSimulation2014(
481  const int ne,
482  const unsigned int seed,
483  const char* fzdfile,
484  const char* tagfile,
485  const float mult,
486  const int pid,
487  const float ptmn,
488  const float ptmx,
489  const float etamn,
490  const float etamx,
491  const char* dbv = 0)
492 {
493  _fpart = mult;
494  _ptmn = ptmn;
495  _ptmx = ptmx;
496  _etamn = etamn;
497  _etamx = etamx;
498  if ( dbv ) DBV = dbv;
499 
500  int indx = 0;
501  for (indx=0;indx<50;indx++) {
502  if ( _pid[indx] == pid ) {
503  _part = _pnm[indx];
504  break;
505  }
506  }
507  if ( indx == 50 ) cout<<"WRONG GeantID: "<<pid<<" !!!"<<endl;
508 
509  //TFile* mHistFile = new TFile("Input/D0_weight_pt.root");
510  //If Particle name has D0, use distributions
511  if ( _part.Contains("D0") )
512  {
513  ptDist = new TF1("ptDist","(1/[0])/(1+(x/[0])^2)^[1]",_ptmn,_ptmx);
514  ptDist->SetParameter(0, 3.0);
515  ptDist->SetParameter(1, 5.0);
516  //ptDist = (TF1 *) mHistFile->Get("D0_weight_pt");
517  etaDist = new TF1("etaDist","1.0",_etamn,_etamx);
518 
519  ptDist->Print();
520  etaDist->Print();
521 
522  }
523  runEmbeddingSimulation2014( tagfile, fzdfile, seed, ne );
524 
525  //mHistFile->Close();
526 }
527 //______________________________________________________________________________________
528 
void SetSigma(Double_t sx, Double_t sy, Double_t sz, Double_t rho=0)
void Print(const Option_t *opts="head") const
Star Simple Kinematics Generator.
Yet another particle class.
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StChain.cxx:77
Float_t GetEnergy()
Get the energy.
static StarParticleData & instance()
Returns a reference to the single instance of this class.
TParticlePDG * GetParticle(const Char_t *name)
Get a particle by name.
Interface to PDG information.
void AddGenerator(StarGenerator *gener)
Connects VMC to class(es) which handle particle decays.
void Dist(Int_t ntrack, const Char_t *type, TF1 *pt, TF1 *y, TF1 *phi=0)
void AddParticle(const Char_t *name, TParticlePDG *particle)
Add a particle to the database.
Int_t Init()
Initialize generator.
virtual Long_t GetNRows() const
Returns the number of the used rows for the wrapped table.
Definition: TTable.cxx:1388
virtual Int_t Make()
Definition: StChain.cxx:110
static void seed(UInt_t s)
Definition: StarRandom.cxx:119
Base class for event records.
Definition: StarGenEvent.h:81
StarGenEvent * event()
Return a pointer to the event.
Main steering class for event generation.
void AddDecayer(Int_t pdgid, TVirtualMCDecayer *decayer)
Definition: TTable.h:48
void SetDebug(int dbg=1)
Set the debug level.
void Print(const Option_t *opts="") const
Print the particle.
const void * At(Int_t i) const
Returns a pointer to the i-th row of the table.
Definition: TTable.cxx:303
static void capture()
Capture gRandom random number generator.
Definition: StarRandom.cxx:57
Sparse class to hold track kinematics.
void Kine(Int_t ntrack, const Char_t *type="pi+,pi-,K+,K-,proton,antiproton", Double_t ptlow=0.0, Double_t pthigh=500.0, Double_t ylow=-10.0, Double_t yhigh=+10.0, Double_t philow=0.0, Double_t phihigh=TMath::TwoPi())
TParticlePDG * AddParticleToG3(const char *name, const double mass, const double lifetime, const double charge, const int tracktype, const int pdgcode, const int g3code, const double *bratio=0, const int *mode=0)
void Set(const char *cmd)
Modify pythia8 behavior.
void SetVertex(Double_t x, Double_t y, Double_t z)
Set the x, y and z vertex position.