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