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