PhotosppInterface.cc

Go to the documentation of this file.

#include <iostream>
#include "GeneratorInterface/PhotosInterface/interface/PhotosppInterface.h"
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h"

#include "FWCore/ServiceRegistry/interface/Service.h"

#include "HepMC/GenEvent.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/HEPEVT_Wrapper.h"

using namespace gen;
using namespace edm;
using namespace std;

#include "Photos/Photos.h"
#include "Photos/PhotosHepMCEvent.h"

CLHEP::HepRandomEngine* PhotosppInterface::fRandomEngine = nullptr;

PhotosppInterface::PhotosppInterface( const edm::ParameterSet& pset)
  : fOnlyPDG(-1),
    fAvoidTauLeptonicDecays(false),
    fIsInitialized(false),
    fPSet(0)
{
  // add ability to keep brem from hadronizer and only modify specific channels 10/27/2014
  bool UseHadronizerQEDBrem=false;
  fPSet = new ParameterSet(pset);
  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
  for (unsigned int ip=0; ip<par.size(); ++ip ){
    std::string curSet = par[ip];
    // Physics settings
    if(curSet=="UseHadronizerQEDBrem") UseHadronizerQEDBrem=true;
  }
  if(!UseHadronizerQEDBrem)fSpecialSettings.push_back("QED-brem-off:all");
}

void PhotosppInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine){fRandomEngine=decayRandomEngine;}

void PhotosppInterface::configureOnlyFor( int ipdg ){
   fOnlyPDG = ipdg;
   fSpecialSettings.clear();
   return;
}

void PhotosppInterface::init(){
  if ( fIsInitialized ) return; // do init only once
  Photospp::Photos::initialize();
  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
  for (unsigned int ip=0; ip<par.size(); ++ip ){
    std::string curSet = par[ip];
    
    // Physics settings
    if(curSet=="maxWtInterference")                Photospp::Photos::maxWtInterference(fPSet->getParameter<double>(curSet));
    if(curSet=="setInfraredCutOff")                Photospp::Photos::setInfraredCutOff(fPSet->getParameter<double>(curSet));
    if(curSet=="setAlphaQED")                      Photospp::Photos::setAlphaQED(fPSet->getParameter<double>(curSet));
    if(curSet=="setInterference")                  Photospp::Photos::setInterference(fPSet->getParameter<bool>(curSet));
    if(curSet=="setDoubleBrem")                    Photospp::Photos::setDoubleBrem(fPSet->getParameter<bool>(curSet));
    if(curSet=="setQuatroBrem")                    Photospp::Photos::setQuatroBrem(fPSet->getParameter<bool>(curSet));
    if(curSet=="setExponentiation")                Photospp::Photos::setExponentiation(fPSet->getParameter<bool>(curSet));
    if(curSet=="setCorrectionWtForW")              Photospp::Photos::setCorrectionWtForW(fPSet->getParameter<bool>(curSet));
    if(curSet=="setMeCorrectionWtForScalar")       Photospp::Photos::setMeCorrectionWtForScalar(fPSet->getParameter<bool>(curSet));
    if(curSet=="setMeCorrectionWtForW")            Photospp::Photos::setMeCorrectionWtForW(fPSet->getParameter<bool>(curSet));
    if(curSet=="setMeCorrectionWtForZ")            Photospp::Photos::setMeCorrectionWtForZ(fPSet->getParameter<bool>(curSet));
    if(curSet=="initializeKinematicCorrections")   Photospp::Photos::initializeKinematicCorrections(fPSet->getParameter<int>(curSet));
    if(curSet=="forceMassFrom4Vector")             Photospp::Photos::forceMassFrom4Vector(fPSet->getParameter<bool>(curSet));
    if(curSet=="forceMassFromEventRecord")         Photospp::Photos::forceMassFromEventRecord(fPSet->getParameter<int>(curSet));
    if(curSet=="ignoreParticlesOfStatus")          Photospp::Photos::ignoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
    if(curSet=="deIgnoreParticlesOfStatus")        Photospp::Photos::deIgnoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
    if(curSet=="setMomentumConservationThreshold") Photospp::Photos::setMomentumConservationThreshold(fPSet->getParameter<double>(curSet));
    if(curSet=="suppressAll")                      if(fPSet->getParameter<bool>(curSet)==true)Photospp::Photos::suppressAll();
  
  // Now setup more complicated radiation/mass supression and forcing.
    if(curSet=="suppressBremForBranch"){
      edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
      std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
      for(unsigned int i=0;i<v.size();i++){
    std::string vs = v[i];
    std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
    if(vpar.size()==1)  Photospp::Photos::suppressBremForBranch(0,vpar[0]);
    if(vpar.size()==2)  Photospp::Photos::suppressBremForBranch(0/*vpar[0]*/,vpar[1]);
    if(vpar.size()==3)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2]);
    if(vpar.size()==4)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
    if(vpar.size()==5)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
    if(vpar.size()==6)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
    if(vpar.size()==7)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
    if(vpar.size()==8)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
    if(vpar.size()==9)  Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8]);
    if(vpar.size()==10) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8],vpar[9]);
      }
    }
    if(curSet=="suppressBremForDecay"){
      edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
      std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
      for(unsigned int i=0;i<v.size();i++){
    std::string vs = v[i];
    std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
        if(vpar.size()==1)  Photospp::Photos::suppressBremForDecay(0,vpar[0]);
        if(vpar.size()==2)  Photospp::Photos::suppressBremForDecay(0/*vpar[0]*/,vpar[1]);
        if(vpar.size()==3)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2]);
        if(vpar.size()==4)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
        if(vpar.size()==5)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
        if(vpar.size()==6)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
        if(vpar.size()==7)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
        if(vpar.size()==8)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
        if(vpar.size()==9)  Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8]);
        if(vpar.size()==10) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8],vpar[9]);
      }
    }

    if(curSet=="forceBremForBranch"){
      edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
      std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
      for(unsigned int i=0;i<v.size();i++){
    std::string vs = v[i];
    std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
    if(vpar.size()==1)  Photospp::Photos::forceBremForBranch(0,vpar[0]);
    if(vpar.size()==2)  Photospp::Photos::forceBremForBranch(0/*vpar[0]*/,vpar[1]);
    if(vpar.size()==3)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2]);
    if(vpar.size()==4)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
    if(vpar.size()==5)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
    if(vpar.size()==6)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
    if(vpar.size()==7)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
    if(vpar.size()==8)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
    if(vpar.size()==9)  Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8]);
    if(vpar.size()==10) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8],vpar[9]);
      }
    }
    if(curSet=="forceBremForDecay"){
      edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
      std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
      for(unsigned int i=0;i<v.size();i++){
    std::string vs = v[i];
    std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
        if(vpar.size()==1)  Photospp::Photos::forceBremForDecay(0,vpar[0]);
        if(vpar.size()==2)  Photospp::Photos::forceBremForDecay(0/*vpar[0]*/,vpar[1]);
        if(vpar.size()==3)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2]);
        if(vpar.size()==4)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
        if(vpar.size()==5)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
        if(vpar.size()==6)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
        if(vpar.size()==7)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
        if(vpar.size()==8)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
        if(vpar.size()==9)  Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8]);
        if(vpar.size()==10) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7],vpar[8],vpar[9]);
      }
    }

    if(curSet=="forceMass"){
      edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
      std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
      for(unsigned int i=0;i<v.size();i++){
    std::string vs = v[i];
    std::vector<double> vpar = cfg.getParameter<std::vector<double> >(vs);
    if(vpar.size()==2) Photospp::Photos::forceMass((int)vpar[0], vpar[1]);
      }
    }
  }
  // implement options set via c++ 
  if(fOnlyPDG>=0){
    Photospp::Photos::suppressAll();
    Photospp::Photos::forceBremForBranch(0,fOnlyPDG);
    Photospp::Photos::forceBremForBranch(0,-1*fOnlyPDG);
  }
  if(fAvoidTauLeptonicDecays){
    Photospp::Photos::suppressBremForDecay(3, 15, 16, 11, -12);
    Photospp::Photos::suppressBremForDecay(3, -15, -16, -11, 12);
    Photospp::Photos::suppressBremForDecay(3, 15, 16, 13, -14);
    Photospp::Photos::suppressBremForDecay(3, -15, -16, -13, -14);
  }
  Photospp::Photos::iniInfo();
  fIsInitialized = true; 
  return;
}

