| 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/pythia8decayer
28 /// \brief Implementation of the Py8Decayer cl
29 ///
30 /// \author J. Yarba; FNAL
31
32 #include "Py8Decayer.hh"
33
34 #include "Py8DecayerEngine.hh"
35
36 #include "G4DynamicParticle.hh"
37 #include "G4ParticleTable.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4Track.hh"
40
41 using namespace Pythia8;
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 Py8Decayer::Py8Decayer() : G4VExtDecayer("Py8D
46 {
47 // use default path to the xml/configs but d
48 //
49 fDecayer = new Pythia("../share/Pythia8/xmld
50
51 fDecayer->setRndmEnginePtr(std::make_shared<
52
53 // this is the trick to make Pythia8 work as
54 //
55 fDecayer->readString("ProcessLevel:all = off
56
57 fDecayer->readString("ProcessLevel:resonance
58
59 // shut off Pythia8 (default) verbosty
60 //
61 fDecayer->readString("Init:showAllSettings=f
62 fDecayer->readString("Init:showChangedSettin
63 fDecayer->readString("Init:showAllParticleDa
64 fDecayer->readString("Init:showChangedPartic
65 //
66 // specify how many Py8 events to print out,
67 // in this particular case print out a maxim
68 //
69 fDecayer->readString("Next:numberShowProcess
70 fDecayer->readString("Next:numberShowEvent =
71
72 fDecayer->init();
73
74 // shut off decays of pi0's as we want Geant
75 // if other immediate decay products should
76 // their respective decay modes should be sh
77 //
78 fDecayer->readString("111:onMode = off");
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 Py8Decayer::~Py8Decayer()
84 {
85 delete fDecayer;
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oo
89
90 G4DecayProducts* Py8Decayer::ImportDecayProduc
91 {
92 fDecayer->event.reset();
93
94 G4DecayProducts* dproducts = nullptr;
95
96 G4ParticleDefinition* pd = track.GetDefiniti
97 int pdgid = pd->GetPDGEncoding();
98
99 // check if pdgid is consistent with Pythia8
100 //
101 if (!fDecayer->particleData.findParticle(pdg
102 G4cout << " can NOT find pdgid = " << pdgi
103 return dproducts;
104 }
105
106 if (!fDecayer->particleData.canDecay(pdgid))
107 G4cout << " Particle of pdgid = " << pdgid
108 return dproducts;
109 }
110
111 // NOTE: Energy should be in GeV
112
113 fDecayer->event.append(pdgid, 1, 0, 0, track
114 track.GetMomentum().y
115 track.GetDynamicParti
116 pd->GetPDGMass() / CL
117
118 // specify polarization, if any
119
120 // NOTE: while in Py8 polarization is a doub
121 // in reality it's expected to be -1,
122 // similar to LHA SPINUP; see Particle
123 // https://pythia.org/manuals/pythia83
124 // so it's not able to handle anything
125 fDecayer->event.back().pol(
126 round(std::cos(track.GetPolarization().ang
127
128 int npart_before_decay = fDecayer->event.siz
129
130 fDecayer->next();
131
132 int npart_after_decay = fDecayer->event.size
133
134 // create & fill up decay products
135 //
136 dproducts = new G4DecayProducts(*(track.GetD
137
138 // create G4DynamicParticle out of each fDec
139 // and push into dproducts
140
141 for (int ip = npart_before_decay; ip < npart
142 // only select final state decay products
143 // skip all others
144 //
145 // NOTE: in general, final state decays pr
146 // positive status code between 91 a
147 // (in case such information could b
148 //
149 if (fDecayer->event[ip].status() < 0) cont
150
151 // should we also skip neutrinos ???
152 // if so, skip products with abs(fDecayer-
153
154 G4ParticleDefinition* pddec =
155 G4ParticleTable::GetParticleTable()->Fin
156 if (!pddec) continue; // maybe we should
157 G4ThreeVector momentum =
158 G4ThreeVector(fDecayer->event[ip].px() *
159 fDecayer->event[ip].pz() *
160 dproducts->PushProducts(new G4DynamicParti
161 }
162
163 return dproducts;
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oo
167