| Geant4 Cross Reference |
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 //
27 /// \file eventgenerator/pythia/pythia8decayer/src/Py8Decayer.cc
28 /// \brief Implementation of the Py8Decayer class
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........oooOO0OOooo........oooOO0OOooo......
44
45 Py8Decayer::Py8Decayer() : G4VExtDecayer("Py8Decayer"), fDecayer(nullptr)
46 {
47 // use default path to the xml/configs but do NOT print banner
48 //
49 fDecayer = new Pythia("../share/Pythia8/xmldoc", false);
50
51 fDecayer->setRndmEnginePtr(std::make_shared<Py8DecayerEngine>());
52
53 // this is the trick to make Pythia8 work as "decayer"
54 //
55 fDecayer->readString("ProcessLevel:all = off");
56
57 fDecayer->readString("ProcessLevel:resonanceDecays=on");
58
59 // shut off Pythia8 (default) verbosty
60 //
61 fDecayer->readString("Init:showAllSettings=false");
62 fDecayer->readString("Init:showChangedSettings=false");
63 fDecayer->readString("Init:showAllParticleData=false");
64 fDecayer->readString("Init:showChangedParticleData=false");
65 //
66 // specify how many Py8 events to print out, at either level
67 // in this particular case print out a maximum of 10 events
68 //
69 fDecayer->readString("Next:numberShowProcess = 0");
70 fDecayer->readString("Next:numberShowEvent = 10");
71
72 fDecayer->init();
73
74 // shut off decays of pi0's as we want Geant4 to handle them
75 // if other immediate decay products should be handled by Geant4,
76 // their respective decay modes should be shut off as well
77 //
78 fDecayer->readString("111:onMode = off");
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 Py8Decayer::~Py8Decayer()
84 {
85 delete fDecayer;
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89
90 G4DecayProducts* Py8Decayer::ImportDecayProducts(const G4Track& track)
91 {
92 fDecayer->event.reset();
93
94 G4DecayProducts* dproducts = nullptr;
95
96 G4ParticleDefinition* pd = track.GetDefinition();
97 int pdgid = pd->GetPDGEncoding();
98
99 // check if pdgid is consistent with Pythia8 particle table
100 //
101 if (!fDecayer->particleData.findParticle(pdgid)) {
102 G4cout << " can NOT find pdgid = " << pdgid << " in Pythia8::ParticleData" << G4endl;
103 return dproducts;
104 }
105
106 if (!fDecayer->particleData.canDecay(pdgid)) {
107 G4cout << " Particle of pdgid = " << pdgid << " can NOT be decayed by Pythia8" << G4endl;
108 return dproducts;
109 }
110
111 // NOTE: Energy should be in GeV
112
113 fDecayer->event.append(pdgid, 1, 0, 0, track.GetMomentum().x() / CLHEP::GeV,
114 track.GetMomentum().y() / CLHEP::GeV, track.GetMomentum().z() / CLHEP::GeV,
115 track.GetDynamicParticle()->GetTotalEnergy() / CLHEP::GeV,
116 pd->GetPDGMass() / CLHEP::GeV);
117
118 // specify polarization, if any
119
120 // NOTE: while in Py8 polarization is a double variable ,
121 // in reality it's expected to be -1, 0., or 1 in case of "external" tau's,
122 // similar to LHA SPINUP; see Particle Decays, Hadron and Tau Decays in docs at
123 // https://pythia.org/manuals/pythia8305/Welcome.html
124 // so it's not able to handle anything like 0.99, thus we're rounding off
125 fDecayer->event.back().pol(
126 round(std::cos(track.GetPolarization().angle(track.GetMomentumDirection()))));
127
128 int npart_before_decay = fDecayer->event.size();
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.GetDynamicParticle()));
137
138 // create G4DynamicParticle out of each fDecayer->event entry (except the 1st one)
139 // and push into dproducts
140
141 for (int ip = npart_before_decay; ip < npart_after_decay; ++ip) {
142 // only select final state decay products (direct or via subsequent decays);
143 // skip all others
144 //
145 // NOTE: in general, final state decays products will have
146 // positive status code between 91 and 99
147 // (in case such information could be of interest in the future)
148 //
149 if (fDecayer->event[ip].status() < 0) continue;
150
151 // should we also skip neutrinos ???
152 // if so, skip products with abs(fDecayer->event[ip].id()) of 12, 14, or 16
153
154 G4ParticleDefinition* pddec =
155 G4ParticleTable::GetParticleTable()->FindParticle(fDecayer->event[ip].id());
156 if (!pddec) continue; // maybe we should print out a warning !
157 G4ThreeVector momentum =
158 G4ThreeVector(fDecayer->event[ip].px() * CLHEP::GeV, fDecayer->event[ip].py() * CLHEP::GeV,
159 fDecayer->event[ip].pz() * CLHEP::GeV);
160 dproducts->PushProducts(new G4DynamicParticle(pddec, momentum));
161 }
162
163 return dproducts;
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167