Pythia6Interface/plugins/Pythia6Gun.cc

0001 /*
0002  *  \author Julia Yarba
0003  */
0004
0005 #include <iostream>
0006
0007 #include "Pythia6Gun.h"
0008
0009 //#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/LuminosityBlock.h"
0017 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0018
0019 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0020
0021 using namespace edm;
0022 using namespace gen;
0023
0024 Pythia6Gun::Pythia6Gun(const ParameterSet& pset)
0025     : fPy6Service(new Pythia6Service(pset)),
0026       fEvt(nullptr)
0027 // fPDGTable( new DefaultConfig::ParticleDataTable("PDG Table") )
0028 {
0029   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0030
0031   // although there's the method ParameterSet::empty(),
0032   // it looks like it's NOT even necessary to check if it is,
0033   // before trying to extract parameters - if it is empty,
0034   // the default values seem to be taken
0035   //
0036   fMinPhi = pgun_params.getParameter<double>("MinPhi");  // ,-3.14159265358979323846);
0037   fMaxPhi = pgun_params.getParameter<double>("MaxPhi");  // , 3.14159265358979323846);
0038
0039   fHepMCVerbosity = pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false);
0040   fPylistVerbosity = pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
0041   fMaxEventsToPrint = pset.getUntrackedParameter<int>("maxEventsToPrint", 0);
0042
0043   // Turn off banner printout
0044   if (!call_pygive("MSTU(12)=12345")) {
0045     throw edm::Exception(edm::errors::Configuration, "PythiaError") << " pythia did not accept MSTU(12)=12345";
0046   }
0047
0048   usesResource(edm::SharedResourceNames::kPythia6);
0049   produces<HepMCProduct>("unsmeared");
0050 }
0051
0052 Pythia6Gun::~Pythia6Gun() {
0053   if (fPy6Service)
0054     delete fPy6Service;
0055   //
0056   // note that GenEvent or any undelaying (GenVertex, GenParticle) do NOT
0057   // need to be cleaned, as it'll be done automatically by HepMCProduct
0058   //
0059 }
0060
0061 void Pythia6Gun::beginJob() {
0062   // es.getData( fPDGTable ) ;
0063   return;
0064 }
0065
0066 void Pythia6Gun::beginRun(Run const&, EventSetup const& es) {
0067   std::cout << " FYI: MSTU(10)=1 is ENFORCED in Py6-PGuns, for technical reasons" << std::endl;
0068   return;
0069 }
0070
0071 void Pythia6Gun::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const&) {
0072   assert(fPy6Service);
0073
0074   RandomEngineSentry<Pythia6Service> sentry(fPy6Service, lumi.index());
0075
0076   Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0077
0078   fPy6Service->setGeneralParams();
0079   fPy6Service->setCSAParams();
0080   fPy6Service->setSLHAParams();
0081
0082   call_pygive("MSTU(10)=1");
0083
0084   call_pyinit("NONE", "", "", 0.0);
0085 }
0086
0087 void Pythia6Gun::endLuminosityBlock(LuminosityBlock const& lumi, EventSetup const&) {}
0088 void Pythia6Gun::endRun(Run const&, EventSetup const& es) {
0089   // here put in GenRunInfoProduct
0090
0091   call_pystat(1);
0092
0093   return;
0094 }
0095
0096 void Pythia6Gun::attachPy6DecaysToGenEvent() {
0097   for (int iprt = fPartIDs.size(); iprt < pyjets.n; iprt++)  // the pointer is shifted by -1, c++ style
0098   {
0099     int parent = pyjets.k[2][iprt];
0100     if (parent != 0) {
0101       // pull up parent particle
0102       //
0103       HepMC::GenParticle* parentPart = fEvt->barcode_to_particle(parent);
0104       parentPart->set_status(2);  // reset status, to mark that it's decayed
0105
0106       HepMC::GenVertex* DecVtx = new HepMC::GenVertex(
0107           HepMC::FourVector(pyjets.v[0][iprt], pyjets.v[1][iprt], pyjets.v[2][iprt], pyjets.v[3][iprt]));
0108       DecVtx->add_particle_in(parentPart);  // this will cleanup end_vertex if exists,
0109                                             // and replace with the new one
0110                                             // I presume barcode will be given automatically
0111
0112       HepMC::FourVector pmom(pyjets.p[0][iprt], pyjets.p[1][iprt], pyjets.p[2][iprt], pyjets.p[3][iprt]);
0113
0114       int dstatus = 0;
0115       if (pyjets.k[0][iprt] >= 1 && pyjets.k[0][iprt] <= 10) {
0116         dstatus = 1;
0117       } else if (pyjets.k[0][iprt] >= 11 && pyjets.k[0][iprt] <= 20) {
0118         dstatus = 2;
0119       } else if (pyjets.k[0][iprt] >= 21 && pyjets.k[0][iprt] <= 30) {
0120         dstatus = 3;
0121       } else if (pyjets.k[0][iprt] >= 31 && pyjets.k[0][iprt] <= 100) {
0122         dstatus = pyjets.k[0][iprt];
0123       }
0124       HepMC::GenParticle* daughter =
0125           new HepMC::GenParticle(pmom, HepPID::translatePythiatoPDT(pyjets.k[1][iprt]), dstatus);
0126       daughter->suggest_barcode(iprt + 1);
0127       DecVtx->add_particle_out(daughter);
0128       // give particle barcode as well !
0129
0130       int iprt1;
0131       for (iprt1 = iprt + 1; iprt1 < pyjets.n; iprt1++)  // the pointer is shifted by -1, c++ style
0132       {
0133         if (pyjets.k[2][iprt1] != parent)
0134           break;  // another parent particle, break the loop
0135
0136         HepMC::FourVector pmomN(pyjets.p[0][iprt1], pyjets.p[1][iprt1], pyjets.p[2][iprt1], pyjets.p[3][iprt1]);
0137
0138         dstatus = 0;
0139         if (pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10) {
0140           dstatus = 1;
0141         } else if (pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20) {
0142           dstatus = 2;
0143         } else if (pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30) {
0144           dstatus = 3;
0145         } else if (pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100) {
0146           dstatus = pyjets.k[0][iprt1];
0147         }
0148         HepMC::GenParticle* daughterN =
0149             new HepMC::GenParticle(pmomN, HepPID::translatePythiatoPDT(pyjets.k[1][iprt1]), dstatus);
0150         daughterN->suggest_barcode(iprt1 + 1);
0151         DecVtx->add_particle_out(daughterN);
0152       }
0153
0154       iprt = iprt1 - 1;  // reset counter such that it doesn't go over the same child more than once
0155                          // don't forget to offset back into c++ counting, as it's already +1 forward
0156
0157       fEvt->add_vertex(DecVtx);
0158     }
0159   }
0160
0161   return;
0162 }
0163
0164 void Pythia6Gun::produce(edm::Event& evt, const edm::EventSetup&) {
0165   RandomEngineSentry<Pythia6Service> sentry(fPy6Service, evt.streamID());
0166
0167   generateEvent(fPy6Service->randomEngine());
0168
0169   fEvt->set_beam_particles(nullptr, nullptr);
0170   fEvt->set_event_number(evt.id().event());
0171   fEvt->set_signal_process_id(pypars.msti[0]);
0172
0173   attachPy6DecaysToGenEvent();
0174
0175   int evtN = evt.id().event();
0176   if (evtN <= fMaxEventsToPrint) {
0177     if (fPylistVerbosity) {
0178       call_pylist(fPylistVerbosity);
0179     }
0180     if (fHepMCVerbosity) {
0181       if (fEvt)
0182         fEvt->print();
0183     }
0184   }
0185
0186   loadEvent(evt);
0187 }
0188
0189 void Pythia6Gun::loadEvent(edm::Event& evt) {
0190   std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0191
0192   if (fEvt)
0193     bare_product->addHepMCData(fEvt);
0194
0195   evt.put(std::move(bare_product), "unsmeared");
0196
0197   return;
0198 }
0199
0200 HepMC::GenParticle* Pythia6Gun::addAntiParticle(int& ip, int& particleID, double& ee, double& eta, double& phi) {
0201   if (ip < 2)
0202     return nullptr;
0203
0204   // translate PDG to Py6
0205   int py6PID = HepPID::translatePDTtoPythia(particleID);
0206   // Check if particle is its own anti-particle.
0207   int pythiaCode = pycomp_(py6PID);  // this is py6 internal validity check, it takes Pythia6 pid
0208                                      // so actually I'll need to convert
0209   int has_antipart = pydat2.kchg[3 - 1][pythiaCode - 1];
0210   int particleID2 = has_antipart ? -1 * particleID : particleID;  // this is PDG, for HepMC::GenEvent
0211   int py6PID2 = has_antipart ? -1 * py6PID : py6PID;              // this py6 id, for py1ent
0212   double the = 2. * atan(exp(eta));
0213   phi = phi + M_PI;
0214   if (phi > 2. * M_PI) {
0215     phi = phi - 2. * M_PI;
0216   }
0217
0218   // copy over mass of the previous one, because then py6 will pick it up
0219   pyjets.p[4][ip - 1] = pyjets.p[4][ip - 2];
0220
0221   py1ent_(ip, py6PID2, ee, the, phi);
0222
0223   double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
0224   double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
0225   double pz = pyjets.p[2][ip - 1];  // mom*cos(the) ;
0226   HepMC::FourVector ap(px, py, pz, ee);
0227   HepMC::GenParticle* APart = new HepMC::GenParticle(ap, particleID2, 1);
0228   APart->suggest_barcode(ip);
0229
0230   return APart;
0231 }