Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-21 01:39:49 

0001 #include <iostream>
0002 #include "GeneratorInterface/PhotosInterface/interface/PhotosppInterface.h"
0003 #include "FWCore/PluginManager/interface/ModuleDef.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h"
0006 
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 
0009 #include "HepMC/GenEvent.h"
0010 #include "HepMC/IO_HEPEVT.h"
0011 #include "HepMC/HEPEVT_Wrapper.h"
0012 
0013 using namespace gen;
0014 using namespace edm;
0015 using namespace std;
0016 
0017 #include "Photos/Photos.h"
0018 #include "Photos/PhotosHepMCEvent.h"
0019 #include "Photos/PhotosHepMC3Event.h"
0020 
0021 CLHEP::HepRandomEngine* PhotosppInterface::fRandomEngine = nullptr;
0022 
0023 PhotosppInterface::PhotosppInterface(const edm::ParameterSet& pset)
0024     : fOnlyPDG(-1), fAvoidTauLeptonicDecays(false), fIsInitialized(false), fPSet(nullptr) {
0025   // add ability to keep brem from hadronizer and only modify specific channels 10/27/2014
0026   bool UseHadronizerQEDBrem = false;
0027   fPSet = new ParameterSet(pset);
0028   std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
0029   for (unsigned int ip = 0; ip < par.size(); ++ip) {
0030     std::string curSet = par[ip];
0031     // Physics settings
0032     if (curSet == "UseHadronizerQEDBrem")
0033       UseHadronizerQEDBrem = true;
0034   }
0035   if (!UseHadronizerQEDBrem)
0036     fSpecialSettings.push_back("QED-brem-off:all");
0037 }
0038 
0039 void PhotosppInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine) {
0040   fRandomEngine = decayRandomEngine;
0041 }
0042 
0043 void PhotosppInterface::configureOnlyFor(int ipdg) {
0044   fOnlyPDG = ipdg;
0045   fSpecialSettings.clear();
0046   return;
0047 }
0048 
0049 void PhotosppInterface::init() {
0050   if (fIsInitialized)
0051     return;  // do init only once
0052   Photospp::Photos::initialize();
0053   Photospp::Photos::createHistoryEntries(true, 746);  // P-H-O
0054   std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
0055   for (unsigned int ip = 0; ip < par.size(); ++ip) {
0056     std::string curSet = par[ip];
0057 
0058     // Physics settings
0059     if (curSet == "maxWtInterference")
0060       Photospp::Photos::maxWtInterference(fPSet->getParameter<double>(curSet));
0061     if (curSet == "setInfraredCutOff")
0062       Photospp::Photos::setInfraredCutOff(fPSet->getParameter<double>(curSet));
0063     if (curSet == "setAlphaQED")
0064       Photospp::Photos::setAlphaQED(fPSet->getParameter<double>(curSet));
0065     if (curSet == "setInterference")
0066       Photospp::Photos::setInterference(fPSet->getParameter<bool>(curSet));
0067     if (curSet == "setDoubleBrem")
0068       Photospp::Photos::setDoubleBrem(fPSet->getParameter<bool>(curSet));
0069     if (curSet == "setQuatroBrem")
0070       Photospp::Photos::setQuatroBrem(fPSet->getParameter<bool>(curSet));
0071     if (curSet == "setExponentiation")
0072       Photospp::Photos::setExponentiation(fPSet->getParameter<bool>(curSet));
0073     if (curSet == "setCorrectionWtForW")
0074       Photospp::Photos::setCorrectionWtForW(fPSet->getParameter<bool>(curSet));
0075     if (curSet == "setMeCorrectionWtForScalar")
0076       Photospp::Photos::setMeCorrectionWtForScalar(fPSet->getParameter<bool>(curSet));
0077     if (curSet == "setMeCorrectionWtForW")
0078       Photospp::Photos::setMeCorrectionWtForW(fPSet->getParameter<bool>(curSet));
0079     if (curSet == "setMeCorrectionWtForZ")
0080       Photospp::Photos::setMeCorrectionWtForZ(fPSet->getParameter<bool>(curSet));
0081     if (curSet == "initializeKinematicCorrections")
0082       Photospp::Photos::initializeKinematicCorrections(fPSet->getParameter<int>(curSet));
0083     if (curSet == "forceMassFrom4Vector")
0084       Photospp::Photos::forceMassFrom4Vector(fPSet->getParameter<bool>(curSet));
0085     if (curSet == "forceMassFromEventRecord")
0086       Photospp::Photos::forceMassFromEventRecord(fPSet->getParameter<int>(curSet));
0087     if (curSet == "ignoreParticlesOfStatus")
0088       Photospp::Photos::ignoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
0089     if (curSet == "deIgnoreParticlesOfStatus")
0090       Photospp::Photos::deIgnoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
0091     if (curSet == "setMomentumConservationThreshold")
0092       Photospp::Photos::setMomentumConservationThreshold(fPSet->getParameter<double>(curSet));
0093     if (curSet == "suppressAll")
0094       if (fPSet->getParameter<bool>(curSet) == true)
0095         Photospp::Photos::suppressAll();
0096     if (curSet == "setPairEmission")
0097       Photospp::Photos::setPairEmission(fPSet->getParameter<bool>(curSet));
0098     if (curSet == "setPhotonEmission")
0099       Photospp::Photos::setPhotonEmission(fPSet->getParameter<bool>(curSet));
0100     if (curSet == "setStopAtCriticalError")
0101       Photospp::Photos::setStopAtCriticalError(fPSet->getParameter<bool>(curSet));
0102     if (curSet == "createHistoryEntries")
0103       Photospp::Photos::createHistoryEntries(fPSet->getParameter<bool>(curSet), 746);
0104 
0105     // Now setup more complicated radiation/mass supression and forcing.
0106     if (curSet == "suppressBremForBranch") {
0107       edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0108       std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
0109       for (unsigned int i = 0; i < v.size(); i++) {
0110         std::string vs = v[i];
0111         std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
0112         if (vpar.size() == 1)
0113           Photospp::Photos::suppressBremForBranch(0, vpar[0]);
0114         if (vpar.size() == 2)
0115           Photospp::Photos::suppressBremForBranch(0 /*vpar[0]*/, vpar[1]);
0116         if (vpar.size() == 3)
0117           Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2]);
0118         if (vpar.size() == 4)
0119           Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3]);
0120         if (vpar.size() == 5)
0121           Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
0122         if (vpar.size() == 6)
0123           Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
0124         if (vpar.size() == 7)
0125           Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
0126         if (vpar.size() == 8)
0127           Photospp::Photos::suppressBremForBranch(
0128               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
0129         if (vpar.size() == 9)
0130           Photospp::Photos::suppressBremForBranch(
0131               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
0132         if (vpar.size() == 10)
0133           Photospp::Photos::suppressBremForBranch(
0134               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
0135       }
0136     }
0137     if (curSet == "suppressBremForDecay") {
0138       edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0139       std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
0140       for (unsigned int i = 0; i < v.size(); i++) {
0141         std::string vs = v[i];
0142         std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
0143         if (vpar.size() == 1)
0144           Photospp::Photos::suppressBremForDecay(0, vpar[0]);
0145         if (vpar.size() == 2)
0146           Photospp::Photos::suppressBremForDecay(0 /*vpar[0]*/, vpar[1]);
0147         if (vpar.size() == 3)
0148           Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2]);
0149         if (vpar.size() == 4)
0150           Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3]);
0151         if (vpar.size() == 5)
0152           Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
0153         if (vpar.size() == 6)
0154           Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
0155         if (vpar.size() == 7)
0156           Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
0157         if (vpar.size() == 8)
0158           Photospp::Photos::suppressBremForDecay(
0159               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
0160         if (vpar.size() == 9)
0161           Photospp::Photos::suppressBremForDecay(
0162               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
0163         if (vpar.size() == 10)
0164           Photospp::Photos::suppressBremForDecay(
0165               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
0166       }
0167     }
0168 
0169     if (curSet == "forceBremForBranch") {
0170       edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0171       std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
0172       for (unsigned int i = 0; i < v.size(); i++) {
0173         std::string vs = v[i];
0174         std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
0175         if (vpar.size() == 1)
0176           Photospp::Photos::forceBremForBranch(0, vpar[0]);
0177         if (vpar.size() == 2)
0178           Photospp::Photos::forceBremForBranch(0 /*vpar[0]*/, vpar[1]);
0179         if (vpar.size() == 3)
0180           Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2]);
0181         if (vpar.size() == 4)
0182           Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3]);
0183         if (vpar.size() == 5)
0184           Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
0185         if (vpar.size() == 6)
0186           Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
0187         if (vpar.size() == 7)
0188           Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
0189         if (vpar.size() == 8)
0190           Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
0191         if (vpar.size() == 9)
0192           Photospp::Photos::forceBremForBranch(
0193               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
0194         if (vpar.size() == 10)
0195           Photospp::Photos::forceBremForBranch(
0196               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
0197       }
0198     }
0199     if (curSet == "forceBremForDecay") {
0200       edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0201       std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
0202       for (unsigned int i = 0; i < v.size(); i++) {
0203         std::string vs = v[i];
0204         std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
0205         if (vpar.size() == 1)
0206           Photospp::Photos::forceBremForDecay(0, vpar[0]);
0207         if (vpar.size() == 2)
0208           Photospp::Photos::forceBremForDecay(0 /*vpar[0]*/, vpar[1]);
0209         if (vpar.size() == 3)
0210           Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2]);
0211         if (vpar.size() == 4)
0212           Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3]);
0213         if (vpar.size() == 5)
0214           Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
0215         if (vpar.size() == 6)
0216           Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
0217         if (vpar.size() == 7)
0218           Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
0219         if (vpar.size() == 8)
0220           Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
0221         if (vpar.size() == 9)
0222           Photospp::Photos::forceBremForDecay(
0223               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
0224         if (vpar.size() == 10)
0225           Photospp::Photos::forceBremForDecay(
0226               vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
0227       }
0228     }
0229 
0230     if (curSet == "forceMass") {
0231       edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0232       std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
0233       for (unsigned int i = 0; i < v.size(); i++) {
0234         std::string vs = v[i];
0235         std::vector<double> vpar = cfg.getParameter<std::vector<double> >(vs);
0236         if (vpar.size() == 2)
0237           Photospp::Photos::forceMass((int)vpar[0], vpar[1]);
0238       }
0239     }
0240   }
0241   // implement options set via c++
0242   if (fOnlyPDG >= 0) {
0243     Photospp::Photos::suppressAll();
0244     Photospp::Photos::forceBremForBranch(0, fOnlyPDG);
0245     Photospp::Photos::forceBremForBranch(0, -1 * fOnlyPDG);
0246   }
0247   if (fAvoidTauLeptonicDecays) {
0248     Photospp::Photos::suppressBremForDecay(3, 15, 16, 11, -12);
0249     Photospp::Photos::suppressBremForDecay(3, -15, -16, -11, 12);
0250     Photospp::Photos::suppressBremForDecay(3, 15, 16, 13, -14);
0251     Photospp::Photos::suppressBremForDecay(3, -15, -16, -13, -14);
0252   }
0253   Photospp::Photos::iniInfo();
0254   fIsInitialized = true;
0255   return;
0256 }
0257 
0258 HepMC::GenEvent* PhotosppInterface::apply(HepMC::GenEvent* evt) {
0259   Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
0260   if (!fIsInitialized)
0261     return evt;
0262   int NPartBefore = evt->particles_size();
0263   Photospp::PhotosHepMCEvent PhotosEvt(evt);
0264   PhotosEvt.process();
0265   //Fix the vertices and barcodes based on Julia Yarba's solution from TauolaInterface
0266   for (HepMC::GenEvent::vertex_const_iterator vtx = evt->vertices_begin(); vtx != evt->vertices_end(); vtx++) {
0267     std::vector<int> BCodes;
0268     BCodes.clear();
0269     if (*vtx) {
0270       for (HepMC::GenVertex::particle_iterator pitr = (*vtx)->particles_begin(HepMC::children);
0271            pitr != (*vtx)->particles_end(HepMC::children);
0272            ++pitr) {
0273         if ((*pitr)->barcode() > 10000) {
0274           BCodes.push_back((*pitr)->barcode());
0275         }
0276       }
0277     }
0278     if (!BCodes.empty()) {
0279       for (size_t ibc = 0; ibc < BCodes.size(); ibc++) {
0280         HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
0281         int nbc = p1->barcode() - 10000 + NPartBefore;
0282         p1->suggest_barcode(nbc);
0283       }
0284     }
0285   }
0286   return evt;
0287 }
0288 
0289 HepMC3::GenEvent* PhotosppInterface::apply(HepMC3::GenEvent* evt) {
0290   Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
0291   if (!fIsInitialized)
0292     return evt;
0293   Photospp::PhotosHepMC3Event PhotosEvt(evt);
0294   PhotosEvt.process();
0295   return evt;
0296 }
0297 
0298 double PhotosppInterface::flat() {
0299   if (!fRandomEngine) {
0300     throw cms::Exception("LogicError")
0301         << "PhotosppInterface::flat: Attempt to generate random number when engine pointer is null\n"
0302         << "This might mean that the code was modified to generate a random number outside the\n"
0303         << "event and beginLuminosityBlock methods, which is not allowed.\n";
0304   }
0305   return fRandomEngine->flat();
0306 }
0307 
0308 void PhotosppInterface::statistics() { Photospp::Photos::iniInfo(); }
0309 
0310 DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosppInterface, "Photospp356");
[ 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!