Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Pythia6Decayer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Pythia6Decayer.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4Pythia6Decayer.cc,v 1.2 2014/10/07 03:06:54 mccumber Exp $
27 //
30 
31 // ----------------------------------------------------------------------------
32 // According to TPythia6Decayer class in Root:
33 // http://root.cern.ch/
34 // see http://root.cern.ch/root/License.html
35 // ----------------------------------------------------------------------------
36 
37 #include "G4Pythia6Decayer.hh"
38 #include "Pythia6.hh"
39 
40 #include <Geant4/G4DecayProducts.hh>
41 #include <Geant4/G4DynamicParticle.hh>
42 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
43 #include <Geant4/G4ParticleTable.hh>
44 #include <Geant4/G4String.hh> // for G4String
45 #include <Geant4/G4SystemOfUnits.hh>
46 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
47 #include <Geant4/G4Track.hh>
48 #include <Geant4/G4Types.hh> // for G4int, G4bool, G4double
49 #include <Geant4/G4VExtDecayer.hh> // for G4VExtDecayer
50 
51 #include <CLHEP/Vector/LorentzVector.h>
52 
53 #include <cstdlib> // for abs
54 #include <iostream> // for operator<<, basic_ostream:...
55 #include <string> // for operator<<
56 
57 const EDecayType G4Pythia6Decayer::fgkDefaultDecayType = kAll;
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 G4Pythia6Decayer::G4Pythia6Decayer()
62  : G4VExtDecayer("G4Pythia6Decayer")
63  , fMessenger(this)
64  , fVerboseLevel(0)
65  , fDecayType(fgkDefaultDecayType)
66  , fDecayProductsArray(nullptr)
67 {
69 
70  fDecayProductsArray = new ParticleVector();
71 
72  ForceDecay(fDecayType);
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 
77 G4Pythia6Decayer::~G4Pythia6Decayer()
78 {
80 
81  delete fDecayProductsArray;
82 }
83 
84 //
85 // private methods
86 //
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 
90 G4ParticleDefinition* G4Pythia6Decayer::
91  GetParticleDefinition(const Pythia6Particle* particle, G4bool warn) const
92 {
94 
95  // get particle definition from G4ParticleTable
96  G4int pdgEncoding = particle->fKF;
97  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
98  G4ParticleDefinition* particleDefinition = nullptr;
99  if (pdgEncoding != 0)
100  {
101  particleDefinition = particleTable->FindParticle(pdgEncoding);
102  }
103 
104  if (particleDefinition == nullptr && warn)
105  {
106  std::cerr
107  << "G4Pythia6Decayer: GetParticleDefinition: " << std::endl
108  << "G4ParticleTable::FindParticle() for particle with PDG = "
109  << pdgEncoding
110  << " failed." << std::endl;
111  }
112 
113  return particleDefinition;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 G4DynamicParticle*
119 G4Pythia6Decayer::CreateDynamicParticle(const Pythia6Particle* particle) const
120 {
122 
123  // get particle properties
124  const G4ParticleDefinition* particleDefinition = GetParticleDefinition(particle);
125  if (!particleDefinition)
126  {
127  return nullptr;
128  }
129 
130  G4ThreeVector momentum = GetParticleMomentum(particle);
131 
132  // create G4DynamicParticle
133  G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, momentum);
134 
135  return dynamicParticle;
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
139 
140 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(
141  const Pythia6Particle* particle) const
142 {
144 
145  G4ThreeVector position = G4ThreeVector(particle->fVx * cm,
146  particle->fVy * cm,
147  particle->fVz * cm);
148  return position;
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 
153 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(
154  const Pythia6Particle* particle) const
155 {
157 
158  G4ThreeVector momentum = G4ThreeVector(particle->fPx * GeV,
159  particle->fPy * GeV,
160  particle->fPz * GeV);
161  return momentum;
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165 
166 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
167 {
169 
170  G4int np = 0;
171  for (G4int i = 1; i <= 5; i++)
172  {
173  if (std::abs(Pythia6::Instance()->GetKFDP(channel, i)) == particle)
174  {
175  np++;
176  }
177  }
178  return np;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
184 {
186 
187  Pythia6* pythia6 = Pythia6::Instance();
188 
189  G4int kc = pythia6->Pycomp(particle);
190  pythia6->SetMDCY(kc, 1, 1);
191 
192  G4int ifirst = pythia6->GetMDCY(kc, 2);
193  G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
194 
195  //
196  // Loop over decay channels
197  for (G4int channel = ifirst; channel <= ilast; channel++)
198  {
199  if (CountProducts(channel, product) >= mult)
200  {
201  pythia6->SetMDME(channel, 1, 1);
202  }
203  else
204  {
205  pythia6->SetMDME(channel, 1, 0);
206  }
207  }
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
211 
212 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products,
213  G4int* mult, G4int npart)
214 {
216 
217  Pythia6* pythia6 = Pythia6::Instance();
218 
219  G4int kc = pythia6->Pycomp(particle);
220  pythia6->SetMDCY(kc, 1, 1);
221  G4int ifirst = pythia6->GetMDCY(kc, 2);
222  G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
223  //
224  // Loop over decay channels
225  for (G4int channel = ifirst; channel <= ilast; channel++)
226  {
227  G4int nprod = 0;
228  for (G4int i = 0; i < npart; i++)
229  {
230  nprod += (CountProducts(channel, products[i]) >= mult[i]);
231  }
232  if (nprod)
233  {
234  pythia6->SetMDME(channel, 1, 1);
235  }
236  else
237  {
238  pythia6->SetMDME(channel, 1, 0);
239  }
240  }
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244 
245 void G4Pythia6Decayer::ForceHadronicD()
246 {
248 
249  const G4int kNHadrons = 4;
250  G4int channel;
251  G4int hadron[kNHadrons] = {411, 421, 431, 4112};
252 
253  // for D+ -> K0* (-> K- pi+) pi+
254  G4int iKstar0 = 313;
255  G4int iKstarbar0 = -313;
256  G4int iKPlus = 321;
257  G4int iKMinus = -321;
258  G4int iPiPlus = 211;
259  G4int iPiMinus = -211;
260 
261  G4int products[2] = {iKPlus, iPiMinus}, mult[2] = {1, 1};
262  ForceParticleDecay(iKstar0, products, mult, 2);
263 
264  // for Ds -> Phi pi+
265  G4int iPhi = 333;
266  ForceParticleDecay(iPhi, iKPlus, 2); // Phi->K+K-
267 
268  G4int decayP1[kNHadrons][3] = {
269  {iKMinus, iPiPlus, iPiPlus},
270  {iKMinus, iPiPlus, 0},
271  {iKPlus, iKstarbar0, 0},
272  {-1, -1, -1}};
273  G4int decayP2[kNHadrons][3] = {
274  {iKstarbar0, iPiPlus, 0},
275  {-1, -1, -1},
276  {iPhi, iPiPlus, 0},
277  {-1, -1, -1}};
278 
279  Pythia6* pythia6 = Pythia6::Instance();
280  for (G4int ihadron = 0; ihadron < kNHadrons; ihadron++)
281  {
282  G4int kc = pythia6->Pycomp(hadron[ihadron]);
283  pythia6->SetMDCY(kc, 1, 1);
284  G4int ifirst = pythia6->GetMDCY(kc, 2);
285  G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
286 
287  for (channel = ifirst; channel <= ilast; channel++)
288  {
289  if ((pythia6->GetKFDP(channel, 1) == decayP1[ihadron][0] &&
290  pythia6->GetKFDP(channel, 2) == decayP1[ihadron][1] &&
291  pythia6->GetKFDP(channel, 3) == decayP1[ihadron][2] &&
292  pythia6->GetKFDP(channel, 4) == 0) ||
293  (pythia6->GetKFDP(channel, 1) == decayP2[ihadron][0] &&
294  pythia6->GetKFDP(channel, 2) == decayP2[ihadron][1] &&
295  pythia6->GetKFDP(channel, 3) == decayP2[ihadron][2] &&
296  pythia6->GetKFDP(channel, 4) == 0))
297  {
298  pythia6->SetMDME(channel, 1, 1);
299  }
300  else
301  {
302  pythia6->SetMDME(channel, 1, 0);
303  } // selected channel ?
304  } // decay channels
305  } // hadrons
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 
310 void G4Pythia6Decayer::ForceOmega()
311 {
313 
314  Pythia6* pythia6 = Pythia6::Instance();
315 
316  G4int iLambda0 = 3122;
317  G4int iKMinus = -321;
318 
319  G4int kc = pythia6->Pycomp(3334);
320  pythia6->SetMDCY(kc, 1, 1);
321  G4int ifirst = pythia6->GetMDCY(kc, 2);
322  G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
323 
324  for (G4int channel = ifirst; channel <= ilast; channel++)
325  {
326  if (pythia6->GetKFDP(channel, 1) == iLambda0 &&
327  pythia6->GetKFDP(channel, 2) == iKMinus &&
328  pythia6->GetKFDP(channel, 3) == 0)
329  {
330  pythia6->SetMDME(channel, 1, 1);
331  }
332  else
333  {
334  pythia6->SetMDME(channel, 1, 0);
335  }
336  // selected channel ?
337  } // decay channels
338 }
339 
340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341 
342 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
343 {
345 
346  Pythia6::Instance()->SetMSTJ(21, 2);
347 
348  if (fDecayType == kNoDecayHeavy)
349  {
350  return;
351  }
352 
353  //
354  // select mode
355  G4int products[3];
356  G4int mult[3];
357 
358  switch (decayType)
359  {
360  case kHardMuons:
361  products[0] = 13;
362  products[1] = 443;
363  products[2] = 100443;
364  mult[0] = 1;
365  mult[1] = 1;
366  mult[2] = 1;
367  ForceParticleDecay(511, products, mult, 3);
368  ForceParticleDecay(521, products, mult, 3);
369  ForceParticleDecay(531, products, mult, 3);
370  ForceParticleDecay(5122, products, mult, 3);
371  ForceParticleDecay(5132, products, mult, 3);
372  ForceParticleDecay(5232, products, mult, 3);
373  ForceParticleDecay(5332, products, mult, 3);
374  ForceParticleDecay(100443, 443, 1); // Psi' -> J/Psi X
375  ForceParticleDecay(443, 13, 2); // J/Psi -> mu+ mu-
376 
377  ForceParticleDecay(411, 13, 1); // D+/-
378  ForceParticleDecay(421, 13, 1); // D0
379  ForceParticleDecay(431, 13, 1); // D_s
380  ForceParticleDecay(4122, 13, 1); // Lambda_c
381  ForceParticleDecay(4132, 13, 1); // Xsi_c
382  ForceParticleDecay(4232, 13, 1); // Sigma_c
383  ForceParticleDecay(4332, 13, 1); // Omega_c
384  break;
385 
386  case kSemiMuonic:
387  ForceParticleDecay(411, 13, 1); // D+/-
388  ForceParticleDecay(421, 13, 1); // D0
389  ForceParticleDecay(431, 13, 1); // D_s
390  ForceParticleDecay(4122, 13, 1); // Lambda_c
391  ForceParticleDecay(4132, 13, 1); // Xsi_c
392  ForceParticleDecay(4232, 13, 1); // Sigma_c
393  ForceParticleDecay(4332, 13, 1); // Omega_c
394  ForceParticleDecay(511, 13, 1); // B0
395  ForceParticleDecay(521, 13, 1); // B+/-
396  ForceParticleDecay(531, 13, 1); // B_s
397  ForceParticleDecay(5122, 13, 1); // Lambda_b
398  ForceParticleDecay(5132, 13, 1); // Xsi_b
399  ForceParticleDecay(5232, 13, 1); // Sigma_b
400  ForceParticleDecay(5332, 13, 1); // Omega_b
401  break;
402 
403  case kDiMuon:
404  ForceParticleDecay(113, 13, 2); // rho
405  ForceParticleDecay(221, 13, 2); // eta
406  ForceParticleDecay(223, 13, 2); // omega
407  ForceParticleDecay(333, 13, 2); // phi
408  ForceParticleDecay(443, 13, 2); // J/Psi
409  ForceParticleDecay(100443, 13, 2); // Psi'
410  ForceParticleDecay(553, 13, 2); // Upsilon
411  ForceParticleDecay(100553, 13, 2); // Upsilon'
412  ForceParticleDecay(200553, 13, 2); // Upsilon''
413  break;
414 
415  case kSemiElectronic:
416  ForceParticleDecay(411, 11, 1); // D+/-
417  ForceParticleDecay(421, 11, 1); // D0
418  ForceParticleDecay(431, 11, 1); // D_s
419  ForceParticleDecay(4122, 11, 1); // Lambda_c
420  ForceParticleDecay(4132, 11, 1); // Xsi_c
421  ForceParticleDecay(4232, 11, 1); // Sigma_c
422  ForceParticleDecay(4332, 11, 1); // Omega_c
423  ForceParticleDecay(511, 11, 1); // B0
424  ForceParticleDecay(521, 11, 1); // B+/-
425  ForceParticleDecay(531, 11, 1); // B_s
426  ForceParticleDecay(5122, 11, 1); // Lambda_b
427  ForceParticleDecay(5132, 11, 1); // Xsi_b
428  ForceParticleDecay(5232, 11, 1); // Sigma_b
429  ForceParticleDecay(5332, 11, 1); // Omega_b
430  break;
431 
432  case kDiElectron:
433  ForceParticleDecay(113, 11, 2); // rho
434  ForceParticleDecay(333, 11, 2); // phi
435  ForceParticleDecay(221, 11, 2); // eta
436  ForceParticleDecay(223, 11, 2); // omega
437  ForceParticleDecay(443, 11, 2); // J/Psi
438  ForceParticleDecay(100443, 11, 2); // Psi'
439  ForceParticleDecay(553, 11, 2); // Upsilon
440  ForceParticleDecay(100553, 11, 2); // Upsilon'
441  ForceParticleDecay(200553, 11, 2); // Upsilon''
442  break;
443 
444  case kBJpsiDiMuon:
445 
446  products[0] = 443;
447  products[1] = 100443;
448  mult[0] = 1;
449  mult[1] = 1;
450 
451  ForceParticleDecay(511, products, mult, 2); // B0 -> J/Psi (Psi') X
452  ForceParticleDecay(521, products, mult, 2); // B+/- -> J/Psi (Psi') X
453  ForceParticleDecay(531, products, mult, 2); // B_s -> J/Psi (Psi') X
454  ForceParticleDecay(5122, products, mult, 2); // Lambda_b -> J/Psi (Psi')X
455  ForceParticleDecay(100443, 443, 1); // Psi' -> J/Psi X
456  ForceParticleDecay(443, 13, 2); // J/Psi -> mu+ mu-
457  break;
458 
459  case kBPsiPrimeDiMuon:
460  ForceParticleDecay(511, 100443, 1); // B0
461  ForceParticleDecay(521, 100443, 1); // B+/-
462  ForceParticleDecay(531, 100443, 1); // B_s
463  ForceParticleDecay(5122, 100443, 1); // Lambda_b
464  ForceParticleDecay(100443, 13, 2); // Psi'
465  break;
466 
467  case kBJpsiDiElectron:
468  ForceParticleDecay(511, 443, 1); // B0
469  ForceParticleDecay(521, 443, 1); // B+/-
470  ForceParticleDecay(531, 443, 1); // B_s
471  ForceParticleDecay(5122, 443, 1); // Lambda_b
472  ForceParticleDecay(443, 11, 2); // J/Psi
473  break;
474 
475  case kBJpsi:
476  ForceParticleDecay(511, 443, 1); // B0
477  ForceParticleDecay(521, 443, 1); // B+/-
478  ForceParticleDecay(531, 443, 1); // B_s
479  ForceParticleDecay(5122, 443, 1); // Lambda_b
480  break;
481 
482  case kBPsiPrimeDiElectron:
483  ForceParticleDecay(511, 100443, 1); // B0
484  ForceParticleDecay(521, 100443, 1); // B+/-
485  ForceParticleDecay(531, 100443, 1); // B_s
486  ForceParticleDecay(5122, 100443, 1); // Lambda_b
487  ForceParticleDecay(100443, 11, 2); // Psi'
488  break;
489 
490  case kPiToMu:
491  ForceParticleDecay(211, 13, 1); // pi->mu
492  break;
493 
494  case kKaToMu:
495  ForceParticleDecay(321, 13, 1); // K->mu
496  break;
497 
498  case kWToMuon:
499  ForceParticleDecay(24, 13, 1); // W -> mu
500  break;
501 
502  case kWToCharm:
503  ForceParticleDecay(24, 4, 1); // W -> c
504  break;
505 
506  case kWToCharmToMuon:
507  ForceParticleDecay(24, 4, 1); // W -> c
508  ForceParticleDecay(411, 13, 1); // D+/- -> mu
509  ForceParticleDecay(421, 13, 1); // D0 -> mu
510  ForceParticleDecay(431, 13, 1); // D_s -> mu
511  ForceParticleDecay(4122, 13, 1); // Lambda_c
512  ForceParticleDecay(4132, 13, 1); // Xsi_c
513  ForceParticleDecay(4232, 13, 1); // Sigma_c
514  ForceParticleDecay(4332, 13, 1); // Omega_c
515  break;
516 
517  case kZDiMuon:
518  ForceParticleDecay(23, 13, 2); // Z -> mu+ mu-
519  break;
520 
521  case kHadronicD:
522  ForceHadronicD();
523  break;
524 
525  case kPhiKK:
526  ForceParticleDecay(333, 321, 2); // Phi->K+K-
527  break;
528 
529  case kOmega:
530  ForceOmega();
531 
532  case kAll:
533  break;
534 
535  case kNoDecay:
536  Pythia6::Instance()->SetMSTJ(21, 0);
537  break;
538 
539  case kNoDecayHeavy:
540  break;
541 
542  case kMaxDecay:
543  break;
544  }
545 }
546 
547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
548 
549 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
550 {
552 
553  Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
554 }
555 
556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
557 
558 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
559 {
561 
562  return Pythia6::Instance()->ImportParticles(particles, "All");
563 }
564 
565 //
566 // public methods
567 //
568 
569 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
570 
571 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
572 {
574 
575  // get particle momentum
576  G4ThreeVector momentum = track.GetMomentum();
577  G4double etot = track.GetDynamicParticle()->GetTotalEnergy();
578  ;
579  CLHEP::HepLorentzVector p;
580  p[0] = momentum.x() / GeV;
581  p[1] = momentum.y() / GeV;
582  p[2] = momentum.z() / GeV;
583  p[3] = etot / GeV;
584 
585  // get particle PDG
586  // ask G4Pythia6Decayer to get PDG encoding
587  // (in order to get PDG from extended TDatabasePDG
588  // in case the standard PDG code is not defined)
589  G4ParticleDefinition* particleDef = track.GetDefinition();
590  G4int pdgEncoding = particleDef->GetPDGEncoding();
591 
592  // let Pythia6Decayer decay the particle
593  // and import the decay products
594  Decay(pdgEncoding, p);
595  G4int nofParticles = ImportParticles(fDecayProductsArray);
596 
597  if (fVerboseLevel > 0)
598  {
599  std::cout << "nofParticles: " << nofParticles << std::endl;
600  }
601 
602  // convert decay products Pythia6Particle type
603  // to G4DecayProducts
604  G4DecayProducts* decayProducts = new G4DecayProducts(*(track.GetDynamicParticle()));
605 
606  G4int counter = 0;
607  for (G4int i = 0; i < nofParticles; i++)
608  {
609  // get particle from ParticleVector
610  Pythia6Particle* particle = (*fDecayProductsArray)[i];
611 
612  G4int status = particle->fKS;
613  G4int pdg = particle->fKF;
614  if (status > 0 && status < 11 &&
615  std::abs(pdg) != 12 && std::abs(pdg) != 14 && std::abs(pdg) != 16)
616  {
617  // pass to tracking final particles only;
618  // skip neutrinos
619 
620  if (fVerboseLevel > 0)
621  {
622  std::cout << " " << i << "th particle PDG: " << pdg << " ";
623  }
624 
625  // create G4DynamicParticle
626  G4DynamicParticle* dynamicParticle = CreateDynamicParticle(particle);
627 
628  if (dynamicParticle)
629  {
630  if (fVerboseLevel > 0)
631  {
632  std::cout << " G4 particle name: "
633  << dynamicParticle->GetDefinition()->GetParticleName()
634  << std::endl;
635  }
636 
637  // add dynamicParticle to decayProducts
638  decayProducts->PushProducts(dynamicParticle);
639 
640  counter++;
641  }
642  }
643  }
644  if (fVerboseLevel > 0)
645  {
646  std::cout << "nofParticles for tracking: " << counter << std::endl;
647  }
648 
649  return decayProducts;
650 }
651 
652 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
653 
654 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
655 {
657 
658  // Do nothing if the decay type is not different from current one
659  if (decayType == fDecayType)
660  {
661  return;
662  }
663 
664  fDecayType = decayType;
665  ForceDecay(fDecayType);
666 }