StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
starsim.hijing.evtgen.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 //
6 // By Y. Zhang 07/29/2014
7 // Modified from Jason 's macro
8 // Added real distributions for pT, y
9 
10 
11 class St_geant_Maker;
12 St_geant_Maker *geant_maker = 0;
13 
14 class StarGenEvent;
15 StarGenEvent *event = 0;
16 
17 class StarPrimaryMaker;
18 StarPrimaryMaker *_primary = 0;
19 
20 class StarKinematics;
22 
23 TF1 *ptDist = 0;
24 TF1 *yDist = 0;
25 
26 //Initialize the settings:
27 Float_t vx = 0.;
28 Float_t vy = 0.;
29 Float_t vz = 0.;
30 Float_t vx_sig = 0.01;
31 Float_t vy_sig = 0.01;
32 Float_t vz_sig = 2.0;
33 //Float_t minVz = -5.0;
34 //Float_t maxVz = +5.0;
35 Float_t minPt = 0.0;
36 Float_t maxPt = +20;
37 Float_t minY = -1.0;
38 Float_t maxY = +1.0;
39 
40 // ----------------------------------------------------------------------------
41 void geometry( TString tag, Bool_t agml=true )
42 {
43  TString cmd = "DETP GEOM "; cmd += tag;
44  if ( !geant_maker ) geant_maker = (St_geant_Maker *)chain->GetMaker("geant");
45  geant_maker -> LoadGeometry(cmd);
46  // if ( agml ) command("gexec $STAR_LIB/libxgeometry.so");
47 }
48 // ----------------------------------------------------------------------------
49 void command( TString cmd )
50 {
51  if ( !geant_maker ) geant_maker = (St_geant_Maker *)chain->GetMaker("geant");
52  geant_maker -> Do( cmd );
53 }
54 // ----------------------------------------------------------------------------
55 void trig( Int_t n=0 )
56 {
57  for ( Int_t i=0; i<n+1; i++ ) {
58  chain->Clear();
59  // Generate 5 D0 according to a PT and Y distribution
60  if (kinematics) kinematics->Kine( 1, "B0", 0, 48, -1.5, 1.5 );
61  //if (kinematics) kinematics->Kine( 2, "B0_bar", 0, 48, -1.5, 1.5 );
62  //if (kinematics) kinematics->Kine( 2, "B_plus", 0, 48, -1.5, 1.5 );
63  //if (kinematics) kinematics->Kine( 2, "B_minus", 0, 48, -1.5, 1.5 );
64  //if (kinematics) kinematics->Kine( 2, "D0", 0, 24, -1.5, 1.5 );
65  //if (kinematics) kinematics->Kine( 2, "D0_bar", 0, 24, -1.5, 1.5 );
66  //if (kinematics) kinematics->Kine( 2, "D_plus", 0, 24, -1.5, 1.5 );
67  //if (kinematics) kinematics->Kine( 2, "D_minus", 0, 24, -1.5, 1.5 );
68  chain->Make();
69  }
70 }
71 
72 void Kinematics()
73 {
74  // The Dalitz particle is defined in starsim in gstar_part.g --
75  //
76  // Particle Dalitz code=10007 TrkTyp=4 mass=0.135 charge=0 tlife=8.4e-17, | // pdg=100111 bratio= { 1.0,} mode= {10203,} | //
77  // The particle database does not know about this particle, so we have to add it. It is
78  // important that we do not overwrite PDG id = 111, the ID of the standard pi0. Otherwise,
79  // we will have all pi0's in the event decaying by dalitz.
80 
82 
83  data.AddParticleToG3("Jpsi", 3.096, 7.48e-21, 0., 4, 443, 160, 0, 0 );
84  data.AddParticleToG3("rho", 0.770, 4.35e-24, 0., 3, 113, 152, 0, 0 );
85  data.AddParticleToG3("rho_plus", 0.767, 4.35e-24, 1., 4, 213, 153, 0, 0 );
86  data.AddParticleToG3("rho_minus", 0.767, 4.35e-24,-1., 4, -213, 154, 0, 0 );
87  data.AddParticleToG3("D_star_plus", 2.01027,6.86e-22, 1., 4, 413, 60, 0, 0 );
88  data.AddParticleToG3("D_star_minus",2.01027,6.86e-22,-1., 4, -413, 61, 0, 0 );
89  data.AddParticleToG3("D_star_0", 2.007, 3.13e-22, 0., 4, 423, 62, 0, 0 );
90  data.AddParticleToG3("D_star_0_bar",2.007, 3.13e-22, 0., 4, -423, 63, 0, 0 );
91 
92  data.AddParticleToG3("B0", 5.2790, 1.536e-12, 0., 4, 511, 72, 0, 0 );
93  data.AddParticleToG3("B0_bar", 5.2790, 1.536e-12, 0., 4, -511, 73, 0, 0 );
94  data.AddParticleToG3("B_plus", 5.2790, 1.671e-12, 1., 4, 521, 70, 0, 0 );
95  data.AddParticleToG3("B_minus",5.2790, 1.671e-12, -1., 4, -521, 71, 0, 0 );
96  data.AddParticleToG3("D0", 1.86484,0.415e-12, 0., 3, 421, 37, 0, 0 );
97  data.AddParticleToG3("D0_bar", 1.86484,1.536e-12, 0., 3, -421, 38, 0, 0 );
98  data.AddParticleToG3("D_plus", 1.869, 1.057e-12, 1., 4, 411, 35, 0, 0 );
99  data.AddParticleToG3("D_minus",1.869, 1.057e-12, -1., 4, -411, 36, 0, 0 );
100 
101  kinematics = new StarKinematics("EvtGen Decay");
102  _primary -> AddGenerator(kinematics);
103 }
104 
105 // ----------------------------------------------------------------------------
106 void Hijing()
107 {
108  StarHijing *hijing = new StarHijing("hijing");
109  hijing->SetTitle("Hijing 1.383");
110 
111  // Setup collision frame, energy and beam species
112  hijing->SetFrame("CMS",200.0);
113  hijing->SetBlue("Au");
114  hijing->SetYell("Au");
115  //hijing->SetImpact(0.0, 30.0); // Impact parameter min/max (fm) 0. 30.
116  hijing->SetImpact(0.0, 15.0); // Minimum Bias
117  hijing->hiparnt().ihpr2(4) = 0; // Jet quenching (1=yes/0=no) 0
118  hijing->hiparnt().ihpr2(3) = 0; // Hard scattering (1=yes/0=no)
119  hijing->hiparnt().hipr1(10) = 2.0; // pT jet
120  hijing->hiparnt().ihpr2(8) = 10; // Max number of jets / nucleon
121  hijing->hiparnt().ihpr2(11) = 1; // Set baryon production
122  hijing->hiparnt().ihpr2(12) = 1; // Turn on/off decay of particles [1=recommended]
123  hijing->hiparnt().ihpr2(18) = 1; // Turn on/off B production
124  hijing->hiparnt().hipr1(7) = 5.35; // Set B production ???? Not really used... Really ????
125 
126  // For more configuration options, see the HIJING manual
127  // http://ntc0.lbl.gov/~xnwang/hijing/doc.html
128 
129  _primary -> AddGenerator(hijing);
130  _primary -> SetCuts( 1.0E-6 , -1., -2.6, 2.6 );
131 
132 }
133 // ----------------------------------------------------------------------------
134 // ----------------------------------------------------------------------------
135 // ----------------------------------------------------------------------------
136 void starsim( Int_t nevents=1, Int_t Index = 0, Int_t rngSeed=4321 )
137 {
138 
139  gROOT->ProcessLine(".L bfc.C");
140  {
141  TString simple = "y2014 geant gstar usexgeom agml ";
142  bfc(0, simple );
143  }
144 
145  gSystem->Load( "libVMC.so");
146 
147  gSystem->Load( "StarGeneratorUtil.so" );
148  gSystem->Load( "StarGeneratorEvent.so" );
149  gSystem->Load( "StarGeneratorBase.so" );
150  gSystem->Load( "StarGeneratorDecay.so" );
151  gSystem->Load( "libMathMore.so" );
152  gSystem->Load( "libHijing1_383.so");
153  gSystem->Load( "libKinematics.so");
154  gSystem->Load( "xgeometry.so" );
155 
156  gSystem->Load("libHepMC2_06_09.so");
157  gSystem->Load("libPythia8_1_86.so");
158  gSystem->Load("libPhotos3_61.so");
159  gSystem->Load("libTauola1_1_5.so");
160  gSystem->Load("libEvtGen1_06_00.so");
161 
162  // force gstar load/call
163  gSystem->Load( "gstar.so" );
164  command("call gstar");
165 
166  // Setup RNG seed and map all ROOT TRandom here
167  StarRandom::seed( rngSeed );
169 
170  char rootname[100],fzname[100];
171  sprintf(rootname,"st_hijingevtgen_%d.starsim.root",Index);
172  sprintf(fzname,"gfile o st_hijingevtgen_%d.starsim.fzd",Index);
173 
174  //
175  // Create the primary event generator and insert it
176  // before the geant maker
177  //
178  _primary = new StarPrimaryMaker();
179  {
180  _primary -> SetFileName(rootname);
181  chain -> AddBefore( "geant", _primary );
182  }
183 
184  //
185  // These should be adjusted to your best vertex estimates
186  //
187  _primary -> SetVertex( vx,vy,vz );
188  _primary -> SetSigma( vx_sig,vy_sig,vz_sig );
189 
190  //
191  // Setup an event generator
192  //
193  Hijing();
194 
195  //
196  // Setup single particle
197  //
198  Kinematics();
199 
200  //
201  // Setup decay manager
202  //
203  StarDecayManager *decayMgr = AgUDecay::Manager();
204  /*
205  StarPythia8Decayer *decayPy8 = new StarPythia8Decayer();
206  decayMgr->AddDecayer( 0, decayPy8 ); // Handle any decay requested
207  decayPy8->SetDebug(1);
208  decayPy8->Set("WeakSingleBoson:all = on");
209  */
210  StarEvtGenDecayer *decayEvt = new StarEvtGenDecayer();
211  //decayEvt->SetDecayTable("StRoot/StSimulationMaker/Decay_Table/Jpsi.DEC");
212  decayMgr->AddDecayer( 0, decayEvt ); // Handle any decay requested
213  decayEvt->SetDebug(1);
214 
215  //
216  // Initialize primary event generator and all sub makers
217  //
218  _primary -> Init();
219 
220  //
221  // Setup geometry and set starsim to use agusread for input
222  //
223  //geometry("y2014");
224  command("gkine -4 0");
225  command(fzname);
226 
227  //Double_t pt0 = 3.0;
228  //ptDist = new TF1("ptDist","(x/[0])/(1+(x/[0])^2)^6",0.0,10.0);
229  //ptDist->SetParameter(0, pt0);
230  //ptDist->Draw();
231  ptDist = new TF1("ptDist","[0]*x*TMath::Exp(-x/[1])",minPt,maxPt); //dN/pT/dpT is exp
232  ptDist->SetParameters(1.,1.);//slope = 1.;
233  //yDist = new TF1("yDist","-TMath::Erf(x+2.6)*TMath::Erf(x-2.6)",minY,maxY);
234  yDist = new TF1("yDist","pol0",minY,maxY);
235  yDist->SetParameter(0,1.);
236  //phi, default 0 ~ TMath::TwoPi() flat
237 
238  //
239  // Trigger on nevents
240  //
241  trig( nevents );
242 
243  _primary->event()->Print();
244 
245  command("call agexit"); // Make sure that STARSIM exits properly
246  command("gprint kine");
247 }
248 // ----------------------------------------------------------------------------
249 
void SetFrame(const Char_t *frame, const Double_t val)
void Print(const Option_t *opts="head") const
Star Simple Kinematics Generator.
void SetImpact(Double_t bmin, Double_t bmax)
Set the minimum and maximum impact parameters for the collision (HI generators)
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StChain.cxx:77
static StarParticleData & instance()
Returns a reference to the single instance of this class.
Interface to PDG information.
void SetBlue(const Char_t *b)
Sets the particle species for the blue beam.
Connects VMC to class(es) which handle particle decays.
int & ihpr2(int i)
Definition: Hijing.h:94
Int_t Init()
Initialize generator.
HiParnt_t & hiparnt()
Returns a reference to the hijing parameters.
Definition: StarHijing.h:58
void SetDebug(Int_t dbg=1)
Set the debug level.
float & hipr1(int i)
Definition: Hijing.h:90
Interface to the HIJING event generator.
Definition: StarHijing.h:48
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
STAR wrapper for EvtGen Decayer.
StarGenEvent * event()
Return a pointer to the event.
Main steering class for event generation.
void AddDecayer(Int_t pdgid, TVirtualMCDecayer *decayer)
void SetYell(const Char_t *y)
Sets the particle species for the yellow beam.
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)