EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FairPrimaryGenerator.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairPrimaryGenerator.cxx
1 #include "FairPrimaryGenerator.h"
2 #include "FairGenerator.h"
3 #include "FairMCEventHeader.h"
4 #include "FairGenericStack.h"
5 
6 #include "TDatabasePDG.h"
7 #include "TParticlePDG.h"
8 #include "TRandom.h"
9 #include "TMath.h"
10 #include "TF1.h"
11 
12 #include <iostream>
13 
14 Int_t FairPrimaryGenerator::fTotPrim=0;
15 using std::cout;
16 using std::cerr;
17 using std::endl;
18 
19 
20 // ----- Default constructor -------------------------------------------
21 FairPrimaryGenerator::FairPrimaryGenerator()
22  :TNamed(),
23  fBeamX0(0),
24  fBeamY0(0),
25  fBeamSigmaX(0),
26  fBeamSigmaY(0),
27  fTargetZ (new Double_t[1]),
28  fNrTargets(1),
29  fTargetDz(0),
30  fVertex(TVector3(0.,0.,0.)),
31  fNTracks(0),
32  fSmearVertexZ(kFALSE),
33  fSmearGausVertexZ(kFALSE),
34  fSmearVertexXY(kFALSE),
35  fStack(NULL),
36  fGenList(new TObjArray()),
37  fListIter(fGenList->MakeIterator()),
38  fEvent(NULL),
39  fdoTracking(kTRUE),
40  fEventTimeMin(0),
41  fEventTimeMax(0),
42  fEventTime(0),
43  fEventMeanTime(0),
44  fTimeProb(0),
45  fMCIndexOffset(0),
46  fLogger(FairLogger::GetLogger()),
47  fEventNr(0)
48 {
49  fTargetZ[0] = 0.;
50 }
51 // -------------------------------------------------------------------------
52 
53 
54 
55 // ----- Constructor with title and name -------------------------------
56 FairPrimaryGenerator::FairPrimaryGenerator(const char* name, const char* title)
57  :TNamed(name,title),
58  fBeamX0(0),
59  fBeamY0(0),
60  fBeamSigmaX(0),
61  fBeamSigmaY(0),
62  fTargetZ (new Double_t[1]),
63  fNrTargets(1),
64  fTargetDz(0),
65  fVertex(TVector3(0.,0.,0.)),
66  fNTracks(0),
67  fSmearVertexZ(kFALSE),
68  fSmearGausVertexZ(kFALSE),
69  fSmearVertexXY(kFALSE),
70  fStack(NULL),
71  fGenList(new TObjArray()),
72  fListIter(fGenList->MakeIterator()),
73  fEvent(NULL),
74  fdoTracking(kTRUE),
75  fEventTimeMin(0),
76  fEventTimeMax(0),
77  fEventTime(0),
78  fEventMeanTime(0),
79  fTimeProb(NULL),
80  fMCIndexOffset(0),
81  fLogger(FairLogger::GetLogger()),
82  fEventNr(0)
83 {
84  fTargetZ[0] = 0.;
85 }
86 // -------------------------------------------------------------------------
87 Bool_t FairPrimaryGenerator::Init()
88 {
90  for(Int_t i=0; i<fGenList->GetEntries(); i++ ) {
91  FairGenerator* gen= (FairGenerator*) fGenList->At(i);
92  if(gen) { gen->Init(); }
93  }
94  return kTRUE;
95 }
96 
97 
98 // ----- Destructor ----------------------------------------------------
99 FairPrimaryGenerator::~FairPrimaryGenerator()
100 {
101  // cout<<"Enter Destructor of FairPrimaryGenerator"<<endl;
102  // the stack is deleted by FairMCApplication
103  if (1 == fNrTargets) {
104  delete fTargetZ;
105  } else {
106  delete [] fTargetZ;
107  }
108  fGenList->Delete();
109  delete fGenList;
110  delete fListIter;
111 
112  delete fTimeProb;
113  // cout<<"Leave Destructor of FairPrimaryGenerator"<<endl;
114 }
115 // -------------------------------------------------------------------------
116 
117 
118 
119 // ----- Public method GenerateEvent -----------------------------------
120 Bool_t FairPrimaryGenerator::GenerateEvent(FairGenericStack* pStack)
121 {
122  // Check for MCEventHeader
123  if ( ! fEvent) {
124  fLogger->Fatal(MESSAGE_ORIGIN,"FairPrimaryGenerator::GenerateEvent: No MCEventHeader branch!");
125  Fatal("GenerateEvent", "No MCEventHeader branch");
126  }
127 
128  // Initialise
129  fStack = pStack;
130  fNTracks = 0;
131  fEvent->Reset();
132 
133  // Create event vertex
134  MakeVertex();
135  fEvent->SetVertex(fVertex);
136 
137  // Call the ReadEvent methods from all registered generators
138  fListIter->Reset();
139  TObject* obj = 0;
140  FairGenerator* gen = 0;
141  while( (obj = fListIter->Next()) ) {
142  gen = dynamic_cast<FairGenerator*> (obj);
143  if ( ! gen ) { return kFALSE; }
144  const char* genName = gen->GetName();
145  fMCIndexOffset = fNTracks;// number tracks before generator is called
146  Bool_t test = gen->ReadEvent(this);
147  if ( ! test ) {
148  fLogger->Error(MESSAGE_ORIGIN,"FairPrimaryGenerator: ReadEvent failed for generator ", genName );
149  return kFALSE;
150  }
151  }
152 
153  if(fTimeProb!=0) {
154  fEventTime = fEventTime + fTimeProb->GetRandom();
155  } else {
156  fEventTime = fEventTime + gRandom->Uniform( fEventTimeMin, fEventTimeMax);
157  }
158 
159  fEvent->SetTime(fEventTime);
160 
161  fTotPrim += fNTracks;
162  // Screen output
163 
164  // Set the event number if not set already by one of the dedicated generators
165  if (-1 == fEvent->GetEventID()) {
166  fEventNr++;
167  fEvent->SetEventID(fEventNr);
168  }
169 
170  fLogger->Info(MESSAGE_ORIGIN,"FairPrimaryGenerator: (Event %i) %i primary tracks from vertex (%f, %f, %f ) Event Time = %f (ns)" ,fEvent->GetEventID(), fNTracks, fVertex.X(), fVertex.Y(), fVertex.Z(), fEventTime);
171 
172  fEvent->SetNPrim(fNTracks);
173 
174 
175  return kTRUE;
176 }
177 // -------------------------------------------------------------------------
178 
179 
180 
181 // ----- Public method AddTrack ----------------------------------------
182 void FairPrimaryGenerator::AddTrack(Int_t pdgid, Double_t px, Double_t py,
183  Double_t pz, Double_t vx, Double_t vy,
184  Double_t vz, Int_t parent,Bool_t wanttracking,Double_t e)
185 {
186 
187  // ---> Add event vertex to track vertex
188  vx += fVertex.X();
189  vy += fVertex.Y();
190  vz += fVertex.Z();
191 
192  // ---> Convert K0 and AntiK0 into K0s and K0L
193  if ( pdgid == 311 || pdgid == -311 ) {
194  Double_t test = gRandom->Uniform(0.,1.);
195  if (test >= 0.5) { pdgid = 310; } // K0S
196  else { pdgid = 130; } // K0L
197  }
198 
199  // ---> Check whether particle type is in PDG Database
200  TDatabasePDG* pdgBase = TDatabasePDG::Instance();
201  if ( ! pdgBase ) {
202  Fatal("FairPrimaryGenerator",
203  "No TDatabasePDG instantiated");
204  } else {
205  TParticlePDG* pdgPart = pdgBase->GetParticle(pdgid);
206  if ( ! pdgPart ) {
207  if( e<0) {
208  cerr << "\033[5m\033[31m -E FairPrimaryGenerator: PDG code " << pdgid << " not found in database.\033[0m " << endl;
209  cerr << "\033[5m\033[31m -E FairPrimaryGenerator: Discarding particle \033[0m " << endl;
210  cerr << "\033[5m\033[31m -E FairPrimaryGenerator: now MC Index is corrupted \033[0m " << endl;
211  return;
212  } else {
213  cout << "\033[5m\033[31m -W FairPrimaryGenerator: PDG code " << pdgid << " not found in database. This warning can be savely ignored.\033[0m " << endl;
214  }
215  }
216  }
217  // ---> Get mass and calculate energy of particle
218  if(e<0) {
219  Double_t mass = pdgBase->GetParticle(pdgid)->Mass();
220  e = TMath::Sqrt( px*px + py*py + pz*pz + mass*mass );
221  }// else, use the value of e given to the function
222 
223  // ---> Set all other parameters required by PushTrack
224  Int_t doTracking = 0; // Go to tracking
225  if(fdoTracking && wanttracking) { doTracking = 1 ; }
226  Int_t dummyparent = -1; // Primary particle (now the value is -1 by default)
227  Double_t tof = 0.; // Time of flight
228  Double_t polx = 0.; // Polarisation
229  Double_t poly = 0.;
230  Double_t polz = 0.;
231  Int_t ntr = 0; // Track number; to be filled by the stack
232  Double_t weight = 1.; // Weight
233  Int_t status = 0; // Generation status
234 
235  if( parent!=-1) { parent+=fMCIndexOffset; }// correct for tracks which are in list before generator is called
236  // Add track to stack
237  fStack->PushTrack(doTracking, dummyparent, pdgid, px, py, pz, e, vx, vy, vz,
238  tof, polx, poly, polz, kPPrimary, ntr, weight, status, parent);
239  fNTracks++;
240 
241 }
242 // -------------------------------------------------------------------------
243 
244 
245 
246 // ----- Public method SetBeam -----------------------------------------
247 void FairPrimaryGenerator::SetBeam(Double_t x0, Double_t y0,
248  Double_t sigmaX, Double_t sigmaY)
249 {
250  fBeamX0 = x0;
251  fBeamY0 = y0;
252  fBeamSigmaX = sigmaX;
253  fBeamSigmaY = sigmaY;
254 }
255 // -------------------------------------------------------------------------
256 
257 
258 
259 // ----- Public method SetTarget ---------------------------------------
260 void FairPrimaryGenerator::SetTarget(Double_t z, Double_t dz)
261 {
262 
263  fTargetZ[0] = z;
264  fTargetDz = dz;
265 }
266 // -------------------------------------------------------------------------
267 
268 // ----- Public method SetMultTarget -----------------------------------
269 void FairPrimaryGenerator::SetMultTarget(Int_t nroftargets, Double_t* targetpos, Double_t dz)
270 {
271 
272  if (1 == fNrTargets) {
273  delete fTargetZ;
274  } else {
275  delete[] fTargetZ;
276  }
277 
278  fNrTargets = nroftargets;
279 
280  fTargetZ = new Double_t[nroftargets];
281  for (Int_t i=0; i<nroftargets; i++) {
282  fTargetZ[i] = targetpos[i];
283  }
284  fTargetDz = dz;
285 
286 }
287 // -------------------------------------------------------------------------
288 
289 
290 
291 // ----- Private method MakeVertex -------------------------------------
292 void FairPrimaryGenerator::MakeVertex()
293 {
294  Double_t vx = fBeamX0;
295  Double_t vy = fBeamY0;
296  Double_t vz;
297  if (1 == fNrTargets) {
298  vz = fTargetZ[0];
299  } else {
300  Int_t Target = (Int_t)gRandom->Uniform(fNrTargets);
301  vz = fTargetZ[Target];
302  }
303 
304  if (fSmearVertexZ) vz = gRandom->Uniform(vz - fTargetDz/2.,
305  vz + fTargetDz/2.);
306 
307  if (fSmearGausVertexZ) { vz = gRandom->Gaus(vz, fTargetDz); }
308  if (fSmearVertexXY) {
309  if (fBeamSigmaX != 0.) { vx = gRandom->Gaus(fBeamX0, fBeamSigmaX); }
310  if (fBeamSigmaY != 0.) { vy = gRandom->Gaus(fBeamY0, fBeamSigmaY); }
311  }
312 
313  fVertex = TVector3(vx,vy,vz);
314 }
315 // -------------------------------------------------------------------------
316 
317 void FairPrimaryGenerator::SetEventTimeInterval(Double_t min, Double_t max)
318 {
319  fEventTimeMin=min;
320  fEventTimeMax=max;
321 
322 }
323 // -------------------------------------------------------------------------
324 
325 void FairPrimaryGenerator::SetEventMeanTime(Double_t mean)
326 {
327  fEventMeanTime =mean;
328  TString form="(1/";
329  form+= mean;
330  form+=")*exp(-x/";
331  form+=mean;
332  form+=")";
333  fTimeProb= new TF1("TimeProb.", form.Data(), 0., mean*10);
334 
335 }
336 // -------------------------------------------------------------------------
337 
338 
339 void FairPrimaryGenerator::SetEventTime(TF1* timeProb)
340 {
341  if (timeProb!=0) {
342  fTimeProb= timeProb;
343  } else {
344  cout << " \033[5m\033[31m -E FairPrimaryGenerator: invalid time function, Event time is not SET \033[0m " << endl;
345  }
346 }
347 
348 
349 
350 
351 
352 
353 ClassImp(FairPrimaryGenerator)
354