ECCE @ EIC 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 //
29 
30 // ----------------------------------------------------------------------------
31 // According to TPythia6Decayer class in Root:
32 // http://root.cern.ch/
33 // see http://root.cern.ch/root/License.html
34 // ----------------------------------------------------------------------------
35 
36 #include "G4Pythia6Decayer.hh"
37 #include "Pythia6.hh"
38 
39 #include "G4DynamicParticle.hh"
40 #include "G4DecayProducts.hh"
41 #include "G4DecayTable.hh"
42 #include "G4ParticleTable.hh"
43 #include "G4Track.hh"
44 #include "G4SystemOfUnits.hh"
45 
46 #include <CLHEP/Vector/LorentzVector.h>
47 
48 #include <cmath>
49 
50 const EDecayType G4Pythia6Decayer::fgkDefaultDecayType = kAll;
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 
54 G4Pythia6Decayer::G4Pythia6Decayer()
55  : G4VExtDecayer("G4Pythia6Decayer"),
56  fMessenger(this),
57  fVerboseLevel(0),
58  fDecayType(fgkDefaultDecayType),
59  fDecayProductsArray(0)
60 {
62 
63  fDecayProductsArray = new ParticleVector();
64 
65  ForceDecay(fDecayType);
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
70 G4Pythia6Decayer::~G4Pythia6Decayer()
71 {
73 
74  delete fDecayProductsArray;
75 }
76 
77 //
78 // private methods
79 //
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
83 G4ParticleDefinition* G4Pythia6Decayer::
84 GetParticleDefinition(const Pythia6Particle* particle, G4bool warn) const
85 {
87 
88  // get particle definition from G4ParticleTable
89  G4int pdgEncoding = particle->fKF;
90  G4ParticleTable* particleTable
91  = G4ParticleTable::GetParticleTable();
92  G4ParticleDefinition* particleDefinition = 0;
93  if (pdgEncoding != 0)
94  particleDefinition = particleTable->FindParticle(pdgEncoding);
95 
96  if ( particleDefinition == 0 && warn) {
97  G4cerr
98  << "G4Pythia6Decayer: GetParticleDefinition: " << std::endl
99  << "G4ParticleTable::FindParticle() for particle with PDG = "
100  << pdgEncoding
101  << " failed." << std::endl;
102  }
103 
104  return particleDefinition;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 G4DynamicParticle*
110 G4Pythia6Decayer::CreateDynamicParticle(const Pythia6Particle* particle) const
111 {
113 
114  // get particle properties
115  const G4ParticleDefinition* particleDefinition
116  = GetParticleDefinition(particle);
117  if ( ! particleDefinition ) return 0;
118 
119  G4ThreeVector momentum = GetParticleMomentum(particle);
120 
121  // create G4DynamicParticle
122  G4DynamicParticle* dynamicParticle
123  = new G4DynamicParticle(particleDefinition, momentum);
124 
125  return dynamicParticle;
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 
130 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(
131  const Pythia6Particle* particle) const
132 {
134 
135  G4ThreeVector position
136  = G4ThreeVector(particle->fVx * cm,
137  particle->fVy * cm,
138  particle->fVz * cm);
139  return position;
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(
145  const Pythia6Particle* particle) const
146 {
148 
149  G4ThreeVector momentum
150  = G4ThreeVector(particle->fPx * GeV,
151  particle->fPy * GeV,
152  particle->fPz * GeV);
153  return momentum;
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157 
158 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
159 {
161 
162  G4int np = 0;
163  for ( G4int i=1; i<=5; i++ )
164  if ( std::abs(Pythia6::Instance()->GetKFDP(channel,i) ) == particle )
165  np++;
166  return np;
167 }
168 
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170 
171 void
172 G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
173 {
175 
176  Pythia6* pythia6 = Pythia6::Instance();
177 
178  G4int kc = pythia6->Pycomp(particle);
179  pythia6->SetMDCY(kc,1,1);
180 
181  G4int ifirst = pythia6->GetMDCY(kc,2);
182  G4int ilast = ifirst + pythia6->GetMDCY(kc,3)-1;
183 
184  //
185  // Loop over decay channels
186  for (G4int channel= ifirst; channel <= ilast; channel++) {
187  if (CountProducts(channel,product) >= mult) {
188  pythia6->SetMDME(channel,1,1);
189  } else {
190  pythia6->SetMDME(channel,1,0);
191  }
192  }
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 
197 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products,
198  G4int* mult, G4int npart)
199 {
201 
202  Pythia6* pythia6 = Pythia6::Instance();
203 
204  G4int kc = pythia6->Pycomp(particle);
205  pythia6->SetMDCY(kc,1,1);
206  G4int ifirst = pythia6->GetMDCY(kc,2);
207  G4int ilast = ifirst+pythia6->GetMDCY(kc,3)-1;
208  //
209  // Loop over decay channels
210  for (G4int channel = ifirst; channel <= ilast; channel++) {
211  G4int nprod = 0;
212  for (G4int i = 0; i < npart; i++)
213  nprod += (CountProducts(channel, products[i]) >= mult[i]);
214  if (nprod)
215  pythia6->SetMDME(channel,1,1);
216  else {
217  pythia6->SetMDME(channel,1,0);
218  }
219  }
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
223 
224 void G4Pythia6Decayer::ForceHadronicD()
225 {
227 
228  const G4int kNHadrons = 4;
229  G4int channel;
230  G4int hadron[kNHadrons] = {411, 421, 431, 4112};
231 
232  // for D+ -> K0* (-> K- pi+) pi+
233  G4int iKstar0 = 313;
234  G4int iKstarbar0 = -313;
235  G4int iKPlus = 321;
236  G4int iKMinus = -321;
237  G4int iPiPlus = 211;
238  G4int iPiMinus = -211;
239 
240  G4int products[2] = {iKPlus, iPiMinus}, mult[2] = {1, 1};
241  ForceParticleDecay(iKstar0, products, mult, 2);
242 
243  // for Ds -> Phi pi+
244  G4int iPhi = 333;
245  ForceParticleDecay(iPhi,iKPlus,2); // Phi->K+K-
246 
247  G4int decayP1[kNHadrons][3] = {
248  {iKMinus, iPiPlus, iPiPlus},
249  {iKMinus, iPiPlus, 0 },
250  {iKPlus , iKstarbar0, 0 },
251  {-1 , -1 , -1 }
252  };
253  G4int decayP2[kNHadrons][3] = {
254  {iKstarbar0, iPiPlus, 0 },
255  {-1 , -1 , -1 },
256  {iPhi , iPiPlus, 0 },
257  {-1 , -1 , -1 }
258  };
259 
260  Pythia6* pythia6 = Pythia6::Instance();
261  for ( G4int ihadron = 0; ihadron < kNHadrons; ihadron++ ) {
262  G4int kc = pythia6->Pycomp(hadron[ihadron]);
263  pythia6->SetMDCY(kc,1,1);
264  G4int ifirst = pythia6->GetMDCY(kc,2);
265  G4int ilast = ifirst + pythia6->GetMDCY(kc,3)-1;
266 
267  for (channel = ifirst; channel <= ilast; channel++) {
268  if ((pythia6->GetKFDP(channel,1) == decayP1[ihadron][0] &&
269  pythia6->GetKFDP(channel,2) == decayP1[ihadron][1] &&
270  pythia6->GetKFDP(channel,3) == decayP1[ihadron][2] &&
271  pythia6->GetKFDP(channel,4) == 0) ||
272  (pythia6->GetKFDP(channel,1) == decayP2[ihadron][0] &&
273  pythia6->GetKFDP(channel,2) == decayP2[ihadron][1] &&
274  pythia6->GetKFDP(channel,3) == decayP2[ihadron][2] &&
275  pythia6->GetKFDP(channel,4) == 0)) {
276  pythia6->SetMDME(channel,1,1);
277  } else {
278  pythia6->SetMDME(channel,1,0);
279  } // selected channel ?
280  } // decay channels
281  } // hadrons
282 }
283 
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285 
286 void G4Pythia6Decayer::ForceOmega()
287 {
289 
290  Pythia6* pythia6 = Pythia6::Instance();
291 
292  G4int iLambda0 = 3122;
293  G4int iKMinus = -321;
294 
295  G4int kc = pythia6->Pycomp(3334);
296  pythia6->SetMDCY(kc,1,1);
297  G4int ifirst = pythia6->GetMDCY(kc,2);
298  G4int ilast = ifirst + pythia6->GetMDCY(kc,3)-1;
299 
300  for (G4int channel = ifirst; channel <= ilast; channel++) {
301  if (pythia6->GetKFDP(channel,1) == iLambda0 &&
302  pythia6->GetKFDP(channel,2) == iKMinus &&
303  pythia6->GetKFDP(channel,3) == 0)
304  pythia6->SetMDME(channel,1,1);
305  else
306  pythia6->SetMDME(channel,1,0);
307  // selected channel ?
308  } // decay channels
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 
313 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
314 {
316 
317  Pythia6::Instance()->SetMSTJ(21,2);
318 
319  if ( fDecayType == kNoDecayHeavy ) return;
320 
321  //
322  // select mode
323  G4int products[3];
324  G4int mult[3];
325 
326  switch ( decayType ) {
327 
328  case kHardMuons:
329  products[0] = 13;
330  products[1] = 443;
331  products[2] = 100443;
332  mult[0] = 1;
333  mult[1] = 1;
334  mult[2] = 1;
335  ForceParticleDecay( 511, products, mult, 3);
336  ForceParticleDecay( 521, products, mult, 3);
337  ForceParticleDecay( 531, products, mult, 3);
338  ForceParticleDecay( 5122, products, mult, 3);
339  ForceParticleDecay( 5132, products, mult, 3);
340  ForceParticleDecay( 5232, products, mult, 3);
341  ForceParticleDecay( 5332, products, mult, 3);
342  ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
343  ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
344 
345  ForceParticleDecay( 411,13,1); // D+/-
346  ForceParticleDecay( 421,13,1); // D0
347  ForceParticleDecay( 431,13,1); // D_s
348  ForceParticleDecay( 4122,13,1); // Lambda_c
349  ForceParticleDecay( 4132,13,1); // Xsi_c
350  ForceParticleDecay( 4232,13,1); // Sigma_c
351  ForceParticleDecay( 4332,13,1); // Omega_c
352  break;
353 
354  case kSemiMuonic:
355  ForceParticleDecay( 411,13,1); // D+/-
356  ForceParticleDecay( 421,13,1); // D0
357  ForceParticleDecay( 431,13,1); // D_s
358  ForceParticleDecay( 4122,13,1); // Lambda_c
359  ForceParticleDecay( 4132,13,1); // Xsi_c
360  ForceParticleDecay( 4232,13,1); // Sigma_c
361  ForceParticleDecay( 4332,13,1); // Omega_c
362  ForceParticleDecay( 511,13,1); // B0
363  ForceParticleDecay( 521,13,1); // B+/-
364  ForceParticleDecay( 531,13,1); // B_s
365  ForceParticleDecay( 5122,13,1); // Lambda_b
366  ForceParticleDecay( 5132,13,1); // Xsi_b
367  ForceParticleDecay( 5232,13,1); // Sigma_b
368  ForceParticleDecay( 5332,13,1); // Omega_b
369  break;
370 
371  case kDiMuon:
372  ForceParticleDecay( 113,13,2); // rho
373  ForceParticleDecay( 221,13,2); // eta
374  ForceParticleDecay( 223,13,2); // omega
375  ForceParticleDecay( 333,13,2); // phi
376  ForceParticleDecay( 443,13,2); // J/Psi
377  ForceParticleDecay(100443,13,2);// Psi'
378  ForceParticleDecay( 553,13,2); // Upsilon
379  ForceParticleDecay(100553,13,2);// Upsilon'
380  ForceParticleDecay(200553,13,2);// Upsilon''
381  break;
382 
383  case kSemiElectronic:
384  ForceParticleDecay( 411,11,1); // D+/-
385  ForceParticleDecay( 421,11,1); // D0
386  ForceParticleDecay( 431,11,1); // D_s
387  ForceParticleDecay( 4122,11,1); // Lambda_c
388  ForceParticleDecay( 4132,11,1); // Xsi_c
389  ForceParticleDecay( 4232,11,1); // Sigma_c
390  ForceParticleDecay( 4332,11,1); // Omega_c
391  ForceParticleDecay( 511,11,1); // B0
392  ForceParticleDecay( 521,11,1); // B+/-
393  ForceParticleDecay( 531,11,1); // B_s
394  ForceParticleDecay( 5122,11,1); // Lambda_b
395  ForceParticleDecay( 5132,11,1); // Xsi_b
396  ForceParticleDecay( 5232,11,1); // Sigma_b
397  ForceParticleDecay( 5332,11,1); // Omega_b
398  break;
399 
400  case kDiElectron:
401  ForceParticleDecay( 113,11,2); // rho
402  ForceParticleDecay( 333,11,2); // phi
403  ForceParticleDecay( 221,11,2); // eta
404  ForceParticleDecay( 223,11,2); // omega
405  ForceParticleDecay( 443,11,2); // J/Psi
406  ForceParticleDecay(100443,11,2);// Psi'
407  ForceParticleDecay( 553,11,2); // Upsilon
408  ForceParticleDecay(100553,11,2);// Upsilon'
409  ForceParticleDecay(200553,11,2);// Upsilon''
410  break;
411 
412  case kBJpsiDiMuon:
413 
414  products[0] = 443;
415  products[1] = 100443;
416  mult[0] = 1;
417  mult[1] = 1;
418 
419  ForceParticleDecay( 511, products, mult, 2); // B0 -> J/Psi (Psi') X
420  ForceParticleDecay( 521, products, mult, 2); // B+/- -> J/Psi (Psi') X
421  ForceParticleDecay( 531, products, mult, 2); // B_s -> J/Psi (Psi') X
422  ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi')X
423  ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
424  ForceParticleDecay( 443,13,2); // J/Psi -> mu+ mu-
425  break;
426 
427  case kBPsiPrimeDiMuon:
428  ForceParticleDecay( 511,100443,1); // B0
429  ForceParticleDecay( 521,100443,1); // B+/-
430  ForceParticleDecay( 531,100443,1); // B_s
431  ForceParticleDecay( 5122,100443,1); // Lambda_b
432  ForceParticleDecay(100443,13,2); // Psi'
433  break;
434 
435  case kBJpsiDiElectron:
436  ForceParticleDecay( 511,443,1); // B0
437  ForceParticleDecay( 521,443,1); // B+/-
438  ForceParticleDecay( 531,443,1); // B_s
439  ForceParticleDecay( 5122,443,1); // Lambda_b
440  ForceParticleDecay( 443,11,2); // J/Psi
441  break;
442 
443  case kBJpsi:
444  ForceParticleDecay( 511,443,1); // B0
445  ForceParticleDecay( 521,443,1); // B+/-
446  ForceParticleDecay( 531,443,1); // B_s
447  ForceParticleDecay( 5122,443,1); // Lambda_b
448  break;
449 
450  case kBPsiPrimeDiElectron:
451  ForceParticleDecay( 511,100443,1); // B0
452  ForceParticleDecay( 521,100443,1); // B+/-
453  ForceParticleDecay( 531,100443,1); // B_s
454  ForceParticleDecay( 5122,100443,1); // Lambda_b
455  ForceParticleDecay(100443,11,2); // Psi'
456  break;
457 
458  case kPiToMu:
459  ForceParticleDecay(211,13,1); // pi->mu
460  break;
461 
462  case kKaToMu:
463  ForceParticleDecay(321,13,1); // K->mu
464  break;
465 
466  case kWToMuon:
467  ForceParticleDecay( 24, 13,1); // W -> mu
468  break;
469 
470  case kWToCharm:
471  ForceParticleDecay( 24, 4,1); // W -> c
472  break;
473 
474  case kWToCharmToMuon:
475  ForceParticleDecay( 24, 4,1); // W -> c
476  ForceParticleDecay( 411,13,1); // D+/- -> mu
477  ForceParticleDecay( 421,13,1); // D0 -> mu
478  ForceParticleDecay( 431,13,1); // D_s -> mu
479  ForceParticleDecay( 4122,13,1); // Lambda_c
480  ForceParticleDecay( 4132,13,1); // Xsi_c
481  ForceParticleDecay( 4232,13,1); // Sigma_c
482  ForceParticleDecay( 4332,13,1); // Omega_c
483  break;
484 
485  case kZDiMuon:
486  ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
487  break;
488 
489  case kHadronicD:
490  ForceHadronicD();
491  break;
492 
493  case kPhiKK:
494  ForceParticleDecay(333,321,2); // Phi->K+K-
495  break;
496 
497  case kOmega:
498  ForceOmega();
499 
500  case kAll:
501  break;
502 
503  case kNoDecay:
504  Pythia6::Instance()->SetMSTJ(21,0);
505  break;
506 
507  case kNoDecayHeavy: break;
508 
509  case kMaxDecay: break;
510  }
511 }
512 
513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
514 
515 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
516 {
518 
519  Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
520 }
521 
522 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
523 
524 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
525 {
527 
528  return Pythia6::Instance()->ImportParticles(particles,"All");
529 }
530 
531 //
532 // public methods
533 //
534 
535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
536 
537 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
538 {
540 
541  // get particle momentum
542  G4ThreeVector momentum = track.GetMomentum();
543  G4double etot = track.GetDynamicParticle()->GetTotalEnergy();;
544  CLHEP::HepLorentzVector p;
545  p[0] = momentum.x() / GeV;
546  p[1] = momentum.y() / GeV;
547  p[2] = momentum.z() / GeV;
548  p[3] = etot / GeV;
549 
550  // get particle PDG
551  // ask G4Pythia6Decayer to get PDG encoding
552  // (in order to get PDG from extended TDatabasePDG
553  // in case the standard PDG code is not defined)
554  G4ParticleDefinition* particleDef = track.GetDefinition();
555  G4int pdgEncoding = particleDef->GetPDGEncoding();
556 
557  // let Pythia6Decayer decay the particle
558  // and import the decay products
559  Decay(pdgEncoding, p);
560  G4int nofParticles = ImportParticles(fDecayProductsArray);
561 
562  if ( fVerboseLevel > 0 ) {
563  G4cout << "nofParticles: " << nofParticles << G4endl;
564  }
565 
566  // convert decay products Pythia6Particle type
567  // to G4DecayProducts
568  G4DecayProducts* decayProducts
569  = new G4DecayProducts(*(track.GetDynamicParticle()));
570 
571  G4int counter = 0;
572  for (G4int i=0; i<nofParticles; i++) {
573 
574  // get particle from ParticleVector
575  Pythia6Particle* particle = (*fDecayProductsArray)[i];
576 
577  G4int status = particle->fKS;
578  G4int pdg = particle->fKF;
579  if ( status>0 && status<11 &&
580  std::abs(pdg)!=12 && std::abs(pdg)!=14 && std::abs(pdg)!=16 ) {
581  // pass to tracking final particles only;
582  // skip neutrinos
583 
584  if ( fVerboseLevel > 0 ) {
585  G4cout << " " << i << "th particle PDG: " << pdg << " ";
586  }
587 
588  // create G4DynamicParticle
589  G4DynamicParticle* dynamicParticle
590  = CreateDynamicParticle(particle);
591 
592  if (dynamicParticle) {
593 
594  if ( fVerboseLevel > 0 ) {
595  G4cout << " G4 particle name: "
596  << dynamicParticle->GetDefinition()->GetParticleName()
597  << G4endl;
598  }
599 
600  // add dynamicParticle to decayProducts
601  decayProducts->PushProducts(dynamicParticle);
602 
603  counter++;
604  }
605  }
606  }
607  if ( fVerboseLevel > 0 ) {
608  G4cout << "nofParticles for tracking: " << counter << G4endl;
609  }
610 
611  return decayProducts;
612 }
613 
614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
615 
616 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
617 {
619 
620  // Do nothing if the decay type is not different from current one
621  if ( decayType == fDecayType ) return;
622 
623  fDecayType = decayType;
624  ForceDecay(fDecayType);
625 }