Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-28 03:54:33 

0001 // Class Based on EvtGenInterface(LHC).
0002 //
0003 // Created March 2014
0004 //
0005 // This class is a modification of the original EvtGenInterface which was developed for EvtGenLHC 9.1.
0006 // The modifications for EvtGen 1.3.0 are implemented by Ian M. Nugent
0007 // I would like to thank the EvtGen developers, in particular John Black, and Mikhail Kirsanov for their assistance.
0008 //
0009 // January 2015: Setting of coherent or incoherent B mixing included by Eduard Burelo
0010 // January 2015: Adding new feature to allow users to provide new evtgen models
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   // fill up default list of particles to be declared stable in the "master generator"
0062   // these are assumed to be PDG ID's
0063   //
0064   // Note: Pythia6's kc=43, 44, and 84 commented out because they're obsolete (per S.Mrenna)
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);  //B_s0
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);  // PDG ID = 9000111, Pythia6 ID = 10111
0131   m_PDGs.push_back(9000211);  // PDG ID = 9000211, Pythia6 ID = 10211
0132   m_PDGs.push_back(9010221);  // PDG ID = 9010211, Pythia6 ID = ???
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);  // PDG ID = 9020211, Pythia6 ID = ???
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   // m_PDGs.push_back( 84 ); // obsolete
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   //m_PDGs.push_back( 13241 ); unknown particle -typo?
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   //m_PDGs.push_back( 43 ); // obsolete (?)
0231   //m_PDGs.push_back( 44 ); // obsolete (?)
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   // flags for pythia8
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   //Setup evtGen following instructions on http://evtgen.warwick.ac.uk/docs/external/
0302   bool convertPythiaCodes = fPSet->getUntrackedParameter<bool>(
0303       "convertPythiaCodes", true);  // Specify if we want to use Pythia 6 physics codes for decays
0304   char* tmp =
0305       std::getenv("PYTHIA8DATA");  // Specify the pythia xml data directory to use the default PYTHIA8DATA location
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");  // Specify the photon type for Photos
0315   bool useEvtGenRandom(true);       // Specify if we want to use the EvtGen random number engine for these generators
0316 
0317   // Set up the default external generator list: Photos, Pythia and/or Tauola
0318   EvtExternalGenList genList(convertPythiaCodes, pythiaDir, photonType, useEvtGenRandom);
0319   EvtAbsRadCorr* radCorrEngine = nullptr;
0320   if (usePhotos)
0321     radCorrEngine = genList.getPhotosModel();                        // Get interface to radiative correction engine
0322   std::list<EvtDecayBase*> extraModels = genList.getListOfModels();  // get interface to Pythia and Tauola
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   //Set up user evtgen models
0335 
0336   EvtModelUserReg userList;
0337   std::list<EvtDecayBase*> userModels = userList.getUserModels();  // get interface to user models
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   // Set up the incoherent (1) or coherent (0) B mixing option
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   // Create the EvtGen generator object, passing the external generators
0355   m_EvtGen = new EvtGen(
0356       decay_table.fullPath().c_str(), pdt.fullPath().c_str(), the_engine, radCorrEngine, &myExtraModels, BmixingOption);
0357 
0358   // Add additional user information
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   // setup pdgid which the generator/hadronizer should not decay
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   // Obtain information to set polarization of particles
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   // Forced decays are particles that are aliased and forced to be decayed by EvtGen
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));  // force_pdgids is the list of stdhep codes
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   // Ignore decays are particles that are not to be decayed by EvtGen
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   // decay all request unforced particles and store the forced decays to later decay one per event
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   // notice this is a dynamic loop
0471   for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0472     if ((*p)->status() == 1) {  // all particles to be decays are set to status 1 by generator.hadronizer
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);  // storing and counting forced decays.
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);  // generate decay and remove daugther if they are forced
0501         }
0502       }
0503     }
0504   }
0505 
0506   // decay all forced particles (only 1/event is forced)... with no mixing allowed
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());  // "standard" decay Id
0515       if (idx == which) {
0516         idEvt = forced_id[i];  // force decay Id
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   // add code to ensure all particles have an end vertex and if they are undecayed with no end vertes set to status 1
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!");  //(*p)->set_status(2);
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 // Add particle to MC
0539 bool EvtGenInterface::addToHepMC(HepMC::GenParticle* partHep,
0540                                  const EvtId& idEvt,
0541                                  HepMC::GenEvent* theEvent,
0542                                  bool del_daug) {
0543   // Set up the parent particle from the HepMC GenEvent tree.
0544   //EvtVector4R pInit(EvtPDL::getMass(idEvt),partHep->momentum().px(),partHep->momentum().py(),partHep->momentum().pz());
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   // Reset polarization if requested....
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     // Particle is spin 1/2, so we can polarize it. Check polarizations map for particle, grab its polarization if it exists
0556     // and make the spin density matrix
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     // Generate the event
0572     m_EvtGen->generateDecay(parent);
0573     if (del_daug)
0574       go_through_daughters(parent);  // will delete all daugthers which are listed as forced, to allow decay them later
0575 
0576     // create HepMCTree
0577     EvtHepMCEvent evtHepMCEvent;
0578     evtHepMCEvent.constructEvent(parent, vInit);
0579     HepMC::GenEvent* evtGenHepMCTree = evtHepMCEvent.getEvent();
0580     parent->deleteTree();
0581 
0582     // update the event using a recursive function
0583     if (!evtGenHepMCTree->particles_empty())
0584       update_particles(partHep, theEvent, (*evtGenHepMCTree->particles_begin()));
0585 
0586   } else {
0587     // this should never hapend, in such case, speak out.
0588     return false;
0589   }
0590   return true;
0591 }
0592 
0593 //Recursivley add EvtGen decay to to Event Decy tree
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));  // if daugthers add them as well
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");
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!