| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \file eventgenerator/pythia/decayer6/src/G
28 /// \brief Implementation of the G4Pythia6Deca
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
38 #include "Pythia6.hh"
39
40 #include "G4DecayProducts.hh"
41 #include "G4DecayTable.hh"
42 #include "G4DynamicParticle.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Track.hh"
46
47 #include <CLHEP/Vector/LorentzVector.h>
48 #include <cmath>
49
50 const EDecayType G4Pythia6Decayer::fgkDefaultD
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 G4Pythia6Decayer::G4Pythia6Decayer()
55 : G4VExtDecayer("G4Pythia6Decayer"),
56 fMessenger(this),
57 fVerboseLevel(0),
58 fDecayType(fgkDefaultDecayType),
59 fDecayProductsArray(0)
60 {
61 /// Standard constructor
62
63 fDecayProductsArray = new ParticleVector();
64
65 ForceDecay(fDecayType);
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 G4Pythia6Decayer::~G4Pythia6Decayer()
71 {
72 /// Destructor
73
74 delete fDecayProductsArray;
75 }
76
77 //
78 // private methods
79 //
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 G4ParticleDefinition* G4Pythia6Decayer::GetPar
84
85 {
86 /// Return G4 particle definition for given
87
88 // get particle definition from G4ParticleTa
89 G4int pdgEncoding = particle->fKF;
90 G4ParticleTable* particleTable = G4ParticleT
91 G4ParticleDefinition* particleDefinition = 0
92 if (pdgEncoding != 0) particleDefinition = p
93
94 if (particleDefinition == 0 && warn) {
95 G4cerr << "G4Pythia6Decayer: GetParticleDe
96 << "G4ParticleTable::FindParticle()
97 << " failed." << std::endl;
98 }
99
100 return particleDefinition;
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oo
104
105 G4DynamicParticle* G4Pythia6Decayer::CreateDyn
106 {
107 /// Create G4DynamicParticle.
108
109 // get particle properties
110 const G4ParticleDefinition* particleDefiniti
111 if (!particleDefinition) return 0;
112
113 G4ThreeVector momentum = GetParticleMomentum
114
115 // create G4DynamicParticle
116 G4DynamicParticle* dynamicParticle = new G4D
117
118 return dynamicParticle;
119 }
120
121 //....oooOO0OOooo........oooOO0OOooo........oo
122
123 G4ThreeVector G4Pythia6Decayer::GetParticlePos
124 {
125 /// Return particle vertex position.
126
127 G4ThreeVector position =
128 G4ThreeVector(particle->fVx * cm, particle
129 return position;
130 }
131
132 //....oooOO0OOooo........oooOO0OOooo........oo
133
134 G4ThreeVector G4Pythia6Decayer::GetParticleMom
135 {
136 /// Return particle momentum.
137
138 G4ThreeVector momentum =
139 G4ThreeVector(particle->fPx * GeV, particl
140 return momentum;
141 }
142
143 //....oooOO0OOooo........oooOO0OOooo........oo
144
145 G4int G4Pythia6Decayer::CountProducts(G4int ch
146 {
147 /// Count number of decay products
148
149 G4int np = 0;
150 for (G4int i = 1; i <= 5; i++)
151 if (std::abs(Pythia6::Instance()->GetKFDP(
152 return np;
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oo
156
157 void G4Pythia6Decayer::ForceParticleDecay(G4in
158 {
159 /// Force decay of particle into products wi
160
161 Pythia6* pythia6 = Pythia6::Instance();
162
163 G4int kc = pythia6->Pycomp(particle);
164 pythia6->SetMDCY(kc, 1, 1);
165
166 G4int ifirst = pythia6->GetMDCY(kc, 2);
167 G4int ilast = ifirst + pythia6->GetMDCY(kc,
168
169 //
170 // Loop over decay channels
171 for (G4int channel = ifirst; channel <= ilas
172 if (CountProducts(channel, product) >= mul
173 pythia6->SetMDME(channel, 1, 1);
174 }
175 else {
176 pythia6->SetMDME(channel, 1, 0);
177 }
178 }
179 }
180
181 //....oooOO0OOooo........oooOO0OOooo........oo
182
183 void G4Pythia6Decayer::ForceParticleDecay(G4in
184 {
185 /// Force decay of particle into products wi
186
187 Pythia6* pythia6 = Pythia6::Instance();
188
189 G4int kc = pythia6->Pycomp(particle);
190 pythia6->SetMDCY(kc, 1, 1);
191 G4int ifirst = pythia6->GetMDCY(kc, 2);
192 G4int ilast = ifirst + pythia6->GetMDCY(kc,
193 //
194 // Loop over decay channels
195 for (G4int channel = ifirst; channel <= ilas
196 G4int nprod = 0;
197 for (G4int i = 0; i < npart; i++)
198 nprod += (CountProducts(channel, product
199 if (nprod)
200 pythia6->SetMDME(channel, 1, 1);
201 else {
202 pythia6->SetMDME(channel, 1, 0);
203 }
204 }
205 }
206
207 //....oooOO0OOooo........oooOO0OOooo........oo
208
209 void G4Pythia6Decayer::ForceHadronicD()
210 {
211 /// Force golden D decay modes
212
213 const G4int kNHadrons = 4;
214 G4int channel;
215 G4int hadron[kNHadrons] = {411, 421, 431, 41
216
217 // for D+ -> K0* (-> K- pi+) pi+
218 G4int iKstar0 = 313;
219 G4int iKstarbar0 = -313;
220 G4int iKPlus = 321;
221 G4int iKMinus = -321;
222 G4int iPiPlus = 211;
223 G4int iPiMinus = -211;
224
225 G4int products[2] = {iKPlus, iPiMinus}, mult
226 ForceParticleDecay(iKstar0, products, mult,
227
228 // for Ds -> Phi pi+
229 G4int iPhi = 333;
230 ForceParticleDecay(iPhi, iKPlus, 2); // Phi
231
232 G4int decayP1[kNHadrons][3] = {
233 {iKMinus, iPiPlus, iPiPlus}, {iKMinus, iPi
234 G4int decayP2[kNHadrons][3] = {
235 {iKstarbar0, iPiPlus, 0}, {-1, -1, -1}, {i
236
237 Pythia6* pythia6 = Pythia6::Instance();
238 for (G4int ihadron = 0; ihadron < kNHadrons;
239 G4int kc = pythia6->Pycomp(hadron[ihadron]
240 pythia6->SetMDCY(kc, 1, 1);
241 G4int ifirst = pythia6->GetMDCY(kc, 2);
242 G4int ilast = ifirst + pythia6->GetMDCY(kc
243
244 for (channel = ifirst; channel <= ilast; c
245 if ((pythia6->GetKFDP(channel, 1) == dec
246 && pythia6->GetKFDP(channel, 2) ==
247 && pythia6->GetKFDP(channel, 3) ==
248 && pythia6->GetKFDP(channel, 4) ==
249 || (pythia6->GetKFDP(channel, 1) ==
250 && pythia6->GetKFDP(channel, 2)
251 && pythia6->GetKFDP(channel, 3)
252 && pythia6->GetKFDP(channel, 4)
253 {
254 pythia6->SetMDME(channel, 1, 1);
255 }
256 else {
257 pythia6->SetMDME(channel, 1, 0);
258 } // selected channel ?
259 } // decay channels
260 } // hadrons
261 }
262
263 //....oooOO0OOooo........oooOO0OOooo........oo
264
265 void G4Pythia6Decayer::ForceOmega()
266 {
267 /// Force Omega -> Lambda K- Decay
268
269 Pythia6* pythia6 = Pythia6::Instance();
270
271 G4int iLambda0 = 3122;
272 G4int iKMinus = -321;
273
274 G4int kc = pythia6->Pycomp(3334);
275 pythia6->SetMDCY(kc, 1, 1);
276 G4int ifirst = pythia6->GetMDCY(kc, 2);
277 G4int ilast = ifirst + pythia6->GetMDCY(kc,
278
279 for (G4int channel = ifirst; channel <= ilas
280 if (pythia6->GetKFDP(channel, 1) == iLambd
281 && pythia6->GetKFDP(channel, 3) == 0)
282 pythia6->SetMDME(channel, 1, 1);
283 else
284 pythia6->SetMDME(channel, 1, 0);
285 // selected channel ?
286 } // decay channels
287 }
288
289 //....oooOO0OOooo........oooOO0OOooo........oo
290
291 void G4Pythia6Decayer::ForceDecay(EDecayType d
292 {
293 /// Force a particle decay mode
294
295 Pythia6::Instance()->SetMSTJ(21, 2);
296
297 if (fDecayType == kNoDecayHeavy) return;
298
299 //
300 // select mode
301 G4int products[3];
302 G4int mult[3];
303
304 switch (decayType) {
305 case kHardMuons:
306 products[0] = 13;
307 products[1] = 443;
308 products[2] = 100443;
309 mult[0] = 1;
310 mult[1] = 1;
311 mult[2] = 1;
312 ForceParticleDecay(511, products, mult,
313 ForceParticleDecay(521, products, mult,
314 ForceParticleDecay(531, products, mult,
315 ForceParticleDecay(5122, products, mult,
316 ForceParticleDecay(5132, products, mult,
317 ForceParticleDecay(5232, products, mult,
318 ForceParticleDecay(5332, products, mult,
319 ForceParticleDecay(100443, 443, 1); //
320 ForceParticleDecay(443, 13, 2); // J/Ps
321
322 ForceParticleDecay(411, 13, 1); // D+/-
323 ForceParticleDecay(421, 13, 1); // D0
324 ForceParticleDecay(431, 13, 1); // D_s
325 ForceParticleDecay(4122, 13, 1); // Lam
326 ForceParticleDecay(4132, 13, 1); // Xsi
327 ForceParticleDecay(4232, 13, 1); // Sig
328 ForceParticleDecay(4332, 13, 1); // Ome
329 break;
330
331 case kSemiMuonic:
332 ForceParticleDecay(411, 13, 1); // D+/-
333 ForceParticleDecay(421, 13, 1); // D0
334 ForceParticleDecay(431, 13, 1); // D_s
335 ForceParticleDecay(4122, 13, 1); // Lam
336 ForceParticleDecay(4132, 13, 1); // Xsi
337 ForceParticleDecay(4232, 13, 1); // Sig
338 ForceParticleDecay(4332, 13, 1); // Ome
339 ForceParticleDecay(511, 13, 1); // B0
340 ForceParticleDecay(521, 13, 1); // B+/-
341 ForceParticleDecay(531, 13, 1); // B_s
342 ForceParticleDecay(5122, 13, 1); // Lam
343 ForceParticleDecay(5132, 13, 1); // Xsi
344 ForceParticleDecay(5232, 13, 1); // Sig
345 ForceParticleDecay(5332, 13, 1); // Ome
346 break;
347
348 case kDiMuon:
349 ForceParticleDecay(113, 13, 2); // rho
350 ForceParticleDecay(221, 13, 2); // eta
351 ForceParticleDecay(223, 13, 2); // omeg
352 ForceParticleDecay(333, 13, 2); // phi
353 ForceParticleDecay(443, 13, 2); // J/Ps
354 ForceParticleDecay(100443, 13, 2); // P
355 ForceParticleDecay(553, 13, 2); // Upsi
356 ForceParticleDecay(100553, 13, 2); // U
357 ForceParticleDecay(200553, 13, 2); // U
358 break;
359
360 case kSemiElectronic:
361 ForceParticleDecay(411, 11, 1); // D+/-
362 ForceParticleDecay(421, 11, 1); // D0
363 ForceParticleDecay(431, 11, 1); // D_s
364 ForceParticleDecay(4122, 11, 1); // Lam
365 ForceParticleDecay(4132, 11, 1); // Xsi
366 ForceParticleDecay(4232, 11, 1); // Sig
367 ForceParticleDecay(4332, 11, 1); // Ome
368 ForceParticleDecay(511, 11, 1); // B0
369 ForceParticleDecay(521, 11, 1); // B+/-
370 ForceParticleDecay(531, 11, 1); // B_s
371 ForceParticleDecay(5122, 11, 1); // Lam
372 ForceParticleDecay(5132, 11, 1); // Xsi
373 ForceParticleDecay(5232, 11, 1); // Sig
374 ForceParticleDecay(5332, 11, 1); // Ome
375 break;
376
377 case kDiElectron:
378 ForceParticleDecay(113, 11, 2); // rho
379 ForceParticleDecay(333, 11, 2); // phi
380 ForceParticleDecay(221, 11, 2); // eta
381 ForceParticleDecay(223, 11, 2); // omeg
382 ForceParticleDecay(443, 11, 2); // J/Ps
383 ForceParticleDecay(100443, 11, 2); // P
384 ForceParticleDecay(553, 11, 2); // Upsi
385 ForceParticleDecay(100553, 11, 2); // U
386 ForceParticleDecay(200553, 11, 2); // U
387 break;
388
389 case kBJpsiDiMuon:
390
391 products[0] = 443;
392 products[1] = 100443;
393 mult[0] = 1;
394 mult[1] = 1;
395
396 ForceParticleDecay(511, products, mult,
397 ForceParticleDecay(521, products, mult,
398 ForceParticleDecay(531, products, mult,
399 ForceParticleDecay(5122, products, mult,
400 ForceParticleDecay(100443, 443, 1); //
401 ForceParticleDecay(443, 13, 2); // J/Ps
402 break;
403
404 case kBPsiPrimeDiMuon:
405 ForceParticleDecay(511, 100443, 1); //
406 ForceParticleDecay(521, 100443, 1); //
407 ForceParticleDecay(531, 100443, 1); //
408 ForceParticleDecay(5122, 100443, 1); //
409 ForceParticleDecay(100443, 13, 2); // P
410 break;
411
412 case kBJpsiDiElectron:
413 ForceParticleDecay(511, 443, 1); // B0
414 ForceParticleDecay(521, 443, 1); // B+/
415 ForceParticleDecay(531, 443, 1); // B_s
416 ForceParticleDecay(5122, 443, 1); // La
417 ForceParticleDecay(443, 11, 2); // J/Ps
418 break;
419
420 case kBJpsi:
421 ForceParticleDecay(511, 443, 1); // B0
422 ForceParticleDecay(521, 443, 1); // B+/
423 ForceParticleDecay(531, 443, 1); // B_s
424 ForceParticleDecay(5122, 443, 1); // La
425 break;
426
427 case kBPsiPrimeDiElectron:
428 ForceParticleDecay(511, 100443, 1); //
429 ForceParticleDecay(521, 100443, 1); //
430 ForceParticleDecay(531, 100443, 1); //
431 ForceParticleDecay(5122, 100443, 1); //
432 ForceParticleDecay(100443, 11, 2); // P
433 break;
434
435 case kPiToMu:
436 ForceParticleDecay(211, 13, 1); // pi->
437 break;
438
439 case kKaToMu:
440 ForceParticleDecay(321, 13, 1); // K->m
441 break;
442
443 case kWToMuon:
444 ForceParticleDecay(24, 13, 1); // W ->
445 break;
446
447 case kWToCharm:
448 ForceParticleDecay(24, 4, 1); // W -> c
449 break;
450
451 case kWToCharmToMuon:
452 ForceParticleDecay(24, 4, 1); // W -> c
453 ForceParticleDecay(411, 13, 1); // D+/-
454 ForceParticleDecay(421, 13, 1); // D0
455 ForceParticleDecay(431, 13, 1); // D_s
456 ForceParticleDecay(4122, 13, 1); // Lam
457 ForceParticleDecay(4132, 13, 1); // Xsi
458 ForceParticleDecay(4232, 13, 1); // Sig
459 ForceParticleDecay(4332, 13, 1); // Ome
460 break;
461
462 case kZDiMuon:
463 ForceParticleDecay(23, 13, 2); // Z ->
464 break;
465
466 case kHadronicD:
467 ForceHadronicD();
468 break;
469
470 case kPhiKK:
471 ForceParticleDecay(333, 321, 2); // Phi
472 break;
473
474 case kOmega:
475 ForceOmega();
476
477 case kAll:
478 break;
479
480 case kNoDecay:
481 Pythia6::Instance()->SetMSTJ(21, 0);
482 break;
483
484 case kNoDecayHeavy:
485 break;
486
487 case kMaxDecay:
488 break;
489 }
490 }
491
492 //....oooOO0OOooo........oooOO0OOooo........oo
493
494 void G4Pythia6Decayer::Decay(G4int pdg, const
495 {
496 /// Decay a particle of type IDPART (PDG cod
497
498 Pythia6::Instance()->Py1ent(0, pdg, p.e(), p
499 }
500
501 //....oooOO0OOooo........oooOO0OOooo........oo
502
503 G4int G4Pythia6Decayer::ImportParticles(Partic
504 {
505 /// Get the decay products into the passed P
506
507 return Pythia6::Instance()->ImportParticles(
508 }
509
510 //
511 // public methods
512 //
513
514 //....oooOO0OOooo........oooOO0OOooo........oo
515
516 G4DecayProducts* G4Pythia6Decayer::ImportDecay
517 {
518 /// Import decay products
519
520 // get particle momentum
521 G4ThreeVector momentum = track.GetMomentum()
522 G4double etot = track.GetDynamicParticle()->
523 ;
524 CLHEP::HepLorentzVector p;
525 p[0] = momentum.x() / GeV;
526 p[1] = momentum.y() / GeV;
527 p[2] = momentum.z() / GeV;
528 p[3] = etot / GeV;
529
530 // get particle PDG
531 // ask G4Pythia6Decayer to get PDG encoding
532 // (in order to get PDG from extended TDatab
533 // in case the standard PDG code is not defi
534 G4ParticleDefinition* particleDef = track.Ge
535 G4int pdgEncoding = particleDef->GetPDGEncod
536
537 // let Pythia6Decayer decay the particle
538 // and import the decay products
539 Decay(pdgEncoding, p);
540 G4int nofParticles = ImportParticles(fDecayP
541
542 if (fVerboseLevel > 0) {
543 G4cout << "nofParticles: " << nofParticles
544 }
545
546 // convert decay products Pythia6Particle ty
547 // to G4DecayProducts
548 G4DecayProducts* decayProducts = new G4Decay
549
550 G4int counter = 0;
551 for (G4int i = 0; i < nofParticles; i++) {
552 // get particle from ParticleVector
553 Pythia6Particle* particle = (*fDecayProduc
554
555 G4int status = particle->fKS;
556 G4int pdg = particle->fKF;
557 if (status > 0 && status < 11 && std::abs(
558 && std::abs(pdg) != 16)
559 {
560 // pass to tracking final particles only
561 // skip neutrinos
562
563 if (fVerboseLevel > 0) {
564 G4cout << " " << i << "th particle PD
565 }
566
567 // create G4DynamicParticle
568 G4DynamicParticle* dynamicParticle = Cre
569
570 if (dynamicParticle) {
571 if (fVerboseLevel > 0) {
572 G4cout << " G4 particle name: " <<
573 << G4endl;
574 }
575
576 // add dynamicParticle to decayProduct
577 decayProducts->PushProducts(dynamicPar
578
579 counter++;
580 }
581 }
582 }
583 if (fVerboseLevel > 0) {
584 G4cout << "nofParticles for tracking: " <<
585 }
586
587 return decayProducts;
588 }
589
590 //....oooOO0OOooo........oooOO0OOooo........oo
591
592 void G4Pythia6Decayer::ForceDecayType(EDecayTy
593 {
594 /// Force a given decay type
595
596 // Do nothing if the decay type is not diffe
597 if (decayType == fDecayType) return;
598
599 fDecayType = decayType;
600 ForceDecay(fDecayType);
601 }
602