Pythia6Interface/plugins/Pythia6Hadronizer.cc

0001
0002 // -*- C++ -*-
0003
0004 #include "Pythia6Hadronizer.h"
0005
0006 #include "HepMC/GenEvent.h"
0007 #include "HepMC/PdfInfo.h"
0008 #include "HepMC/PythiaWrapper6_4.h"
0009 #include "HepMC/HEPEVT_Wrapper.h"
0010 #include "HepMC/IO_HEPEVT.h"
0011
0012 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "GeneratorInterface/Core/interface/FortranCallback.h"
0016 #include "GeneratorInterface/Core/interface/FortranInstance.h"
0017
0018 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0019
0020 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0021
0022 #include "GeneratorInterface/PartonShowerVeto/interface/JetMatching.h"
0023
0024 HepMC::IO_HEPEVT pythia6_conv;
0025
0026 #include "HepPID/ParticleIDTranslations.hh"
0027
0028 // NOTE: here a number of Pythia6 routines are declared,
0029 // plus some functionalities to pass around Pythia6 params
0030 //
0031 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
0032 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h"
0033
0034 #include <iomanip>
0035 #include <memory>
0036
0037 namespace gen {
0038
0039   extern "C" {
0040
0041   //
0042   // these two are NOT part of Pythi6 core code but are "custom" add-ons
0043   // we keep them interfaced here, rather than in GenExtensions, because
0044   // they tweak not at the ME level, but a step further, at the framgmentation
0045   //
0046   // stop-hadrons
0047   //
0048   void pystrhad_();    // init stop-hadrons (id's, names, charges...)
0049   void pystfr_(int&);  // tweaks fragmentation, fragments the string near to a stop,
0050                        // to form stop-hadron by producing a new q-qbar pair
0051   // gluino/r-hadrons
0052   void pyglrhad_();
0053   void pyglfr_();  // tweaks fragmentation, fragment the string near to a gluino,
0054                    // to form gluino-hadron, either by producing a new g-g pair,
0055                    // or two new q-qbar ones
0056
0057   }  // extern "C"
0058
0059   class Pythia6ServiceWithCallback : public Pythia6Service {
0060   public:
0061     Pythia6ServiceWithCallback(const edm::ParameterSet& ps) : Pythia6Service(ps) {}
0062
0063   private:
0064     void upInit() override { FortranCallback::getInstance()->fillHeader(); }
0065
0066     void upEvnt() override {
0067       FortranCallback::getInstance()->fillEvent();
0068       if (Pythia6Hadronizer::getJetMatching())
0069         Pythia6Hadronizer::getJetMatching()->beforeHadronisationExec();
0070     }
0071
0072     bool upVeto() override {
0073       if (!Pythia6Hadronizer::getJetMatching())
0074         return false;
0075
0076       if (!hepeup_.nup || Pythia6Hadronizer::getJetMatching()->isMatchingDone())
0077         return true;
0078
0079       bool retValue = Pythia6Hadronizer::getJetMatching()->match(nullptr, nullptr);
0080       // below is old code and a note of it
0081       // NOTE: I'm passing NULL pointers, instead of HepMC::GenEvent, etc.
0082       //retValur = Pythia6Hadronizer::getJetMatching()->match(0, 0, true);
0083       return retValue;
0084     }
0085   };
0086
0087   static struct {
0088     int n, npad, k[5][pyjets_maxn];
0089     double p[5][pyjets_maxn], v[5][pyjets_maxn];
0090   } pyjets_local;
0091
0092   JetMatching* Pythia6Hadronizer::fJetMatching = nullptr;
0093
0094   const std::vector<std::string> Pythia6Hadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6,
0095                                                                           gen::FortranInstance::kFortranInstance};
0096
0097   Pythia6Hadronizer::Pythia6Hadronizer(edm::ParameterSet const& ps)
0098       : BaseHadronizer(ps),
0099         fPy6Service(new Pythia6ServiceWithCallback(ps)),  // this will store py6 params for further settings
0100         fInitialState(PP),
0101         fCOMEnergy(ps.getParameter<double>("comEnergy")),
0102         fHepMCVerbosity(ps.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false)),
0103         fMaxEventsToPrint(ps.getUntrackedParameter<int>("maxEventsToPrint", 0)),
0104         fPythiaListVerbosity(ps.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
0105         fDisplayPythiaBanner(ps.getUntrackedParameter<bool>("displayPythiaBanner", false)),
0106         fDisplayPythiaCards(ps.getUntrackedParameter<bool>("displayPythiaCards", false)) {
0107     // J.Y.: the following 3 parameters are hacked "for a reason"
0108     //
0109     if (ps.exists("PPbarInitialState")) {
0110       if (fInitialState == PP) {
0111         fInitialState = PPbar;
0112         edm::LogInfo("GeneratorInterface|Pythia6Interface")
0113             << "Pythia6 will be initialized for PROTON-ANTIPROTON INITIAL STATE. "
0114             << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0115         std::cout << "Pythia6 will be initialized for PROTON-ANTIPROTON INITIAL STATE." << std::endl;
0116         std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0117       } else {
0118         // probably need to throw on attempt to override ?
0119       }
0120     } else if (ps.exists("ElectronPositronInitialState")) {
0121       if (fInitialState == PP) {
0122         fInitialState = ElectronPositron;
0123         edm::LogInfo("GeneratorInterface|Pythia6Interface")
0124             << "Pythia6 will be initialized for ELECTRON-POSITRON INITIAL STATE. "
0125             << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0126         std::cout << "Pythia6 will be initialized for ELECTRON-POSITRON INITIAL STATE." << std::endl;
0127         std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0128       } else {
0129         // probably need to throw on attempt to override ?
0130       }
0131     } else if (ps.exists("ElectronProtonInitialState")) {
0132       if (fInitialState == PP) {
0133         fInitialState = ElectronProton;
0134         fBeam1PZ =
0135             (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("electronMomentum");
0136         fBeam2PZ =
0137             (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("protonMomentum");
0138         edm::LogInfo("GeneratorInterface|Pythia6Interface")
0139             << "Pythia6 will be initialized for ELECTRON-PROTON INITIAL STATE. "
0140             << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0141         std::cout << "Pythia6 will be initialized for ELECTRON-PROTON INITIAL STATE." << std::endl;
0142         std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0143       } else {
0144         // probably need to throw on attempt to override ?
0145       }
0146     } else if (ps.exists("PositronProtonInitialState")) {
0147       if (fInitialState == PP) {
0148         fInitialState = PositronProton;
0149         fBeam1PZ =
0150             (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("positronMomentum");
0151         fBeam2PZ =
0152             (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("protonMomentum");
0153         edm::LogInfo("GeneratorInterface|Pythia6Interface")
0154             << "Pythia6 will be initialized for POSITRON-PROTON INITIAL STATE. "
0155             << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0156         std::cout << "Pythia6 will be initialized for POSITRON-PROTON INITIAL STATE." << std::endl;
0157         std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
0158       } else {
0159         // throw on unknown initial state !
0160         throw edm::Exception(edm::errors::Configuration, "Pythia6Interface")
0161             << " UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron, "
0162                "ElectronProton, and PositronProton \n";
0163       }
0164     }
0165
0166     // J.Y.: the following 4 params are "hacked", in the sense
0167     // that they're tracked but get in optionally;
0168     // this will be fixed once we update all applications
0169     //
0170
0171     fStopHadronsEnabled = false;
0172     if (ps.exists("stopHadrons"))
0173       fStopHadronsEnabled = ps.getParameter<bool>("stopHadrons");
0174
0175     fGluinoHadronsEnabled = false;
0176     if (ps.exists("gluinoHadrons"))
0177       fGluinoHadronsEnabled = ps.getParameter<bool>("gluinoHadrons");
0178
0179     fImposeProperTime = false;
0180     if (ps.exists("imposeProperTime")) {
0181       fImposeProperTime = ps.getParameter<bool>("imposeProperTime");
0182     }
0183
0184     fConvertToPDG = false;
0185     if (ps.exists("doPDGConvert"))
0186       fConvertToPDG = ps.getParameter<bool>("doPDGConvert");
0187
0188     if (ps.exists("jetMatching")) {
0189       edm::ParameterSet jmParams = ps.getUntrackedParameter<edm::ParameterSet>("jetMatching");
0190
0191       fJetMatching = JetMatching::create(jmParams).release();
0192     }
0193
0194     // first of all, silence Pythia6 banner printout, unless display requested
0195     //
0196     if (!fDisplayPythiaBanner) {
0197       if (!call_pygive("MSTU(12)=12345")) {
0198         throw edm::Exception(edm::errors::Configuration, "PythiaError") << " Pythia did not accept MSTU(12)=12345";
0199       }
0200     }
0201
0202     // silence printouts from PYGIVE, unless display requested
0203     //
0204     if (!fDisplayPythiaCards) {
0205       if (!call_pygive("MSTU(13)=0")) {
0206         throw edm::Exception(edm::errors::Configuration, "PythiaError") << " Pythia did not accept MSTU(13)=0";
0207       }
0208     }
0209
0210     // tmp stuff to deal with EvtGen corrupting pyjets
0211     // NPartsBeforeDecays = 0;
0212     flushTmpStorage();
0213   }
0214
0215   Pythia6Hadronizer::~Pythia6Hadronizer() {
0216     if (fPy6Service != nullptr)
0217       delete fPy6Service;
0218     if (fJetMatching != nullptr)
0219       delete fJetMatching;
0220   }
0221
0222   void Pythia6Hadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { fPy6Service->setRandomEngine(v); }
0223
0224   void Pythia6Hadronizer::flushTmpStorage() {
0225     pyjets_local.n = 0;
0226     pyjets_local.npad = 0;
0227     for (int ip = 0; ip < pyjets_maxn; ip++) {
0228       for (int i = 0; i < 5; i++) {
0229         pyjets_local.k[i][ip] = 0;
0230         pyjets_local.p[i][ip] = 0.;
0231         pyjets_local.v[i][ip] = 0.;
0232       }
0233     }
0234     return;
0235   }
0236
0237   void Pythia6Hadronizer::fillTmpStorage() {
0238     pyjets_local.n = pyjets.n;
0239     pyjets_local.npad = pyjets.npad;
0240     for (int ip = 0; ip < pyjets_maxn; ip++) {
0241       for (int i = 0; i < 5; i++) {
0242         pyjets_local.k[i][ip] = pyjets.k[i][ip];
0243         pyjets_local.p[i][ip] = pyjets.p[i][ip];
0244         pyjets_local.v[i][ip] = pyjets.v[i][ip];
0245       }
0246     }
0247
0248     return;
0249   }
0250
0251   void Pythia6Hadronizer::finalizeEvent() {
0252     bool lhe = lheEvent() != nullptr;
0253
0254     HepMC::PdfInfo pdf;
0255
0256     // if we are in hadronizer mode, we can pass on information from
0257     // the LHE input
0258     if (lhe) {
0259       lheEvent()->fillEventInfo(event().get());
0260       lheEvent()->fillPdfInfo(&pdf);
0261     } else {
0262       // filling in factorization "Q scale" now! pthat moved to binningValues()
0263       //
0264
0265       if (event()->signal_process_id() <= 0)
0266         event()->set_signal_process_id(pypars.msti[0]);
0267       if (event()->event_scale() <= 0)
0268         event()->set_event_scale(pypars.pari[22]);
0269       if (event()->alphaQED() <= 0)
0270         event()->set_alphaQED(pyint1.vint[56]);
0271       if (event()->alphaQCD() <= 0)
0272         event()->set_alphaQCD(pyint1.vint[57]);
0273
0274       // get pdf info directly from Pythia6 and set it up into HepMC::GenEvent
0275       // S. Mrenna: Prefer vint block
0276       //
0277       if (pdf.id1() <= 0)
0278         pdf.set_id1(pyint1.mint[14] == 21 ? 0 : pyint1.mint[14]);
0279       if (pdf.id2() <= 0)
0280         pdf.set_id2(pyint1.mint[15] == 21 ? 0 : pyint1.mint[15]);
0281       if (pdf.x1() <= 0)
0282         pdf.set_x1(pyint1.vint[40]);
0283       if (pdf.x2() <= 0)
0284         pdf.set_x2(pyint1.vint[41]);
0285       if (pdf.pdf1() <= 0)
0286         pdf.set_pdf1(pyint1.vint[38] / pyint1.vint[40]);
0287       if (pdf.pdf2() <= 0)
0288         pdf.set_pdf2(pyint1.vint[39] / pyint1.vint[41]);
0289       if (pdf.scalePDF() <= 0)
0290         pdf.set_scalePDF(pyint1.vint[50]);
0291     }
0292
0293     /* 9/9/2010 - JVY: This is the old piece of code - I can't remember why we implemented it this way.
0294                    However, it's causing problems with pdf1 & pdf2 when processing LHE samples,
0295            specifically, because both are set to -1, it tries to fill with Py6 numbers that
0296            are NOT valid/right at this point !
0297            In general, for LHE/ME event processing we should implement the correct calculation
0298            of the pdf's, rather than using py6 ones.
0299
0300    // filling in factorization "Q scale" now! pthat moved to binningValues()
0301
0302    if (!lhe || event()->signal_process_id() < 0) event()->set_signal_process_id( pypars.msti[0] );
0303    if (!lhe || event()->event_scale() < 0)       event()->set_event_scale( pypars.pari[22] );
0304    if (!lhe || event()->alphaQED() < 0)          event()->set_alphaQED( pyint1.vint[56] );
0305    if (!lhe || event()->alphaQCD() < 0)          event()->set_alphaQCD( pyint1.vint[57] );
0306
0307    // get pdf info directly from Pythia6 and set it up into HepMC::GenEvent
0308    // S. Mrenna: Prefer vint block
0309    //
0310    if (!lhe || pdf.id1() < 0)      pdf.set_id1( pyint1.mint[14] == 21 ? 0 : pyint1.mint[14] );
0311    if (!lhe || pdf.id2() < 0)      pdf.set_id2( pyint1.mint[15] == 21 ? 0 : pyint1.mint[15] );
0312    if (!lhe || pdf.x1() < 0)       pdf.set_x1( pyint1.vint[40] );
0313    if (!lhe || pdf.x2() < 0)       pdf.set_x2( pyint1.vint[41] );
0314    if (!lhe || pdf.pdf1() < 0)     pdf.set_pdf1( pyint1.vint[38] / pyint1.vint[40] );
0315    if (!lhe || pdf.pdf2() < 0)     pdf.set_pdf2( pyint1.vint[39] / pyint1.vint[41] );
0316    if (!lhe || pdf.scalePDF() < 0) pdf.set_scalePDF( pyint1.vint[50] );
0317 */
0318
0319     event()->set_pdf_info(pdf);
0320
0321     HepMC::GenCrossSection xsec;
0322     double cs = pypars.pari[0];  // cross section in mb
0323     cs *= 1.0e9;                 // translate to pb
0324     double cserr = cs / sqrt(pypars.msti[4]);
0325     xsec.set_cross_section(cs, cserr);
0326     event()->set_cross_section(xsec);
0327
0328     // this is "standard" Py6 event weight (corresponds to PYINT1/VINT(97)
0329     //
0330     if (lhe && std::abs(lheRunInfo()->getHEPRUP()->IDWTUP) == 4)
0331       // translate mb to pb (CMS/Gen "convention" as of May 2009)
0332       event()->weights().push_back(pyint1.vint[96] * 1.0e9);
0333     else
0334       event()->weights().push_back(pyint1.vint[96]);
0335     //
0336     // this is event weight as 1./VINT(99) (PYINT1/VINT(99) is returned by the PYEVWT)
0337     //
0338     event()->weights().push_back(1. / (pyint1.vint[98]));
0339
0340     // now create the GenEventInfo product from the GenEvent and fill
0341     // the missing pieces
0342
0343     eventInfo() = std::make_unique<GenEventInfoProduct>(event().get());
0344
0345     // in Pythia6 pthat is used to subdivide samples into different bins
0346     // in LHE mode the binning is done by the external ME generator
0347     // which is likely not pthat, so only filling it for Py6 internal mode
0348     if (!lhe) {
0349       eventInfo()->setBinningValues(std::vector<double>(1, pypars.pari[16]));
0350     }
0351
0352     // here we treat long-lived particles
0353     //
0354     if (fImposeProperTime || pydat1.mstj[21] == 3 || pydat1.mstj[21] == 4)
0355       imposeProperTime();
0356
0357     // convert particle IDs Py6->PDG, if requested
0358     if (fConvertToPDG) {
0359       for (HepMC::GenEvent::particle_iterator part = event()->particles_begin(); part != event()->particles_end();
0360            ++part) {
0361         (*part)->set_pdg_id(HepPID::translatePythiatoPDT((*part)->pdg_id()));
0362       }
0363     }
0364
0365     // service printouts, if requested
0366     //
0367     if (fMaxEventsToPrint > 0) {
0368       fMaxEventsToPrint--;
0369       if (fPythiaListVerbosity)
0370         call_pylist(fPythiaListVerbosity);
0371       if (fHepMCVerbosity) {
0372         std::cout << "Event process = " << pypars.msti[0] << std::endl << "----------------------" << std::endl;
0373         event()->print();
0374       }
0375     }
0376
0377     // dump of all settings after all initializations
0378     if (fDisplayPythiaCards) {
0379       fDisplayPythiaCards = false;
0380       call_pylist(12);
0381       call_pylist(13);
0382       std::cout << "\n PYPARS \n" << std::endl;
0383       std::cout << std::setw(5) << std::fixed << "I" << std::setw(10) << std::fixed << "MSTP(I)" << std::setw(16)
0384                 << std::fixed << "PARP(I)" << std::setw(10) << std::fixed << "MSTI(I)" << std::setw(16) << std::fixed
0385                 << "PARI(I)" << std::endl;
0386       for (unsigned int ind = 0; ind < 200; ind++) {
0387         std::cout << std::setw(5) << std::fixed << ind + 1 << std::setw(10) << std::fixed << pypars.mstp[ind]
0388                   << std::setw(16) << std::fixed << pypars.parp[ind] << std::setw(10) << std::fixed << pypars.msti[ind]
0389                   << std::setw(16) << std::fixed << pypars.pari[ind] << std::endl;
0390       }
0391     }
0392
0393     return;
0394   }
0395
0396   bool Pythia6Hadronizer::generatePartonsAndHadronize() {
0397     Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0398
0399     FortranCallback::getInstance()->resetIterationsPerEvent();
0400
0401     // generate event with Pythia6
0402     //
0403
0404     if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
0405       // call_pygive("MSTJ(1)=-1");
0406       call_pygive("MSTJ(14)=-1");
0407     }
0408
0409     call_pyevnt();
0410
0411     if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
0412       // call_pygive("MSTJ(1)=1");
0413       call_pygive("MSTJ(14)=1");
0414       int ierr = 0;
0415       if (fStopHadronsEnabled) {
0416         pystfr_(ierr);
0417         if (ierr != 0)  // skip failed events
0418         {
0419           event().reset();
0420           return false;
0421         }
0422       }
0423       if (fGluinoHadronsEnabled)
0424         pyglfr_();
0425     }
0426
0427     if (pyint1.mint[50] != 0)  // skip event if py6 considers it bad
0428     {
0429       event().reset();
0430       return false;
0431     }
0432
0433     call_pyhepc(1);
0434     event().reset(pythia6_conv.read_next_event());
0435
0436     // this is to deal with post-gen tools & residualDecay() that may reuse PYJETS
0437     //
0438     flushTmpStorage();
0439     fillTmpStorage();
0440
0441     return true;
0442   }
0443
0444   bool Pythia6Hadronizer::hadronize() {
0445     Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0446
0447     FortranCallback::getInstance()->setLHEEvent(lheEvent());
0448     FortranCallback::getInstance()->resetIterationsPerEvent();
0449     if (fJetMatching) {
0450       fJetMatching->resetMatchingStatus();
0451       fJetMatching->beforeHadronisation(lheEvent());
0452     }
0453
0454     // generate event with Pythia6
0455     //
0456     if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
0457       call_pygive("MSTJ(1)=-1");
0458       call_pygive("MSTJ(14)=-1");
0459     }
0460
0461     call_pyevnt();
0462
0463     if (FortranCallback::getInstance()->getIterationsPerEvent() > 1 || hepeup_.nup <= 0 || pypars.msti[0] == 1) {
0464       // update LHE matching statistics
0465       lheEvent()->count(lhef::LHERunInfo::kSelected);
0466
0467       event().reset();
0468       return false;
0469     }
0470
0471     // update LHE matching statistics
0472     //
0473     lheEvent()->count(lhef::LHERunInfo::kAccepted);
0474
0475     if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
0476       call_pygive("MSTJ(1)=1");
0477       call_pygive("MSTJ(14)=1");
0478       int ierr = 0;
0479       if (fStopHadronsEnabled) {
0480         pystfr_(ierr);
0481         if (ierr != 0)  // skip failed events
0482         {
0483           event().reset();
0484           return false;
0485         }
0486       }
0487
0488       if (fGluinoHadronsEnabled)
0489         pyglfr_();
0490     }
0491
0492     if (pyint1.mint[50] != 0)  // skip event if py6 considers it bad
0493     {
0494       event().reset();
0495       return false;
0496     }
0497
0498     call_pyhepc(1);
0499     event().reset(pythia6_conv.read_next_event());
0500
0501     // this is to deal with post-gen tools and/or residualDecay(), that may reuse PYJETS
0502     //
0503     flushTmpStorage();
0504     fillTmpStorage();
0505
0506     return true;
0507   }
0508
0509   bool Pythia6Hadronizer::decay() { return true; }
0510
0511   bool Pythia6Hadronizer::residualDecay() {
0512     // event().get()->print();
0513
0514     Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0515
0516     // int nDocLines = pypars.msti[3];
0517
0518     int NPartsBeforeDecays = pyjets_local.n;
0519     int NPartsAfterDecays = event().get()->particles_size();
0520     int barcode = NPartsAfterDecays;
0521
0522     // JVY: well, in principle, it's not a 100% fair to go up to NPartsBeforeDecays,
0523     // because Photos will attach gamma's to existing vertexes, i.e. in the middle
0524     // of the event rather than at the end; but this will only shift poiters down,
0525     // so we'll be going again over a few "original" particle...
0526     // in the alternative, we may go all the way up to the beginning of the event
0527     // and re-check if anything remains to decay, that's fine even if it'll take
0528     // some extra CPU...
0529
0530     for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) {
0531       HepMC::GenParticle* part = event().get()->barcode_to_particle(ipart);
0532       int status = part->status();
0533       if (status != 1)
0534         continue;  // check only "stable" particles,
0535                    // as some undecayed may still be marked as such
0536       int pdgid = part->pdg_id();
0537       int py6ptr = pycomp_(pdgid);
0538       if (pydat3.mdcy[0][py6ptr - 1] != 1)
0539         continue;  // particle is not expected to decay
0540       int py6id = HepPID::translatePDTtoPythia(pdgid);
0541       //
0542       // first, will need to zero out, then fill up PYJETS
0543       // I better do it directly (by hands) rather than via py1ent
0544       // - to avoid (additional) mass smearing
0545       //
0546       if (part->momentum().t() <= part->generated_mass()) {
0547         continue;  // e==m --> 0-momentum, nothing to decay...
0548       } else {
0549         pyjets.n = 0;
0550         for (int i = 0; i < 5; i++) {
0551           pyjets.k[i][0] = 0;
0552           pyjets.p[i][0] = 0.;
0553           pyjets.v[i][0] = 0.;
0554         }
0555         pyjets.k[0][0] = 1;
0556         pyjets.k[1][0] = py6id;
0557         pyjets.p[4][0] = part->generated_mass();
0558         pyjets.p[3][0] = part->momentum().t();
0559         pyjets.p[0][0] = part->momentum().x();
0560         pyjets.p[1][0] = part->momentum().y();
0561         pyjets.p[2][0] = part->momentum().z();
0562         HepMC::GenVertex* prod_vtx = part->production_vertex();
0563         if (!prod_vtx)
0564           continue;  // in principle, should never happen but...
0565         pyjets.v[0][0] = prod_vtx->position().x();
0566         pyjets.v[1][0] = prod_vtx->position().y();
0567         pyjets.v[2][0] = prod_vtx->position().z();
0568         pyjets.v[3][0] = prod_vtx->position().t();
0569         pyjets.v[4][0] = 0.;
0570         pyjets.n = 1;
0571         pyjets.npad = pyjets_local.npad;
0572       }
0573
0574       // now call Py6 decay routine
0575       //
0576       int parent = 1;  // since we always pass to Py6 a single particle
0577       pydecy_(parent);
0578
0579       // now attach decay products to mother
0580       //
0581       for (int iprt1 = 1; iprt1 < pyjets.n; iprt1++) {
0582         part->set_status(2);
0583
0584         HepMC::GenVertex* DecVtx = new HepMC::GenVertex(
0585             HepMC::FourVector(pyjets.v[0][iprt1], pyjets.v[1][iprt1], pyjets.v[2][iprt1], pyjets.v[3][iprt1]));
0586         DecVtx->add_particle_in(part);  // this will cleanup end_vertex if exists, replace with the new one
0587                                         // I presume (vtx) barcode will be given automatically
0588
0589         HepMC::FourVector pmom(pyjets.p[0][iprt1], pyjets.p[1][iprt1], pyjets.p[2][iprt1], pyjets.p[3][iprt1]);
0590
0591         int dstatus = 0;
0592         if (pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10) {
0593           dstatus = 1;
0594         } else if (pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20) {
0595           dstatus = 2;
0596         } else if (pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30) {
0597           dstatus = 3;
0598         } else if (pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100) {
0599           dstatus = pyjets.k[0][iprt1];
0600         }
0601         HepMC::GenParticle* daughter =
0602             new HepMC::GenParticle(pmom, HepPID::translatePythiatoPDT(pyjets.k[1][iprt1]), dstatus);
0603         barcode++;
0604         daughter->suggest_barcode(barcode);
0605         DecVtx->add_particle_out(daughter);
0606
0607         int iprt2;
0608         for (iprt2 = iprt1 + 1; iprt2 < pyjets.n; iprt2++)  // the pointer is shifted by -1, c++ style
0609         {
0610           if (pyjets.k[2][iprt2] != parent) {
0611             parent = pyjets.k[2][iprt2];
0612             break;  // another parent particle; reset & break the loop
0613           }
0614
0615           HepMC::FourVector pmomN(pyjets.p[0][iprt2], pyjets.p[1][iprt2], pyjets.p[2][iprt2], pyjets.p[3][iprt2]);
0616
0617           dstatus = 0;
0618           if (pyjets.k[0][iprt2] >= 1 && pyjets.k[0][iprt2] <= 10) {
0619             dstatus = 1;
0620           } else if (pyjets.k[0][iprt2] >= 11 && pyjets.k[0][iprt2] <= 20) {
0621             dstatus = 2;
0622           } else if (pyjets.k[0][iprt2] >= 21 && pyjets.k[0][iprt2] <= 30) {
0623             dstatus = 3;
0624           } else if (pyjets.k[0][iprt2] >= 31 && pyjets.k[0][iprt2] <= 100) {
0625             dstatus = pyjets.k[0][iprt2];
0626           }
0627           HepMC::GenParticle* daughterN =
0628               new HepMC::GenParticle(pmomN, HepPID::translatePythiatoPDT(pyjets.k[1][iprt2]), dstatus);
0629           barcode++;
0630           daughterN->suggest_barcode(barcode);
0631           DecVtx->add_particle_out(daughterN);
0632         }
0633
0634         iprt1 = iprt2 - 1;  // reset counter such that it doesn't go over the same child more than once
0635                             // don't forget to offset back into c++ counting, as it's already +1 forward
0636
0637         event().get()->add_vertex(DecVtx);
0638       }
0639     }
0640
0641     // now restore the very original Py6 event record
0642     //
0643     if (pyjets_local.n != pyjets.n) {
0644       // restore pyjets to its state as it was before external decays -
0645       // might have been jammed by action above or by py1ent calls in EvtGen
0646       pyjets.n = pyjets_local.n;
0647       pyjets.npad = pyjets_local.npad;
0648       for (int ip = 0; ip < pyjets_local.n; ip++) {
0649         for (int i = 0; i < 5; i++) {
0650           pyjets.k[i][ip] = pyjets_local.k[i][ip];
0651           pyjets.p[i][ip] = pyjets_local.p[i][ip];
0652           pyjets.v[i][ip] = pyjets_local.v[i][ip];
0653         }
0654       }
0655     }
0656
0657     return true;
0658   }
0659
0660   bool Pythia6Hadronizer::readSettings(int key) {
0661     Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0662
0663     fPy6Service->setGeneralParams();
0664     if (key == 0)
0665       fPy6Service->setCSAParams();
0666     fPy6Service->setSLHAParams();
0667
0668     return true;
0669   }
0670
0671   bool Pythia6Hadronizer::initializeForExternalPartons() {
0672     Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0673
0674     // note: CSA mode is NOT supposed to work with external partons !!!
0675
0676     // fPy6Service->setGeneralParams();
0677
0678     fPy6Service->setPYUPDAParams(false);
0679
0680     FortranCallback::getInstance()->setLHERunInfo(lheRunInfo());
0681
0682     if (fStopHadronsEnabled) {
0683       // overwrite mstp(111), no matter what
0684       call_pygive("MSTP(111)=0");
0685       pystrhad_();
0686       call_pygive("MWID(302)=0");    // I don't know if this is specific to processing ME/LHE only,
0687       call_pygive("MDCY(302,1)=0");  // or this should also be the case for full event...
0688                                      // anyway, this comes from experience of processing MG events
0689     }
0690
0691     if (fGluinoHadronsEnabled) {
0692       // overwrite mstp(111), no matter what
0693       call_pygive("MSTP(111)=0");
0694       pyglrhad_();
0695       //call_pygive("MWID(309)=0");
0696       //call_pygive("MDCY(309,1)=0");
0697     }
0698
0699     call_pyinit("USER", "", "", 0.0);
0700
0701     fPy6Service->setPYUPDAParams(true);
0702
0703     std::vector<std::string> slha = lheRunInfo()->findHeader("slha");
0704     if (!slha.empty()) {
0705       edm::LogInfo("Generator|LHEInterface") << "Pythia6 hadronisation found an SLHA header, "
0706                                              << "will be passed on to Pythia." << std::endl;
0707       fPy6Service->setSLHAFromHeader(slha);
0708       fPy6Service->closeSLHA();
0709     }
0710
0711     if (fJetMatching) {
0712       fJetMatching->init(lheRunInfo());
0713       // FIXME: the jet matching routine might not be interested in PS callback
0714       call_pygive("MSTP(143)=1");
0715     }
0716
0717     return true;
0718   }
0719
0720   bool Pythia6Hadronizer::initializeForInternalPartons() {
0721     Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0722
0723     fPy6Service->setPYUPDAParams(false);
0724
0725     if (fStopHadronsEnabled) {
0726       // overwrite mstp(111), no matter what
0727       call_pygive("MSTP(111)=0");
0728       pystrhad_();
0729     }
0730
0731     if (fGluinoHadronsEnabled) {
0732       // overwrite mstp(111), no matter what
0733       call_pygive("MSTP(111)=0");
0734       pyglrhad_();
0735     }
0736
0737     if (fInitialState == PP)  // default
0738     {
0739       call_pyinit("CMS", "p", "p", fCOMEnergy);
0740     } else if (fInitialState == PPbar) {
0741       call_pyinit("CMS", "p", "pbar", fCOMEnergy);
0742     } else if (fInitialState == ElectronPositron) {
0743       call_pyinit("CMS", "e+", "e-", fCOMEnergy);
0744     } else if (fInitialState == ElectronProton) {
0745       // set p(1,i) & p(2,i) for the beams in pyjets !
0746       pyjets.p[0][0] = 0.;
0747       pyjets.p[1][0] = 0.;
0748       pyjets.p[2][0] = fBeam1PZ;
0749       pyjets.p[0][1] = 0.;
0750       pyjets.p[1][1] = 0.;
0751       pyjets.p[2][1] = fBeam2PZ;
0752       // call "3mon" frame & 0.0 win
0753       call_pyinit("3mom", "e-", "p", 0.0);
0754     } else if (fInitialState == PositronProton) {
0755       // set p(1,i) & p(2,i) for the beams in pyjets !
0756       pyjets.p[0][0] = 0.;
0757       pyjets.p[1][0] = 0.;
0758       pyjets.p[2][0] = fBeam1PZ;
0759       pyjets.p[0][1] = 0.;
0760       pyjets.p[1][1] = 0.;
0761       pyjets.p[2][1] = fBeam2PZ;
0762       // call "3mon" frame & 0.0 win
0763       call_pyinit("3mom", "e+", "p", 0.0);
0764     } else {
0765       // throw on unknown initial state !
0766       throw edm::Exception(edm::errors::Configuration, "Pythia6Interface")
0767           << " UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron, ElectronProton, "
0768              "and PositronProton \n";
0769     }
0770
0771     fPy6Service->setPYUPDAParams(true);
0772
0773     fPy6Service->closeSLHA();
0774
0775     return true;
0776   }
0777
0778   bool Pythia6Hadronizer::declareStableParticles(const std::vector<int>& pdg) {
0779     for (unsigned int i = 0; i < pdg.size(); i++) {
0780       int PyID = HepPID::translatePDTtoPythia(pdg[i]);
0781       // int PyID = pdg[i];
0782       int pyCode = pycomp_(PyID);
0783       if (pyCode > 0) {
0784         std::ostringstream pyCard;
0785         pyCard << "MDCY(" << pyCode << ",1)=0";
0786         /* this is a test printout...
0787          std::cout << "pdg= " << pdg[i] << " " << pyCard.str() << std::endl;
0788 */
0789         call_pygive(pyCard.str());
0790       }
0791     }
0792
0793     return true;
0794   }
0795
0796   bool Pythia6Hadronizer::declareSpecialSettings(const std::vector<std::string>& settings) {
0797     for (unsigned int iss = 0; iss < settings.size(); iss++) {
0798       if (settings[iss].find("QED-brem-off") == std::string::npos)
0799         continue;
0800       size_t fnd1 = settings[iss].find(':');
0801       if (fnd1 == std::string::npos)
0802         continue;
0803
0804       std::string value = settings[iss].substr(fnd1 + 1);
0805
0806       if (value == "all") {
0807         call_pygive("MSTJ(41)=3");
0808       } else {
0809         int number = atoi(value.c_str());
0810         int PyID = HepPID::translatePDTtoPythia(number);
0811         int pyCode = pycomp_(PyID);
0812         if (pyCode > 0) {
0813           // first of all, check if mstj(39) is 0 or if we're trying to override user's setting
0814           // if so, throw an exception and stop, because otherwise the user will get behaviour
0815           // that's different from what she/he expects !
0816           if (pydat1_.mstj[38] > 0 && pydat1_.mstj[38] != pyCode) {
0817             throw edm::Exception(edm::errors::Configuration, "Pythia6Interface")
0818                 << " Fatal conflict: \n mandatory internal directive to set MSTJ(39)=" << pyCode
0819                 << " overrides user setting MSTJ(39)=" << pydat1_.mstj[38]
0820                 << " - user will not get expected behaviour \n";
0821           }
0822           std::ostringstream pyCard;
0823           pyCard << "MSTJ(39)=" << pyCode;
0824           call_pygive(pyCard.str());
0825         }
0826       }
0827     }
0828
0829     return true;
0830   }
0831
0832   void Pythia6Hadronizer::imposeProperTime() {
0833     // this is practically a copy/paste of the original code by J.Alcaraz,
0834     // taken directly from PythiaSource
0835
0836     int dumm = 0;
0837     HepMC::GenEvent::vertex_const_iterator vbegin = event()->vertices_begin();
0838     HepMC::GenEvent::vertex_const_iterator vend = event()->vertices_end();
0839     HepMC::GenEvent::vertex_const_iterator vitr = vbegin;
0840     for (; vitr != vend; ++vitr) {
0841       HepMC::GenVertex::particle_iterator pbegin = (*vitr)->particles_begin(HepMC::children);
0842       HepMC::GenVertex::particle_iterator pend = (*vitr)->particles_end(HepMC::children);
0843       HepMC::GenVertex::particle_iterator pitr = pbegin;
0844       for (; pitr != pend; ++pitr) {
0845         if ((*pitr)->end_vertex())
0846           continue;
0847         if ((*pitr)->status() != 1)
0848           continue;
0849
0850         int pdgcode = abs((*pitr)->pdg_id());
0851         // Do nothing if the particle is not expected to decay
0852         if (pydat3.mdcy[0][pycomp_(pdgcode) - 1] != 1)
0853           continue;
0854
0855         double ctau = pydat2.pmas[3][pycomp_(pdgcode) - 1];
0856         HepMC::FourVector mom = (*pitr)->momentum();
0857         HepMC::FourVector vin = (*vitr)->position();
0858         double x = 0.;
0859         double y = 0.;
0860         double z = 0.;
0861         double t = 0.;
0862         bool decayInRange = false;
0863         while (!decayInRange) {
0864           double unif_rand = fPy6Service->call(pyr_, &dumm);
0865           // Value of 0 is excluded, so following line is OK
0866           double proper_length = -ctau * log(unif_rand);
0867           double factor = proper_length / mom.m();
0868           x = vin.x() + factor * mom.px();
0869           y = vin.y() + factor * mom.py();
0870           z = vin.z() + factor * mom.pz();
0871           t = vin.t() + factor * mom.e();
0872           // Decay must be happen outside a cylindrical region
0873           if (pydat1.mstj[21] == 4) {
0874             if (std::sqrt(x * x + y * y) > pydat1.parj[72] || fabs(z) > pydat1.parj[73])
0875               decayInRange = true;
0876             // Decay must be happen outside a given sphere
0877           } else if (pydat1.mstj[21] == 3) {
0878             if (std::sqrt(x * x + y * y + z * z) > pydat1.parj[71])
0879               decayInRange = true;
0880           }
0881           // Decay is always OK otherwise
0882           else {
0883             decayInRange = true;
0884           }
0885         }
0886
0887         HepMC::GenVertex* vdec = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t));
0888         event()->add_vertex(vdec);
0889         vdec->add_particle_in((*pitr));
0890       }
0891     }
0892
0893     return;
0894   }
0895
0896   void Pythia6Hadronizer::statistics() {
0897     if (!runInfo().internalXSec()) {
0898       // set xsec if not already done (e.g. from LHE cross section collector)
0899       double cs = pypars.pari[0];  // cross section in mb
0900       cs *= 1.0e9;                 // translate to pb (CMS/Gen "convention" as of May 2009)
0901       runInfo().setInternalXSec(cs);
0902       // FIXME: can we get the xsec statistical error somewhere?
0903     }
0904
0905     call_pystat(1);
0906
0907     return;
0908   }
0909
0910   const char* Pythia6Hadronizer::classname() const { return "gen::Pythia6Hadronizer"; }
0911
0912 }  // namespace gen