Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:53 

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