File indexing completed on 2024-04-06 12:13:54
0001
0002
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
0029
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
0043
0044
0045
0046
0047
0048 void pystrhad_();
0049 void pystfr_(int&);
0050
0051
0052 void pyglrhad_();
0053 void pyglfr_();
0054
0055
0056
0057 }
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
0081
0082
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)),
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
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
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
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
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
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
0167
0168
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
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
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
0211
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
0257
0258 if (lhe) {
0259 lheEvent()->fillEventInfo(event().get());
0260 lheEvent()->fillPdfInfo(&pdf);
0261 } else {
0262
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
0275
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
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319 event()->set_pdf_info(pdf);
0320
0321 HepMC::GenCrossSection xsec;
0322 double cs = pypars.pari[0];
0323 cs *= 1.0e9;
0324 double cserr = cs / sqrt(pypars.msti[4]);
0325 xsec.set_cross_section(cs, cserr);
0326 event()->set_cross_section(xsec);
0327
0328
0329
0330 if (lhe && std::abs(lheRunInfo()->getHEPRUP()->IDWTUP) == 4)
0331
0332 event()->weights().push_back(pyint1.vint[96] * 1.0e9);
0333 else
0334 event()->weights().push_back(pyint1.vint[96]);
0335
0336
0337
0338 event()->weights().push_back(1. / (pyint1.vint[98]));
0339
0340
0341
0342
0343 eventInfo() = std::make_unique<GenEventInfoProduct>(event().get());
0344
0345
0346
0347
0348 if (!lhe) {
0349 eventInfo()->setBinningValues(std::vector<double>(1, pypars.pari[16]));
0350 }
0351
0352
0353
0354 if (fImposeProperTime || pydat1.mstj[21] == 3 || pydat1.mstj[21] == 4)
0355 imposeProperTime();
0356
0357
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
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
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);
0398
0399 FortranCallback::getInstance()->resetIterationsPerEvent();
0400
0401
0402
0403
0404 if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
0405
0406 call_pygive("MSTJ(14)=-1");
0407 }
0408
0409 call_pyevnt();
0410
0411 if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
0412
0413 call_pygive("MSTJ(14)=1");
0414 int ierr = 0;
0415 if (fStopHadronsEnabled) {
0416 pystfr_(ierr);
0417 if (ierr != 0)
0418 {
0419 event().reset();
0420 return false;
0421 }
0422 }
0423 if (fGluinoHadronsEnabled)
0424 pyglfr_();
0425 }
0426
0427 if (pyint1.mint[50] != 0)
0428 {
0429 event().reset();
0430 return false;
0431 }
0432
0433 call_pyhepc(1);
0434 event().reset(pythia6_conv.read_next_event());
0435
0436
0437
0438 flushTmpStorage();
0439 fillTmpStorage();
0440
0441 return true;
0442 }
0443
0444 bool Pythia6Hadronizer::hadronize() {
0445 Pythia6Service::InstanceWrapper guard(fPy6Service);
0446
0447 FortranCallback::getInstance()->setLHEEvent(lheEvent());
0448 FortranCallback::getInstance()->resetIterationsPerEvent();
0449 if (fJetMatching) {
0450 fJetMatching->resetMatchingStatus();
0451 fJetMatching->beforeHadronisation(lheEvent());
0452 }
0453
0454
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
0465 lheEvent()->count(lhef::LHERunInfo::kSelected);
0466
0467 event().reset();
0468 return false;
0469 }
0470
0471
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)
0482 {
0483 event().reset();
0484 return false;
0485 }
0486 }
0487
0488 if (fGluinoHadronsEnabled)
0489 pyglfr_();
0490 }
0491
0492 if (pyint1.mint[50] != 0)
0493 {
0494 event().reset();
0495 return false;
0496 }
0497
0498 call_pyhepc(1);
0499 event().reset(pythia6_conv.read_next_event());
0500
0501
0502
0503 flushTmpStorage();
0504 fillTmpStorage();
0505
0506 return true;
0507 }
0508
0509 bool Pythia6Hadronizer::decay() { return true; }
0510
0511 bool Pythia6Hadronizer::residualDecay() {
0512
0513
0514 Pythia6Service::InstanceWrapper guard(fPy6Service);
0515
0516
0517
0518 int NPartsBeforeDecays = pyjets_local.n;
0519 int NPartsAfterDecays = event().get()->particles_size();
0520 int barcode = NPartsAfterDecays;
0521
0522
0523
0524
0525
0526
0527
0528
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;
0535
0536 int pdgid = part->pdg_id();
0537 int py6ptr = pycomp_(pdgid);
0538 if (pydat3.mdcy[0][py6ptr - 1] != 1)
0539 continue;
0540 int py6id = HepPID::translatePDTtoPythia(pdgid);
0541
0542
0543
0544
0545
0546 if (part->momentum().t() <= part->generated_mass()) {
0547 continue;
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;
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
0575
0576 int parent = 1;
0577 pydecy_(parent);
0578
0579
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);
0587
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++)
0609 {
0610 if (pyjets.k[2][iprt2] != parent) {
0611 parent = pyjets.k[2][iprt2];
0612 break;
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;
0635
0636
0637 event().get()->add_vertex(DecVtx);
0638 }
0639 }
0640
0641
0642
0643 if (pyjets_local.n != pyjets.n) {
0644
0645
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);
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);
0673
0674
0675
0676
0677
0678 fPy6Service->setPYUPDAParams(false);
0679
0680 FortranCallback::getInstance()->setLHERunInfo(lheRunInfo());
0681
0682 if (fStopHadronsEnabled) {
0683
0684 call_pygive("MSTP(111)=0");
0685 pystrhad_();
0686 call_pygive("MWID(302)=0");
0687 call_pygive("MDCY(302,1)=0");
0688
0689 }
0690
0691 if (fGluinoHadronsEnabled) {
0692
0693 call_pygive("MSTP(111)=0");
0694 pyglrhad_();
0695
0696
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
0714 call_pygive("MSTP(143)=1");
0715 }
0716
0717 return true;
0718 }
0719
0720 bool Pythia6Hadronizer::initializeForInternalPartons() {
0721 Pythia6Service::InstanceWrapper guard(fPy6Service);
0722
0723 fPy6Service->setPYUPDAParams(false);
0724
0725 if (fStopHadronsEnabled) {
0726
0727 call_pygive("MSTP(111)=0");
0728 pystrhad_();
0729 }
0730
0731 if (fGluinoHadronsEnabled) {
0732
0733 call_pygive("MSTP(111)=0");
0734 pyglrhad_();
0735 }
0736
0737 if (fInitialState == PP)
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
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
0753 call_pyinit("3mom", "e-", "p", 0.0);
0754 } else if (fInitialState == PositronProton) {
0755
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
0763 call_pyinit("3mom", "e+", "p", 0.0);
0764 } else {
0765
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
0782 int pyCode = pycomp_(PyID);
0783 if (pyCode > 0) {
0784 std::ostringstream pyCard;
0785 pyCard << "MDCY(" << pyCode << ",1)=0";
0786
0787
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
0814
0815
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
0834
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
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
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
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
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
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
0899 double cs = pypars.pari[0];
0900 cs *= 1.0e9;
0901 runInfo().setInternalXSec(cs);
0902
0903 }
0904
0905 call_pystat(1);
0906
0907 return;
0908 }
0909
0910 const char* Pythia6Hadronizer::classname() const { return "gen::Pythia6Hadronizer"; }
0911
0912 }