HepMC::GenEvent* PhotosppInterface::apply( HepMC::GenEvent* evt){
  Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
  if(!fIsInitialized) return evt;
  int NPartBefore = evt->particles_size();
  Photospp::PhotosHepMCEvent PhotosEvt(evt);
  PhotosEvt.process();
  //Fix the vertices and barcodes based on Julia Yarba's solution from TauolaInterface
  for (HepMC::GenEvent::vertex_const_iterator vtx=evt->vertices_begin(); vtx!=evt->vertices_end(); vtx++ ){
    std::vector<int> BCodes;
    BCodes.clear();
    if(*vtx){
      for(HepMC::GenVertex::particle_iterator pitr=(*vtx)->particles_begin(HepMC::children);pitr!=(*vtx)->particles_end(HepMC::children);++pitr){
    if((*pitr)->barcode()>10000){
      BCodes.push_back((*pitr)->barcode());
    }
      }
    }
    if(BCodes.size() > 0){
      for(size_t ibc=0; ibc<BCodes.size(); ibc++){
    HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
    int nbc = p1->barcode() - 10000 + NPartBefore;
    p1->suggest_barcode(nbc);
      }
    }
  }
  return evt;
}

double PhotosppInterface::flat(){
  if ( !fRandomEngine ) {
    throw cms::Exception("LogicError")
      << "PhotosppInterface::flat: Attempt to generate random number when engine pointer is null\n"
      << "This might mean that the code was modified to generate a random number outside the\n"
      << "event and beginLuminosityBlock methods, which is not allowed.\n";
  }
  return fRandomEngine->flat();
}

void PhotosppInterface::statistics(){Photospp::Photos::iniInfo();}

DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosppInterface, "Photospp356");