File indexing completed on 2024-11-28 03:54:33
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include "GeneratorInterface/EvtGenInterface/plugins/EvtGenUserModels/EvtModelUserReg.h"
0013 #include "GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h"
0014
0015 #include "GeneratorInterface/EvtGenInterface/interface/EvtGenFactory.h"
0016 #include "GeneratorInterface/EvtGenInterface/interface/myEvtRandomEngine.h"
0017
0018 #include "FWCore/PluginManager/interface/PluginManager.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/PluginManager/interface/ModuleDef.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024 #include "FWCore/Utilities/interface/EDMException.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/Run.h"
0027 #include "FWCore/ParameterSet/interface/FileInPath.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029
0030 #include "CLHEP/Random/RandFlat.h"
0031 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0032
0033 #include "EvtGen/EvtGen.hh"
0034 #include "EvtGenExternal/EvtExternalGenList.hh"
0035 #include "EvtGenBase/EvtAbsRadCorr.hh"
0036 #include "EvtGenBase/EvtDecayBase.hh"
0037 #include "EvtGenBase/EvtId.hh"
0038 #include "EvtGenBase/EvtDecayTable.hh"
0039 #include "EvtGenBase/EvtParticle.hh"
0040 #include "EvtGenBase/EvtParticleFactory.hh"
0041 #include "EvtGenBase/EvtHepMCEvent.hh"
0042
0043 #include "HepPID/ParticleIDTranslations.hh"
0044
0045 #include "TString.h"
0046 #include <string>
0047 #include <cstdlib>
0048 #include <cstdio>
0049
0050 using namespace gen;
0051 using namespace edm;
0052
0053 CLHEP::HepRandomEngine* EvtGenInterface::fRandomEngine;
0054
0055 EvtGenInterface::EvtGenInterface(const ParameterSet& pset) {
0056 fPSet = new ParameterSet(pset);
0057 the_engine = new myEvtRandomEngine(nullptr);
0058 }
0059
0060 void EvtGenInterface::SetDefault_m_PDGs() {
0061
0062
0063
0064
0065
0066 m_PDGs.push_back(300553);
0067 m_PDGs.push_back(511);
0068 m_PDGs.push_back(521);
0069 m_PDGs.push_back(523);
0070 m_PDGs.push_back(513);
0071 m_PDGs.push_back(533);
0072 m_PDGs.push_back(531);
0073
0074 m_PDGs.push_back(15);
0075
0076 m_PDGs.push_back(413);
0077 m_PDGs.push_back(423);
0078 m_PDGs.push_back(433);
0079 m_PDGs.push_back(411);
0080 m_PDGs.push_back(421);
0081 m_PDGs.push_back(431);
0082 m_PDGs.push_back(10411);
0083 m_PDGs.push_back(10421);
0084 m_PDGs.push_back(10413);
0085 m_PDGs.push_back(10423);
0086 m_PDGs.push_back(20413);
0087 m_PDGs.push_back(20423);
0088
0089 m_PDGs.push_back(415);
0090 m_PDGs.push_back(425);
0091 m_PDGs.push_back(10431);
0092 m_PDGs.push_back(20433);
0093 m_PDGs.push_back(10433);
0094 m_PDGs.push_back(435);
0095
0096 m_PDGs.push_back(310);
0097 m_PDGs.push_back(311);
0098 m_PDGs.push_back(313);
0099 m_PDGs.push_back(323);
0100 m_PDGs.push_back(10321);
0101 m_PDGs.push_back(10311);
0102 m_PDGs.push_back(10313);
0103 m_PDGs.push_back(10323);
0104 m_PDGs.push_back(20323);
0105 m_PDGs.push_back(20313);
0106 m_PDGs.push_back(325);
0107 m_PDGs.push_back(315);
0108
0109 m_PDGs.push_back(100313);
0110 m_PDGs.push_back(100323);
0111 m_PDGs.push_back(30313);
0112 m_PDGs.push_back(30323);
0113 m_PDGs.push_back(30343);
0114 m_PDGs.push_back(30353);
0115 m_PDGs.push_back(30363);
0116
0117 m_PDGs.push_back(111);
0118 m_PDGs.push_back(221);
0119 m_PDGs.push_back(113);
0120 m_PDGs.push_back(213);
0121 m_PDGs.push_back(223);
0122 m_PDGs.push_back(331);
0123 m_PDGs.push_back(333);
0124 m_PDGs.push_back(20213);
0125 m_PDGs.push_back(20113);
0126 m_PDGs.push_back(215);
0127 m_PDGs.push_back(115);
0128 m_PDGs.push_back(10213);
0129 m_PDGs.push_back(10113);
0130 m_PDGs.push_back(9000111);
0131 m_PDGs.push_back(9000211);
0132 m_PDGs.push_back(9010221);
0133 m_PDGs.push_back(10221);
0134 m_PDGs.push_back(20223);
0135 m_PDGs.push_back(20333);
0136 m_PDGs.push_back(225);
0137 m_PDGs.push_back(9020221);
0138 m_PDGs.push_back(335);
0139 m_PDGs.push_back(10223);
0140 m_PDGs.push_back(10333);
0141 m_PDGs.push_back(100213);
0142 m_PDGs.push_back(100113);
0143
0144 m_PDGs.push_back(441);
0145 m_PDGs.push_back(100441);
0146 m_PDGs.push_back(443);
0147 m_PDGs.push_back(100443);
0148 m_PDGs.push_back(9000443);
0149 m_PDGs.push_back(9010443);
0150 m_PDGs.push_back(9020443);
0151 m_PDGs.push_back(10441);
0152 m_PDGs.push_back(20443);
0153 m_PDGs.push_back(445);
0154
0155 m_PDGs.push_back(30443);
0156 m_PDGs.push_back(551);
0157 m_PDGs.push_back(553);
0158 m_PDGs.push_back(100553);
0159 m_PDGs.push_back(200553);
0160 m_PDGs.push_back(10551);
0161 m_PDGs.push_back(20553);
0162 m_PDGs.push_back(555);
0163 m_PDGs.push_back(10553);
0164
0165 m_PDGs.push_back(110551);
0166 m_PDGs.push_back(120553);
0167 m_PDGs.push_back(100555);
0168 m_PDGs.push_back(210551);
0169 m_PDGs.push_back(220553);
0170 m_PDGs.push_back(200555);
0171 m_PDGs.push_back(30553);
0172 m_PDGs.push_back(20555);
0173
0174 m_PDGs.push_back(557);
0175 m_PDGs.push_back(130553);
0176 m_PDGs.push_back(120555);
0177 m_PDGs.push_back(100557);
0178 m_PDGs.push_back(110553);
0179 m_PDGs.push_back(210553);
0180 m_PDGs.push_back(10555);
0181 m_PDGs.push_back(110555);
0182
0183 m_PDGs.push_back(4122);
0184 m_PDGs.push_back(4132);
0185
0186 m_PDGs.push_back(4112);
0187 m_PDGs.push_back(4212);
0188 m_PDGs.push_back(4232);
0189 m_PDGs.push_back(4222);
0190 m_PDGs.push_back(4322);
0191 m_PDGs.push_back(4312);
0192
0193 m_PDGs.push_back(102132);
0194 m_PDGs.push_back(103124);
0195 m_PDGs.push_back(203122);
0196 m_PDGs.push_back(103122);
0197 m_PDGs.push_back(123122);
0198 m_PDGs.push_back(213122);
0199 m_PDGs.push_back(103126);
0200 m_PDGs.push_back(13212);
0201
0202
0203 m_PDGs.push_back(203126);
0204 m_PDGs.push_back(102134);
0205 m_PDGs.push_back(3122);
0206 m_PDGs.push_back(3222);
0207 m_PDGs.push_back(2214);
0208 m_PDGs.push_back(2224);
0209 m_PDGs.push_back(3324);
0210 m_PDGs.push_back(2114);
0211 m_PDGs.push_back(1114);
0212 m_PDGs.push_back(3112);
0213 m_PDGs.push_back(3212);
0214 m_PDGs.push_back(3114);
0215 m_PDGs.push_back(3224);
0216 m_PDGs.push_back(3214);
0217 m_PDGs.push_back(3216);
0218 m_PDGs.push_back(3322);
0219 m_PDGs.push_back(3312);
0220 m_PDGs.push_back(3314);
0221 m_PDGs.push_back(3334);
0222
0223 m_PDGs.push_back(4114);
0224 m_PDGs.push_back(4214);
0225 m_PDGs.push_back(4224);
0226 m_PDGs.push_back(4314);
0227 m_PDGs.push_back(4324);
0228 m_PDGs.push_back(4332);
0229 m_PDGs.push_back(4334);
0230
0231
0232 m_PDGs.push_back(10443);
0233
0234 m_PDGs.push_back(5122);
0235 m_PDGs.push_back(5132);
0236 m_PDGs.push_back(5232);
0237 m_PDGs.push_back(5332);
0238 m_PDGs.push_back(5222);
0239 m_PDGs.push_back(5112);
0240 m_PDGs.push_back(5212);
0241 m_PDGs.push_back(541);
0242 m_PDGs.push_back(14122);
0243 m_PDGs.push_back(14124);
0244 m_PDGs.push_back(5312);
0245 m_PDGs.push_back(5322);
0246 m_PDGs.push_back(10521);
0247 m_PDGs.push_back(20523);
0248 m_PDGs.push_back(10523);
0249
0250 m_PDGs.push_back(525);
0251 m_PDGs.push_back(10511);
0252 m_PDGs.push_back(20513);
0253 m_PDGs.push_back(10513);
0254 m_PDGs.push_back(515);
0255 m_PDGs.push_back(10531);
0256 m_PDGs.push_back(20533);
0257 m_PDGs.push_back(10533);
0258 m_PDGs.push_back(535);
0259 m_PDGs.push_back(543);
0260 m_PDGs.push_back(545);
0261 m_PDGs.push_back(5114);
0262 m_PDGs.push_back(5224);
0263 m_PDGs.push_back(5214);
0264 m_PDGs.push_back(5314);
0265 m_PDGs.push_back(5324);
0266 m_PDGs.push_back(5334);
0267 m_PDGs.push_back(10541);
0268 m_PDGs.push_back(10543);
0269 m_PDGs.push_back(20543);
0270
0271 m_PDGs.push_back(4424);
0272 m_PDGs.push_back(4422);
0273 m_PDGs.push_back(4414);
0274 m_PDGs.push_back(4412);
0275 m_PDGs.push_back(4432);
0276 m_PDGs.push_back(4434);
0277
0278 m_PDGs.push_back(130);
0279
0280 for (unsigned int i = 0; i < m_PDGs.size(); i++) {
0281 int pdt = HepPID::translatePythiatoPDT(m_PDGs.at(i));
0282 if (pdt != m_PDGs.at(i))
0283 m_PDGs.at(i) = pdt;
0284 }
0285 }
0286
0287 EvtGenInterface::~EvtGenInterface() {}
0288
0289 void EvtGenInterface::init() {
0290
0291 fSpecialSettings.push_back("Pythia8:ParticleDecays:mixB = off");
0292
0293
0294 edm::FileInPath decay_table(fPSet->getParameter<std::string>("decay_table"));
0295 edm::FileInPath pdt(fPSet->getParameter<edm::FileInPath>("particle_property_file"));
0296
0297 bool usePythia = fPSet->getUntrackedParameter<bool>("use_internal_pythia", true);
0298 bool useTauola = fPSet->getUntrackedParameter<bool>("use_internal_tauola", true);
0299 bool usePhotos = fPSet->getUntrackedParameter<bool>("use_internal_photos", true);
0300
0301
0302 bool convertPythiaCodes = fPSet->getUntrackedParameter<bool>(
0303 "convertPythiaCodes", true);
0304 char* tmp =
0305 std::getenv("PYTHIA8DATA");
0306
0307 if (tmp == nullptr) {
0308 edm::LogError("EvtGenInterface::~EvtGenInterface")
0309 << "EvtGenInterface::init() PYTHIA8DATA not defined. Terminating program ";
0310 exit(0);
0311 }
0312
0313 std::string pythiaDir(tmp);
0314 std::string photonType("gamma");
0315 bool useEvtGenRandom(true);
0316
0317
0318 EvtExternalGenList genList(convertPythiaCodes, pythiaDir, photonType, useEvtGenRandom);
0319 EvtAbsRadCorr* radCorrEngine = nullptr;
0320 if (usePhotos)
0321 radCorrEngine = genList.getPhotosModel();
0322 std::list<EvtDecayBase*> extraModels = genList.getListOfModels();
0323 std::list<EvtDecayBase*> myExtraModels;
0324 for (unsigned int i = 0; i < extraModels.size(); i++) {
0325 std::list<EvtDecayBase*>::iterator it = extraModels.begin();
0326 std::advance(it, i);
0327 TString name = (*it)->getName();
0328 if (name.Contains("PYTHIA") && usePythia)
0329 myExtraModels.push_back(*it);
0330 if (name.Contains("TAUOLA") && useTauola)
0331 myExtraModels.push_back(*it);
0332 }
0333
0334
0335
0336 EvtModelUserReg userList;
0337 std::list<EvtDecayBase*> userModels = userList.getUserModels();
0338 for (unsigned int i = 0; i < userModels.size(); i++) {
0339 std::list<EvtDecayBase*>::iterator it = userModels.begin();
0340 std::advance(it, i);
0341 TString name = (*it)->getName();
0342 edm::LogInfo("EvtGenInterface::~EvtGenInterface") << "Adding user model: " << name;
0343 myExtraModels.push_back(*it);
0344 }
0345
0346
0347 BmixingOption = fPSet->getUntrackedParameter<int>("B_Mixing", 1);
0348 if (BmixingOption != 0 && BmixingOption != 1) {
0349 throw cms::Exception("Configuration") << "EvtGenProducer requires B_Mixing to be 0 (coherent) or 1 (incoherent) \n"
0350 "Please fix this in your configuration.";
0351 }
0352
0353
0354
0355 m_EvtGen = new EvtGen(
0356 decay_table.fullPath().c_str(), pdt.fullPath().c_str(), the_engine, radCorrEngine, &myExtraModels, BmixingOption);
0357
0358
0359 if (fPSet->exists("user_decay_file")) {
0360 std::vector<std::string> user_decays = fPSet->getParameter<std::vector<std::string> >("user_decay_file");
0361 for (unsigned int i = 0; i < user_decays.size(); i++) {
0362 edm::FileInPath user_decay(user_decays.at(i));
0363 m_EvtGen->readUDecay(user_decay.fullPath().c_str());
0364 }
0365 }
0366
0367 if (fPSet->exists("user_decay_embedded")) {
0368 std::vector<std::string> user_decay_lines = fPSet->getParameter<std::vector<std::string> >("user_decay_embedded");
0369 char user_decay_tmp[] = "user_decay_tmpfileXXXXXX";
0370 int tmp_creation = mkstemp(user_decay_tmp);
0371 FILE* tmpf = std::fopen(user_decay_tmp, "w");
0372 if (!tmpf || (tmp_creation == -1)) {
0373 edm::LogError("EvtGenInterface::~EvtGenInterface")
0374 << "EvtGenInterface::init() fails when trying to open a temporary file for embedded user.dec. Terminating "
0375 "program ";
0376 exit(0);
0377 }
0378 for (unsigned int i = 0; i < user_decay_lines.size(); i++) {
0379 user_decay_lines.at(i) += "\n";
0380 std::fputs(user_decay_lines.at(i).c_str(), tmpf);
0381 }
0382 std::fclose(tmpf);
0383 m_EvtGen->readUDecay(user_decay_tmp);
0384 }
0385
0386
0387 if (fPSet->exists("operates_on_particles")) {
0388 std::vector<int> tmpPIDs = fPSet->getParameter<std::vector<int> >("operates_on_particles");
0389 m_PDGs.clear();
0390 bool goodinput = false;
0391 if (!tmpPIDs.empty()) {
0392 if (tmpPIDs.size() == 1 && tmpPIDs[0] == 0)
0393 goodinput = false;
0394 else
0395 goodinput = true;
0396 } else {
0397 goodinput = false;
0398 }
0399 if (goodinput)
0400 m_PDGs = tmpPIDs;
0401 else
0402 SetDefault_m_PDGs();
0403 } else
0404 SetDefault_m_PDGs();
0405
0406 for (unsigned int i = 0; i < m_PDGs.size(); i++) {
0407 edm::LogInfo("EvtGenInterface::~EvtGenInterface")
0408 << "EvtGenInterface::init() Particles to Operate on: " << m_PDGs[i];
0409 }
0410
0411
0412 polarize_ids = fPSet->getUntrackedParameter<std::vector<int> >("particles_to_polarize", std::vector<int>());
0413 polarize_pol = fPSet->getUntrackedParameter<std::vector<double> >("particle_polarizations", std::vector<double>());
0414 if (polarize_ids.size() != polarize_pol.size()) {
0415 throw cms::Exception("Configuration")
0416 << "EvtGenProducer requires that the particles_to_polarize and particle_polarization\n"
0417 "vectors be the same size. Please fix this in your configuration.";
0418 }
0419 for (unsigned int ndx = 0; ndx < polarize_ids.size(); ndx++) {
0420 if (polarize_pol[ndx] < -1. || polarize_pol[ndx] > 1.) {
0421 throw cms::Exception("Configuration")
0422 << "EvtGenProducer error: particle polarizations must be in the range -1 < P < 1";
0423 }
0424 polarizations.insert(std::pair<int, float>(polarize_ids[ndx], polarize_pol[ndx]));
0425 }
0426
0427
0428 if (fPSet->exists("list_forced_decays")) {
0429 std::vector<std::string> forced_names = fPSet->getParameter<std::vector<std::string> >("list_forced_decays");
0430 for (unsigned int i = 0; i < forced_names.size(); i++) {
0431 EvtId found = EvtPDL::getId(forced_names[i]);
0432 if (found.getId() == -1)
0433 throw cms::Exception("Configuration") << "name in part list for ignored decays not found: " << forced_names[i];
0434 if (found.getId() == found.getAlias())
0435 throw cms::Exception("Configuration") << "name of ignored decays is not an alias: " << forced_names[i];
0436 forced_id.push_back(found);
0437 forced_pdgids.push_back(EvtPDL::getStdHep(found));
0438 }
0439 }
0440 edm::LogInfo("EvtGenInterface::~EvtGenInterface")
0441 << "Number of Forced Paricles is: " << forced_pdgids.size() << std::endl;
0442 for (unsigned int j = 0; j < forced_id.size(); j++) {
0443 edm::LogInfo("EvtGenInterface::~EvtGenInterface")
0444 << "Forced Paricles are: " << forced_pdgids.at(j) << " " << forced_id.at(j) << std::endl;
0445 }
0446
0447 if (fPSet->exists("list_ignored_pdgids")) {
0448 ignore_pdgids = fPSet->getUntrackedParameter<std::vector<int> >("list_ignored_pdgids");
0449 }
0450
0451 return;
0452 }
0453
0454 HepMC::GenEvent* EvtGenInterface::decay(HepMC::GenEvent* evt) {
0455 if (the_engine->engine() == nullptr) {
0456 throw edm::Exception(edm::errors::LogicError)
0457 << "The EvtGen code attempted to use a random number engine while\n"
0458 << "the engine pointer was null in EvtGenInterface::decay. This might\n"
0459 << "mean that the code was modified to generate a random number outside\n"
0460 << "the event and beginLuminosityBlock methods, which is not allowed.\n";
0461 }
0462 CLHEP::RandFlat m_flat(*the_engine->engine(), 0., 1.);
0463
0464
0465 unsigned int nisforced = 0;
0466 std::vector<std::vector<HepMC::GenParticle*> > forcedparticles;
0467 for (unsigned int i = 0; i < forced_pdgids.size(); i++)
0468 forcedparticles.push_back(std::vector<HepMC::GenParticle*>());
0469
0470
0471 for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0472 if ((*p)->status() == 1) {
0473 int idHep = (*p)->pdg_id();
0474 bool isignore = false;
0475 for (unsigned int i = 0; i < ignore_pdgids.size(); i++) {
0476 if (idHep == ignore_pdgids[i])
0477 isignore = true;
0478 }
0479 if (!isignore) {
0480 bool isforced = false;
0481 for (unsigned int i = 0; i < forced_pdgids.size(); i++) {
0482 if (idHep == forced_pdgids[i]) {
0483 forcedparticles.at(i).push_back(*p);
0484 nisforced++;
0485 isforced = true;
0486 break;
0487 }
0488 }
0489 bool isDefaultEvtGen = false;
0490 for (unsigned int i = 0; i < m_PDGs.size(); i++) {
0491 if (abs(idHep) == abs(m_PDGs[i])) {
0492 isDefaultEvtGen = true;
0493 break;
0494 }
0495 }
0496 EvtId idEvt = EvtPDL::evtIdFromStdHep(idHep);
0497 int ipart = idEvt.getId();
0498 EvtDecayTable* evtDecayTable = EvtDecayTable::getInstance();
0499 if (!isforced && isDefaultEvtGen && ipart != -1 && evtDecayTable->getNMode(ipart) != 0) {
0500 addToHepMC(*p, idEvt, evt, true);
0501 }
0502 }
0503 }
0504 }
0505
0506
0507 unsigned int which = (unsigned int)(nisforced * flat());
0508 if (which == nisforced && nisforced > 0)
0509 which = nisforced - 1;
0510
0511 unsigned int idx = 0;
0512 for (unsigned int i = 0; i < forcedparticles.size(); i++) {
0513 for (unsigned int j = 0; j < forcedparticles.at(i).size(); j++) {
0514 EvtId idEvt = EvtPDL::evtIdFromStdHep(forcedparticles.at(i).at(j)->pdg_id());
0515 if (idx == which) {
0516 idEvt = forced_id[i];
0517 edm::LogInfo("EvtGenInterface::decay ")
0518 << EvtPDL::getStdHep(idEvt) << " will force to decay " << idx + 1 << " out of " << nisforced << std::endl;
0519 }
0520 bool decayed = false;
0521 while (!decayed) {
0522 decayed = addToHepMC(forcedparticles.at(i).at(j), idEvt, evt, false);
0523 }
0524 idx++;
0525 }
0526 }
0527
0528
0529 for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0530 if ((*p)->end_vertex() && (*p)->status() == 1)
0531 edm::LogWarning("EvtGenInterface::decay error: incorrect status!");
0532 if ((*p)->end_vertex() && (*p)->end_vertex()->particles_out_size() == 0)
0533 edm::LogWarning("EvtGenInterface::decay error: empty end vertex!");
0534 }
0535 return evt;
0536 }
0537
0538
0539 bool EvtGenInterface::addToHepMC(HepMC::GenParticle* partHep,
0540 const EvtId& idEvt,
0541 HepMC::GenEvent* theEvent,
0542 bool del_daug) {
0543
0544
0545 EvtVector4R pInit(
0546 partHep->momentum().e(), partHep->momentum().px(), partHep->momentum().py(), partHep->momentum().pz());
0547 EvtParticle* parent = EvtParticleFactory::particleFactory(idEvt, pInit);
0548 HepMC::FourVector posHep = (partHep->production_vertex())->position();
0549 EvtVector4R vInit(posHep.t(), posHep.x(), posHep.y(), posHep.z());
0550
0551 if (EvtPDL::getSpinType(idEvt) == EvtSpinType::DIRAC && polarizations.count(partHep->pdg_id()) > 0) {
0552 const HepMC::FourVector& momHep = partHep->momentum();
0553 EvtVector4R momEvt;
0554 momEvt.set(momHep.t(), momHep.x(), momHep.y(), momHep.z());
0555
0556
0557 float pol = polarizations.find(partHep->pdg_id())->second;
0558 GlobalVector pPart(momHep.x(), momHep.y(), momHep.z());
0559 GlobalVector zHat(0., 0., 1.);
0560 GlobalVector zCrossP = zHat.cross(pPart);
0561 GlobalVector polVec = pol * zCrossP.unit();
0562 EvtSpinDensity theSpinDensity;
0563 theSpinDensity.setDim(2);
0564 theSpinDensity.set(0, 0, EvtComplex(1. / 2. + polVec.z() / 2., 0.));
0565 theSpinDensity.set(0, 1, EvtComplex(polVec.x() / 2., -polVec.y() / 2.));
0566 theSpinDensity.set(1, 0, EvtComplex(polVec.x() / 2., polVec.y() / 2.));
0567 theSpinDensity.set(1, 1, EvtComplex(1. / 2. - polVec.z() / 2., 0.));
0568 parent->setSpinDensityForwardHelicityBasis(theSpinDensity);
0569 }
0570 if (parent) {
0571
0572 m_EvtGen->generateDecay(parent);
0573 if (del_daug)
0574 go_through_daughters(parent);
0575
0576
0577 EvtHepMCEvent evtHepMCEvent;
0578 evtHepMCEvent.constructEvent(parent, vInit);
0579 HepMC::GenEvent* evtGenHepMCTree = evtHepMCEvent.getEvent();
0580 parent->deleteTree();
0581
0582
0583 if (!evtGenHepMCTree->particles_empty())
0584 update_particles(partHep, theEvent, (*evtGenHepMCTree->particles_begin()));
0585
0586 } else {
0587
0588 return false;
0589 }
0590 return true;
0591 }
0592
0593
0594 void EvtGenInterface::update_particles(HepMC::GenParticle* partHep, HepMC::GenEvent* theEvent, HepMC::GenParticle* p) {
0595 if (p->end_vertex()) {
0596 if (!partHep->end_vertex()) {
0597 HepMC::GenVertex* vtx = new HepMC::GenVertex(p->end_vertex()->position());
0598 theEvent->add_vertex(vtx);
0599 vtx->add_particle_in(partHep);
0600 }
0601 if (p->end_vertex()->particles_out_size() != 0) {
0602 for (HepMC::GenVertex::particles_out_const_iterator d = p->end_vertex()->particles_out_const_begin();
0603 d != p->end_vertex()->particles_out_const_end();
0604 d++) {
0605 HepMC::GenParticle* daughter = new HepMC::GenParticle((*d)->momentum(), (*d)->pdg_id(), (*d)->status());
0606 daughter->suggest_barcode(theEvent->particles_size() + 1);
0607 partHep->end_vertex()->add_particle_out(daughter);
0608 if ((*d)->end_vertex())
0609 update_particles(daughter, theEvent, (*d));
0610 }
0611 partHep->set_status(p->status());
0612 }
0613 }
0614 }
0615
0616 void EvtGenInterface::setRandomEngine(CLHEP::HepRandomEngine* v) {
0617 the_engine->setRandomEngine(v);
0618 fRandomEngine = v;
0619 }
0620
0621 double EvtGenInterface::flat() {
0622 if (!fRandomEngine) {
0623 throw cms::Exception("LogicError")
0624 << "EvtGenInterface::flat: Attempt to generate random number when engine pointer is null\n"
0625 << "This might mean that the code was modified to generate a random number outside the\n"
0626 << "event and beginLuminosityBlock methods, which is not allowed.\n";
0627 }
0628 return fRandomEngine->flat();
0629 }
0630
0631 bool EvtGenInterface::findLastinChain(HepMC::GenParticle*& p) {
0632 if (p->end_vertex()) {
0633 if (p->end_vertex()->particles_out_size() != 0) {
0634 for (HepMC::GenVertex::particles_out_const_iterator d = p->end_vertex()->particles_out_const_begin();
0635 d != p->end_vertex()->particles_out_const_end();
0636 d++) {
0637 if (abs((*d)->pdg_id()) == abs(p->pdg_id())) {
0638 p = *d;
0639 findLastinChain(p);
0640 return false;
0641 }
0642 }
0643 }
0644 }
0645 return true;
0646 }
0647
0648 bool EvtGenInterface::hasnoDaughter(HepMC::GenParticle* p) {
0649 if (p->end_vertex()) {
0650 if (p->end_vertex()->particles_out_size() != 0) {
0651 return false;
0652 }
0653 }
0654 return true;
0655 }
0656
0657 void EvtGenInterface::go_through_daughters(EvtParticle* part) {
0658 int NDaug = part->getNDaug();
0659 if (NDaug) {
0660 EvtParticle* Daughter;
0661 for (int i = 0; i < NDaug; i++) {
0662 Daughter = part->getDaug(i);
0663 int idHep = Daughter->getPDGId();
0664 int found = 0;
0665 for (unsigned int j = 0; j < forced_pdgids.size(); j++) {
0666 if (idHep == forced_pdgids[j]) {
0667 found = 1;
0668 Daughter->deleteDaughters();
0669 break;
0670 }
0671 }
0672 if (!found)
0673 go_through_daughters(Daughter);
0674 }
0675 }
0676 }
0677
0678 DEFINE_EDM_PLUGIN(EvtGenFactory, gen::EvtGenInterface, "EvtGen130");