py8decayer/src/Py8Decayer.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file eventgenerator/pythia/pythia8decayer/src/Py8Decayer.cc
0028 /// \brief Implementation of the Py8Decayer class
0029 ///
0030 /// \author J. Yarba; FNAL
0031
0032 #include "Py8Decayer.hh"
0033
0034 #include "Py8DecayerEngine.hh"
0035
0036 #include "G4DynamicParticle.hh"
0037 #include "G4ParticleTable.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4Track.hh"
0040
0041 using namespace Pythia8;
0042
0043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0044
0045 Py8Decayer::Py8Decayer() : G4VExtDecayer("Py8Decayer"), fDecayer(nullptr)
0046 {
0047   // use default path to the xml/configs but do NOT print banner
0048   //
0049   fDecayer = new Pythia("../share/Pythia8/xmldoc", false);
0050
0051   fDecayer->setRndmEnginePtr(std::make_shared<Py8DecayerEngine>());
0052
0053   // this is the trick to make Pythia8 work as "decayer"
0054   //
0055   fDecayer->readString("ProcessLevel:all = off");
0056
0057   fDecayer->readString("ProcessLevel:resonanceDecays=on");
0058
0059   // shut off Pythia8 (default) verbosty
0060   //
0061   fDecayer->readString("Init:showAllSettings=false");
0062   fDecayer->readString("Init:showChangedSettings=false");
0063   fDecayer->readString("Init:showAllParticleData=false");
0064   fDecayer->readString("Init:showChangedParticleData=false");
0065   //
0066   // specify how many Py8 events to print out, at either level
0067   // in this particular case print out a maximum of 10 events
0068   //
0069   fDecayer->readString("Next:numberShowProcess = 0");
0070   fDecayer->readString("Next:numberShowEvent = 10");
0071
0072   fDecayer->init();
0073
0074   // shut off decays of pi0's as we want Geant4 to handle them
0075   // if other immediate decay products should be handled by Geant4,
0076   // their respective decay modes should be shut off as well
0077   //
0078   fDecayer->readString("111:onMode = off");
0079 }
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 Py8Decayer::~Py8Decayer()
0084 {
0085   delete fDecayer;
0086 }
0087
0088 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0089
0090 G4DecayProducts* Py8Decayer::ImportDecayProducts(const G4Track& track)
0091 {
0092   fDecayer->event.reset();
0093
0094   G4DecayProducts* dproducts = nullptr;
0095
0096   G4ParticleDefinition* pd = track.GetDefinition();
0097   int pdgid = pd->GetPDGEncoding();
0098
0099   // check if pdgid is consistent with Pythia8 particle table
0100   //
0101   if (!fDecayer->particleData.findParticle(pdgid)) {
0102     G4cout << " can NOT find pdgid = " << pdgid << " in Pythia8::ParticleData" << G4endl;
0103     return dproducts;
0104   }
0105
0106   if (!fDecayer->particleData.canDecay(pdgid)) {
0107     G4cout << " Particle of pdgid = " << pdgid << " can NOT be decayed by Pythia8" << G4endl;
0108     return dproducts;
0109   }
0110
0111   // NOTE: Energy should be in GeV
0112
0113   fDecayer->event.append(pdgid, 1, 0, 0, track.GetMomentum().x() / CLHEP::GeV,
0114                          track.GetMomentum().y() / CLHEP::GeV, track.GetMomentum().z() / CLHEP::GeV,
0115                          track.GetDynamicParticle()->GetTotalEnergy() / CLHEP::GeV,
0116                          pd->GetPDGMass() / CLHEP::GeV);
0117
0118   // specify polarization, if any
0119
0120   // NOTE: while in Py8 polarization is a double variable ,
0121   //       in reality it's expected to be -1, 0., or 1 in case of "external" tau's,
0122   //       similar to LHA SPINUP; see Particle Decays, Hadron and Tau Decays in docs at
0123   //       https://pythia.org/manuals/pythia8305/Welcome.html
0124   //       so it's not able to handle anything like 0.99, thus we're rounding off
0125   fDecayer->event.back().pol(
0126     round(std::cos(track.GetPolarization().angle(track.GetMomentumDirection()))));
0127
0128   int npart_before_decay = fDecayer->event.size();
0129
0130   fDecayer->next();
0131
0132   int npart_after_decay = fDecayer->event.size();
0133
0134   // create & fill up decay products
0135   //
0136   dproducts = new G4DecayProducts(*(track.GetDynamicParticle()));
0137
0138   // create G4DynamicParticle out of each fDecayer->event entry (except the 1st one)
0139   // and push into dproducts
0140
0141   for (int ip = npart_before_decay; ip < npart_after_decay; ++ip) {
0142     // only select final state decay products (direct or via subsequent decays);
0143     // skip all others
0144     //
0145     // NOTE: in general, final state decays products will have
0146     //       positive status code between 91 and 99
0147     //       (in case such information could be of interest in the future)
0148     //
0149     if (fDecayer->event[ip].status() < 0) continue;
0150
0151     // should we also skip neutrinos ???
0152     // if so, skip products with abs(fDecayer->event[ip].id()) of 12, 14, or 16
0153
0154     G4ParticleDefinition* pddec =
0155       G4ParticleTable::GetParticleTable()->FindParticle(fDecayer->event[ip].id());
0156     if (!pddec) continue;  // maybe we should print out a warning !
0157     G4ThreeVector momentum =
0158       G4ThreeVector(fDecayer->event[ip].px() * CLHEP::GeV, fDecayer->event[ip].py() * CLHEP::GeV,
0159                     fDecayer->event[ip].pz() * CLHEP::GeV);
0160     dproducts->PushProducts(new G4DynamicParticle(pddec, momentum));
0161   }
0162
0163   return dproducts;
0164 }
0165
0166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......