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