StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
starsim.decayer.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 class St_geant_Maker;
7 St_geant_Maker *geant_maker = 0;
8 
9 class StarGenEvent;
10 StarGenEvent *event = 0;
11 
12 class StarPrimaryMaker;
13 StarPrimaryMaker *_primary = 0;
14 
15 class StarKinematics;
17 
18 TF1 *ptDist = 0;
19 TF1 *etaDist = 0;
20 
21 // ----------------------------------------------------------------------------
22 void geometry( TString tag, Bool_t agml=true )
23 {
24  TString cmd = "DETP GEOM "; cmd += tag;
25  if ( !geant_maker ) geant_maker = (St_geant_Maker *)chain->GetMaker("geant");
26  geant_maker -> LoadGeometry(cmd);
27  // if ( agml ) command("gexec $STAR_LIB/libxgeometry.so");
28 }
29 // ----------------------------------------------------------------------------
30 void command( TString cmd )
31 {
32  if ( !geant_maker ) geant_maker = (St_geant_Maker *)chain->GetMaker("geant");
33  geant_maker -> Do( cmd );
34 }
35 // ----------------------------------------------------------------------------
36 void trig( Int_t n=1 )
37 {
38  for ( Int_t i=0; i<n; i++ ) {
39 
40  // Clear the chain from the previous event
41  chain->Clear();
42 
43  // Generate 1 mu minus at high pT
44  kinematics->Kine( 1, "D0", 10.0, 100.0, -2.0, 2.0 );
45 
46  // Generate the event
47  chain->Make();
48 
49  // Print the event
50  _primary->event()->Print();
51  command("gprint kine");
52  }
53 }
54 // ----------------------------------------------------------------------------
55 void Kinematics()
56 {
57 
58  // gSystem->Load( "libStarGeneratorPoolPythia6_4_23.so" );
59  gSystem->Load( "libKinematics.so");
60  kinematics = new StarKinematics();
61 
62  _primary->AddGenerator(kinematics);
63 }
64 // ----------------------------------------------------------------------------
65 void starsim( Int_t nevents=1, Int_t rngSeed=1234 )
66 {
67 
68  gROOT->ProcessLine(".L bfc.C");
69  {
70  TString simple = "y2012 geant gstar usexgeom agml ";
71  bfc(0, simple );
72  }
73 
74  gSystem->Load( "libVMC.so");
75 
76  gSystem->Load( "StarGeneratorUtil.so" );
77  gSystem->Load( "StarGeneratorEvent.so" );
78  gSystem->Load( "StarGeneratorBase.so" );
79  gSystem->Load( "StarGeneratorDecay.so" );
80  gSystem->Load( "libPythia8_1_62.so");
81 
82  gSystem->Load( "libMathMore.so" );
83  gSystem->Load( "xgeometry.so" );
84 
85  // Setup RNG seed and map all ROOT TRandom here
86  StarRandom::seed( rngSeed );
88 
89 
90 
91  //
92  // Create the primary event generator and insert it
93  // before the geant maker
94  //
95  // StarPrimaryMaker *
96  _primary = new StarPrimaryMaker();
97  {
98  _primary -> SetFileName( "kinematics.starsim.root");
99  chain -> AddBefore( "geant", _primary );
100  }
101 
102  Kinematics();
103 
104 
105 
106  //
107  // Setup decay manager
108  //
109  StarDecayManager *decayMgr = AgUDecay::Manager();
110  StarPythia8Decayer *decayPy8 = new StarPythia8Decayer();
111  decayMgr->AddDecayer( 0, decayPy8 ); // Handle any decay requested
112  decayPy8->SetDebug(1);
113  decayPy8->Set("WeakSingleBoson:all = on");
114 
115 
116 
117 
118  TString name;
119  double mass, lifetime, charge;
120  int tracktype, pdgcode, g3code;
121 
122  // Particle data
124  // One can add a particle to G3 using...
125  //data.AddParticleToG3( "MyD0", 0.1865E+01, 0.42800E-12, 0., 3, 421, 37, 0, 0 );
126  TParticlePDG* D0 = data.GetParticle("D0");
127  TParticlePDG* rho_pl = data.GetParticle("rho+");
128  TParticlePDG* rho_mn = data.GetParticle("rho-");
129  data.AddParticleToG3( D0, 37 );
130  data.AddParticleToG3( rho_pl, 153 );
131  data.AddParticleToG3( rho_mn, 154 );
132 
133 
134  //
135  // Initialize primary event generator and all sub makers
136  //
137  _primary -> Init();
138 
139  //
140  // Setup geometry and set starsim to use agusread for input
141  //
142  geometry("y2012");
143  command("gkine -4 0");
144  command("gfile o kinematics.starsim.fzd");
145 
146  //
147  // Trigger on nevents
148  //
149  // trig( nevents );
150 
151  // command("call agexit"); // Make sure that STARSIM exits properly
152 
153  trig(1); Validate();
154 
155 }
156 // ----------------------------------------------------------------------------
157 
158 void Validate()
159 {
160  StarGenEvent* event = _primary->event();
161  StarGenParticle* part = (*event)[1];
162  part->Print();
163  double Etotal = part->GetEnergy();
164 
165  // loop over MC particles
166  TTable* gtable = (TTable* ) chain->GetDataSet("g2t_track");
167  assert(gtable);
168  int ntable = gtable->GetNRows();
169 
170  // gtable->Print(0,ntable);
171  return;
172 
173  double Etest = 0;
174  for ( int itable=1; itable<ntable; itable++ )
175  {
176  g2t_track_st* track = (g2t_track_st*)gtable->At(itable);
177  if ( 0 == track->eg_label ) { Etest += track->e; }
178  }
179 
180  cout << "Egener = " << Etotal << endl;
181  cout << "Egeant = " << Etest << endl;
182  cout << "Violation = " << 100*(Etest / Etotal - 1) << "%" << endl;
183 
184 }
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.
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.