LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Py8Decayer.cc
Go to the documentation of this file.
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 //
31 
32 #include "Py8Decayer.hh"
33 
34 #include "Geant4/G4DynamicParticle.hh"
35 #include "Geant4/G4ParticleTable.hh"
36 #include "Geant4/G4Track.hh"
37 #include "Geant4/G4SystemOfUnits.hh"
38 
39 #include <cstdlib>
40 #include <sstream>
41 
42 using namespace Pythia8;
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
47  : G4VExtDecayer("Py8Decayer"),
48  fDecayer(nullptr)
49 {
50 
51  // use default path to the xml/configs but do NOT print banner
52  //
53  fDecayer = new Pythia( "../share/Pythia8/xmldoc", false );
54 
55  // this is the trick to make Pythia8 work as "decayer"
56  //
57  fDecayer->readString("ProcessLevel:all = off");
58 
59  fDecayer->readString("ProcessLevel:resonanceDecays=on");
60 
61  // shut off Pythia8 (default) verbosty
62  //
63  fDecayer->readString("Init:showAllSettings=false");
64  fDecayer->readString("Init:showChangedSettings=false");
65  fDecayer->readString("Init:showAllParticleData=false");
66  fDecayer->readString("Init:showChangedParticleData=false");
67  //
68  // specify how many Py8 events to print out, at either level
69  // in this particular case print out a maximum of 10 events
70  //
71 
72  fDecayer->readString( "Next:numberShowProcess = 0" );
73 
74  //fDecayer->readString( "Next:numberShowEvent = 10" );
75  ostringstream s;
76 
77  int nShowEvent = 10;
78  if (const char* env_nShowEvent = std::getenv("PY8_NSHOWEVENT")) {
79  nShowEvent = std::atoi(env_nShowEvent);
80  G4cout << "=====> set pythia8 Next:numberShowEvent to "
81  << nShowEvent << " (" << env_nShowEvent << ")" << G4endl;
82  }
83  //s.str("");
84  //s.clear();
85  s << "Next:numberShowEvent = " << nShowEvent;
86  //G4cout << "=====> set fDecayer->readString(" << s.str().c_str() << ")" << G4endl;
87  fDecayer->readString( s.str().c_str() );
88 
89 
90  fDecayer->init();
91 
92  // shut off decays of pi0's as we want Geant4 to handle them
93  // if other immediate decay products should be handled by Geant4,
94  // their respective decay modes should be shut off as well
95  //
96  fDecayer->readString("111:onMode = off");
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
102 {
103 
104  delete fDecayer;
105 
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 
110 G4DecayProducts* Py8Decayer::ImportDecayProducts(const G4Track& track)
111 {
112 
113  fDecayer->event.reset();
114 
115  G4DecayProducts* dproducts = nullptr;
116 
117  G4ParticleDefinition* pd = track.GetDefinition();
118  int pdgid = pd->GetPDGEncoding();
119 
120  // check if pdgid is consistent with Pythia8 particle table
121  //
122  if ( !fDecayer->particleData.findParticle( pdgid ) )
123  {
124  G4cout << " can NOT find pdgid = " << pdgid
125  << " in Pythia8::ParticleData" << G4endl;
126  return dproducts;
127  }
128 
129  if ( !fDecayer->particleData.canDecay(pdgid) )
130  {
131  G4cout << " Particle of pdgid = " << pdgid
132  << " can NOT be decayed by Pythia8" << G4endl;
133  return dproducts;
134  }
135 
136  // NOTE: Energy should be in GeV
137 
138  fDecayer->event.append( pdgid, 1, 0, 0,
139  track.GetMomentum().x() / CLHEP::GeV,
140  track.GetMomentum().y() / CLHEP::GeV,
141  track.GetMomentum().z() / CLHEP::GeV,
142  track.GetDynamicParticle()->GetTotalEnergy() / CLHEP::GeV,
143  pd->GetPDGMass() / CLHEP::GeV );
144 
145  // specify polarization, if any
146 
147  // NOTE: while in Py8 polarization is a double variable ,
148  // in reality it's expected to be -1, 0., or 1 in case of "external" tau's,
149  // similar to LHA SPINUP; see Particle Decays, Hadron and Tau Decays in docs at
150  // https://pythia.org/manuals/pythia8305/Welcome.html
151  // so it's not able to handle anything like 0.99, thus we're rounding off
152  fDecayer->event.back().pol(
153  round( std::cos( track.GetPolarization().angle( track.GetMomentumDirection() ) )
154  )
155  );
156 
157  int npart_before_decay = fDecayer->event.size();
158 
159  fDecayer->next();
160 
161  int npart_after_decay = fDecayer->event.size();
162 
163  // create & fill up decay products
164  //
165  dproducts = new G4DecayProducts(*(track.GetDynamicParticle()));
166 
167  // create G4DynamicParticle out of each fDecayer->event entry (except the 1st one)
168  // and push into dproducts
169 
170  for ( int ip=npart_before_decay; ip<npart_after_decay; ++ip )
171  {
172 
173  // only select final state decay products (direct or via subsequent decays);
174  // skip all others
175  //
176  // NOTE: in general, final state decays products will have
177  // positive status code between 91 and 99
178  // (in case such information could be of interest in the future)
179  //
180  if ( fDecayer->event[ip].status() < 0 ) continue;
181 
182  // should we also skip neutrinos ???
183  // if so, skip products with abs(fDecayer->event[ip].id()) of 12, 14, or 16
184 
185  G4ParticleDefinition* pddec =
186  G4ParticleTable::GetParticleTable()->FindParticle( fDecayer->event[ip].id() );
187  if ( !pddec ) continue; // maybe we should print out a warning !
188  G4ThreeVector momentum = G4ThreeVector( fDecayer->event[ip].px() * CLHEP::GeV,
189  fDecayer->event[ip].py() * CLHEP::GeV,
190  fDecayer->event[ip].pz() * CLHEP::GeV );
191  dproducts->PushProducts( new G4DynamicParticle( pddec, momentum) );
192  }
193 
194  return dproducts;
195 
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual G4DecayProducts * ImportDecayProducts(const G4Track &)
Definition: Py8Decayer.cc:110
virtual ~Py8Decayer()
Definition: Py8Decayer.cc:101
Float_t track
Definition: plot.C:35
Pythia8::Pythia * fDecayer
Definition: Py8Decayer.hh:58