CbmRoot
CbmFastDecayer.cxx
Go to the documentation of this file.
1 #include "FairPrimaryGenerator.h"
2 
3 #include "CbmFastDecayer.h"
4 #include "CbmStack.h"
5 
6 #include "TMCProcess.h"
7 #include "TSystem.h"
8 #include "TVirtualMC.h"
9 #include <TParticle.h>
10 
11 #include <algorithm> // std::find
12 #include <iostream> // std::cout
13 #include <sstream> // stream
14 #include <vector> // std::vector
15 
16 using std::string;
17 using std::stringstream;
18 
19 #include <TVirtualMCDecayer.h>
20 
24  : CbmFastDecayer("fastDecayer", "fastDecayer") {
25  //
26  // Default Construction
27  //
28 }
29 
31 CbmFastDecayer::CbmFastDecayer(const char* fileName, TString particle)
32  : FairGenerator(), fDecayPdgCodes(0), fGeantPdgCodes(0) {
33  //
34  // Standrad Construction
35  //
36  particle += fileName; // dummy against unsed variables
37 }
38 
41  //
42  // Standard Destructor
43  //
44  if (fStack) { delete fStack; }
45  fStack = 0;
46  if (fDecayer) { delete fDecayer; }
47  fDecayer = 0;
48 }
49 
52  //
53  // Standard CbmGenerator Initializer - no input
54  // 1) initialize decayer with default decay and particle table
55  // 2) set the decay mode to force particle
56  // 3) set a user decay table if defined
57  //
58  if (!fDecayer) { Fatal("Init", "CbmFastDecayer has no external decayer!!!"); }
59  // fDecayer = new CbmDecayerEvtGen();
60 
62  // fDecayer->Init(); //read the default decay table DECAY.DEC and particle table
63 
64  //if is set a decay mode: default decay mode is kAll
65  // fDecayer->SetForceDecay(fForceDecay);
66  // fDecayer->ForceDecay();
67 
68  // //if is defined a user decay table
69  // if(fUserDecay)
70  // {
71  // fDecayer->SetDecayTablePath(fUserDecayTablePath);
72  // fDecayer->ReadDecayTable();
73  // }
74  return kTRUE;
75 }
76 
78 Bool_t CbmFastDecayer::ReadEvent(FairPrimaryGenerator* primGen) {
79  //
80  //Generate method - Input none - Output none
81  //For each event:
82  //1)return the stack of the previous generator and select particles to be decayed by EvtGen
83  //2)decay particles selected and put the decay products on the stack
84  //
85  //
86 
87  // ---> Check for primary generator
88  if (!primGen) {
89  Error("ReadEvent", "No PrimaryGenerator!");
90  return kFALSE;
91  }
92 
93  // ---> Check for particle stack
94  fStack = dynamic_cast<CbmStack*>(gMC->GetStack());
95  if (!fStack) {
96  Error("ReadEvent", "Error: No stack found!");
97  return kFALSE;
98  }
99 
100  Float_t polar[3] = {0, 0, 0}; // Polarisation of daughter particles
101  Float_t origin0[3]; // Origin of the parent particle
102  Float_t pc[3],
103  och[3]; // Momentum and origin of the children particles from EvtGen
104  Int_t nt;
105  Float_t tof;
106  Int_t nPrimStack;
107  Int_t nTrackStack;
108  TLorentzVector* mom = new TLorentzVector();
109 
110  // create array with new particles/ deacy products
111  static TClonesArray* particles;
112  if (!particles) particles = new TClonesArray("TParticle", 1000);
113 
114  nPrimStack = fStack->GetNprimary();
115  nTrackStack = fStack->GetNtrack();
116  Info("ReadEvent",
117  "nPrimaries = %d and ntracks = %d before fast decayer",
118  nPrimStack,
119  nTrackStack);
120 
122  primGen->DoTracking(kTRUE);
123  fStack->DoTracking(kTRUE);
124 
126  std::vector<int>::iterator it;
127 
128  // loop over all tracks need to decay
129  for (Int_t iTrack = 0; iTrack < nTrackStack; ++iTrack) {
130 
131  // get particle
132  TParticle* part = fStack->GetParticle(iTrack);
133  Int_t pdg = part->GetPdgCode();
134 
136  it =
137  std::find(fDecayPdgCodes.begin(), fDecayPdgCodes.end(), TMath::Abs(pdg));
138  if (it == fDecayPdgCodes.end()) continue;
139 
140  //check if particle is already decayed
142  //TODO: how?
143  if (/*part->GetStatusCode() != 1 ||*/ part->GetNDaughters() > 0) {
144  Info("ReadEvent",
145  "Attention: particle %d is already decayed (code:%d,daughters:%d)!",
146  pdg,
147  part->GetStatusCode(),
148  part->GetNDaughters());
149  continue;
150  }
151 
153  part->SetStatusCode(11); //Set particle as decayed : change the status code
154 
156  mom->SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
157  Int_t np;
158 
160  do {
161  fDecayer->Decay(pdg, mom);
162  np = fDecayer->ImportParticles(particles); // into array
163  } while (np < 0);
164 
165  Int_t* trackIt = new Int_t[np];
166  Int_t* pParent = new Int_t[np];
167 
168  Info(
169  "ReadEvent", "number of products: np = %d for mother: %d", np - 1, pdg);
170  // mom->Print(); // debug
171 
172  for (int i = 0; i < np; i++) {
173  pParent[i] = -1;
174  trackIt[i] = 0;
175  }
176 
177  //select trackable particles (decay products)
178  if (np > 1) {
179  TParticle* iparticle = 0;
180 
181  for (int i = 1; i < np; i++) {
182  iparticle = (TParticle*) particles->At(i);
183  Int_t ks = iparticle->GetStatusCode();
184  Int_t pdgd = iparticle->GetPdgCode();
185 
186  //track last decay products
187  if (ks == 1) {
188 
190  it = std::find(
191  fGeantPdgCodes.begin(), fGeantPdgCodes.end(), TMath::Abs(pdgd));
192  if (it == fGeantPdgCodes.end())
193  trackIt[i] = 1;
194  else
195  continue;
196  }
197 
198  } //decay particles loop
199 
200  } // if decay products
201 
203  origin0[0] = part->Vx(); //[cm]
204  origin0[1] = part->Vy(); //[cm]
205  origin0[2] = part->Vz(); //[cm]
206 
208  // TParticle* mother = (TParticle *) particles->At(0);
209  // TLorentzVector motherL;
210  // motherL.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy());
211  // motherL.Print();
212 
213  //
214  // Put decay products on the stack
215  //
216  for (int i = 1; i < np; i++) {
217  TParticle* iparticle = (TParticle*) particles->At(i);
218  Int_t kf = iparticle->GetPdgCode();
219  Int_t ksc = iparticle->GetStatusCode();
220  Int_t jpa = iparticle->GetFirstMother()
221  - 1; //jpa = 0 for daughters ofmother particles
222 
226  Int_t iparent = (jpa > 0) ? fStack->GetNtrack() - (i - jpa) : iTrack;
227 
228  Info("ReadEvent",
229  "FirstMother = %d (iparent %d), indicePart = %d, pdg = %d ",
230  jpa,
231  iparent,
232  i,
233  kf);
234 
236  och[0] = origin0[0] + iparticle->Vx(); //[cm]
237  och[1] = origin0[1] + iparticle->Vy(); //[cm]
238  och[2] = origin0[2] + iparticle->Vz(); //[cm]
239  pc[0] = iparticle->Px(); //[GeV/c]
240  pc[1] = iparticle->Py(); //[GeV/c]
241  pc[2] = iparticle->Pz(); //[GeV/c]
242  tof = part->T() + iparticle->T();
243  Double_t ee = iparticle->Energy();
244 
245 
248  Info("ReadEvent",
249  " Add track %d to primary Generator: pdg=%d with parent=%d and do "
250  "tracking %d (stack primaries: %d, tracks: %d)",
251  i,
252  kf,
253  iparent,
254  trackIt[i],
255  fStack->GetNprimary(),
256  fStack->GetNtrack());
257 
258  primGen->AddTrack(kf,
259  pc[0],
260  pc[1],
261  pc[2],
262  och[0],
263  och[1],
264  och[2],
265  iparent - nTrackStack,
266  trackIt[i],
267  ee);
268 
269  } // Particle loop
270 
271  // clean up
272  particles->Clear();
273  if (trackIt) delete[] trackIt;
274  if (pParent) delete[] pParent;
275  }
276 
278  primGen->DoTracking(kFALSE);
279  fStack->DoTracking(kFALSE);
280 
281  return kTRUE;
282 }
283 
284 void CbmFastDecayer::SetParticlesForDecay(char const* pdgs) {
285 
286  stringstream ss(pdgs);
287  int n;
288  char ch;
289 
290  while (ss >> n) {
291  if (ss >> ch)
292  fDecayPdgCodes.push_back(n);
293  else
294  fDecayPdgCodes.push_back(n);
295  }
296 
297  // for (auto i = fDecayPdgCodes.begin(); i != fDecayPdgCodes.end(); ++i)
298  // std::cout << *i << ' ';
299 }
300 
301 void CbmFastDecayer::SetParticlesForGeant(char const* pdgs) {
302 
303  stringstream ss(pdgs);
304  int n;
305  char ch;
306 
307  while (ss >> n) {
308  if (ss >> ch)
309  fGeantPdgCodes.push_back(n);
310  else
311  fGeantPdgCodes.push_back(n);
312  }
313 
314  // for (auto i = fGeantPdgCodes.begin(); i != fGeantPdgCodes.end(); ++i)
315  // std::cout << *i << ' ';
316 }
CbmFastDecayer::SetParticlesForDecay
void SetParticlesForDecay(char const *pdgs="")
Definition: CbmFastDecayer.cxx:284
CbmFastDecayer::fDecayer
TVirtualMCDecayer * fDecayer
pointer to CbmStack
Definition: CbmFastDecayer.h:45
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmFastDecayer
Definition: CbmFastDecayer.h:14
CbmFastDecayer::CbmFastDecayer
CbmFastDecayer()
CbmFastDecayer::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: CbmFastDecayer.cxx:78
CbmFastDecayer::~CbmFastDecayer
~CbmFastDecayer()
Definition: CbmFastDecayer.cxx:40
CbmFastDecayer::fStack
CbmStack * fStack
Definition: CbmFastDecayer.h:44
CbmStack.h
CbmFastDecayer::Init
virtual Bool_t Init()
Definition: CbmFastDecayer.cxx:51
CbmStack::GetNprimary
virtual Int_t GetNprimary() const
Definition: CbmStack.h:185
CbmFastDecayer.h
CbmStack::DoTracking
void DoTracking(Bool_t doTracking=kTRUE)
Definition: CbmStack.h:264
CbmFastDecayer::SetParticlesForGeant
void SetParticlesForGeant(char const *pdgs="")
Definition: CbmFastDecayer.cxx:301
CbmFastDecayer::fGeantPdgCodes
std::vector< int > fGeantPdgCodes
Definition: CbmFastDecayer.h:48
CbmFastDecayer::fDecayPdgCodes
std::vector< int > fDecayPdgCodes
pointer to decayer
Definition: CbmFastDecayer.h:47
CbmStack::GetNtrack
virtual Int_t GetNtrack() const
Definition: CbmStack.h:179
CbmStack::GetParticle
TParticle * GetParticle(Int_t trackId) const
Definition: CbmStack.cxx:420
CbmStack
Definition: CbmStack.h:45
ClassImp
ClassImp(CbmFastDecayer) CbmFastDecayer
Definition: CbmFastDecayer.cxx:21