CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PhotosppInterface.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include "GeneratorInterface/PhotosInterface/interface/PhotosppInterface.h"
3 #include "FWCore/PluginManager/interface/ModuleDef.h"
4 #include "FWCore/Framework/interface/MakerMacros.h"
5 #include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h"
6 
7 #include "FWCore/ServiceRegistry/interface/Service.h"
8 
9 #include "HepMC/GenEvent.h"
10 #include "HepMC/IO_HEPEVT.h"
11 #include "HepMC/HEPEVT_Wrapper.h"
12 
13 using namespace gen;
14 using namespace edm;
15 using namespace std;
16 
17 #include "Photos/Photos.h"
18 #include "Photos/PhotosHepMCEvent.h"
19 
20 CLHEP::HepRandomEngine* PhotosppInterface::fRandomEngine = nullptr;
21 
22 PhotosppInterface::PhotosppInterface( const edm::ParameterSet& pset)
23  : fOnlyPDG(-1),
24  fAvoidTauLeptonicDecays(false),
25  fIsInitialized(false),
26  fPSet(0)
27 {
28  // add ability to keep brem from hadronizer and only modify specific channels 10/27/2014
29  bool UseHadronizerQEDBrem=false;
30  fPSet = new ParameterSet(pset);
31  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
32  for (unsigned int ip=0; ip<par.size(); ++ip ){
33  std::string curSet = par[ip];
34  // Physics settings
35  if(curSet=="UseHadronizerQEDBrem") UseHadronizerQEDBrem=true;
36  }
37  if(!UseHadronizerQEDBrem)fSpecialSettings.push_back("QED-brem-off:all");
38 }
39 
40 void PhotosppInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine){fRandomEngine=decayRandomEngine;}
41 
42 void PhotosppInterface::configureOnlyFor( int ipdg ){
43  fOnlyPDG = ipdg;
44  fSpecialSettings.clear();
45  return;
46 }
47 
48 void PhotosppInterface::init(){
49  if ( fIsInitialized ) return; // do init only once
50  Photospp::Photos::initialize();
51  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
52  for (unsigned int ip=0; ip<par.size(); ++ip ){
53  std::string curSet = par[ip];
54 
55  // Physics settings
56  if(curSet=="maxWtInterference") Photospp::Photos::maxWtInterference(fPSet->getParameter<double>(curSet));
57  if(curSet=="setInfraredCutOff") Photospp::Photos::setInfraredCutOff(fPSet->getParameter<double>(curSet));
58  if(curSet=="setAlphaQED") Photospp::Photos::setAlphaQED(fPSet->getParameter<double>(curSet));
59  if(curSet=="setInterference") Photospp::Photos::setInterference(fPSet->getParameter<bool>(curSet));
60  if(curSet=="setDoubleBrem") Photospp::Photos::setDoubleBrem(fPSet->getParameter<bool>(curSet));
61  if(curSet=="setQuatroBrem") Photospp::Photos::setQuatroBrem(fPSet->getParameter<bool>(curSet));
62  if(curSet=="setExponentiation") Photospp::Photos::setExponentiation(fPSet->getParameter<bool>(curSet));
63  if(curSet=="setCorrectionWtForW") Photospp::Photos::setCorrectionWtForW(fPSet->getParameter<bool>(curSet));
64  if(curSet=="setMeCorrectionWtForScalar") Photospp::Photos::setMeCorrectionWtForScalar(fPSet->getParameter<bool>(curSet));
65  if(curSet=="setMeCorrectionWtForW") Photospp::Photos::setMeCorrectionWtForW(fPSet->getParameter<bool>(curSet));
66  if(curSet=="setMeCorrectionWtForZ") Photospp::Photos::setMeCorrectionWtForZ(fPSet->getParameter<bool>(curSet));
67  if(curSet=="initializeKinematicCorrections") Photospp::Photos::initializeKinematicCorrections(fPSet->getParameter<int>(curSet));
68  if(curSet=="forceMassFrom4Vector") Photospp::Photos::forceMassFrom4Vector(fPSet->getParameter<bool>(curSet));
69  if(curSet=="forceMassFromEventRecord") Photospp::Photos::forceMassFromEventRecord(fPSet->getParameter<int>(curSet));
70  if(curSet=="ignoreParticlesOfStatus") Photospp::Photos::ignoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
71  if(curSet=="deIgnoreParticlesOfStatus") Photospp::Photos::deIgnoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
72  if(curSet=="setMomentumConservationThreshold") Photospp::Photos::setMomentumConservationThreshold(fPSet->getParameter<double>(curSet));
73  if(curSet=="suppressAll") if(fPSet->getParameter<bool>(curSet)==true)Photospp::Photos::suppressAll();
74 
75  // Now setup more complicated radiation/mass supression and forcing.
76  if(curSet=="suppressBremForBranch"){
77  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
78  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
79  for(unsigned int i=0;i<v.size();i++){
80  std::string vs = v[i];
81  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
82  if(vpar.size()==1) Photospp::Photos::suppressBremForBranch(0,vpar[0]);
83  if(vpar.size()==2) Photospp::Photos::suppressBremForBranch(0/*vpar[0]*/,vpar[1]);
84  if(vpar.size()==3) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2]);
85  if(vpar.size()==4) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
86  if(vpar.size()==5) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
87  if(vpar.size()==6) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
88  if(vpar.size()==7) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
89  if(vpar.size()==8) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
90  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]);
91  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]);
92  }
93  }
94  if(curSet=="suppressBremForDecay"){
95  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
96  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
97  for(unsigned int i=0;i<v.size();i++){
98  std::string vs = v[i];
99  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
100  if(vpar.size()==1) Photospp::Photos::suppressBremForDecay(0,vpar[0]);
101  if(vpar.size()==2) Photospp::Photos::suppressBremForDecay(0/*vpar[0]*/,vpar[1]);
102  if(vpar.size()==3) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2]);
103  if(vpar.size()==4) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
104  if(vpar.size()==5) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
105  if(vpar.size()==6) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
106  if(vpar.size()==7) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
107  if(vpar.size()==8) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
108  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]);
109  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]);
110  }
111  }
112 
113  if(curSet=="forceBremForBranch"){
114  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
115  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
116  for(unsigned int i=0;i<v.size();i++){
117  std::string vs = v[i];
118  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
119  if(vpar.size()==1) Photospp::Photos::forceBremForBranch(0,vpar[0]);
120  if(vpar.size()==2) Photospp::Photos::forceBremForBranch(0/*vpar[0]*/,vpar[1]);
121  if(vpar.size()==3) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2]);
122  if(vpar.size()==4) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
123  if(vpar.size()==5) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
124  if(vpar.size()==6) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
125  if(vpar.size()==7) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
126  if(vpar.size()==8) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
127  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]);
128  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]);
129  }
130  }
131  if(curSet=="forceBremForDecay"){
132  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
133  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
134  for(unsigned int i=0;i<v.size();i++){
135  std::string vs = v[i];
136  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
137  if(vpar.size()==1) Photospp::Photos::forceBremForDecay(0,vpar[0]);
138  if(vpar.size()==2) Photospp::Photos::forceBremForDecay(0/*vpar[0]*/,vpar[1]);
139  if(vpar.size()==3) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2]);
140  if(vpar.size()==4) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
141  if(vpar.size()==5) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
142  if(vpar.size()==6) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
143  if(vpar.size()==7) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
144  if(vpar.size()==8) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
145  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]);
146  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]);
147  }
148  }
149 
150  if(curSet=="forceMass"){
151  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
152  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
153  for(unsigned int i=0;i<v.size();i++){
154  std::string vs = v[i];
155  std::vector<double> vpar = cfg.getParameter<std::vector<double> >(vs);
156  if(vpar.size()==2) Photospp::Photos::forceMass((int)vpar[0], vpar[1]);
157  }
158  }
159  }
160  // implement options set via c++
161  if(fOnlyPDG>=0){
162  Photospp::Photos::suppressAll();
163  Photospp::Photos::forceBremForBranch(0,fOnlyPDG);
164  Photospp::Photos::forceBremForBranch(0,-1*fOnlyPDG);
165  }
166  if(fAvoidTauLeptonicDecays){
167  Photospp::Photos::suppressBremForDecay(3, 15, 16, 11, -12);
168  Photospp::Photos::suppressBremForDecay(3, -15, -16, -11, 12);
169  Photospp::Photos::suppressBremForDecay(3, 15, 16, 13, -14);
170  Photospp::Photos::suppressBremForDecay(3, -15, -16, -13, -14);
171  }
172  Photospp::Photos::iniInfo();
173  fIsInitialized = true;
174  return;
175 }
176 
177 HepMC::GenEvent* PhotosppInterface::apply( HepMC::GenEvent* evt){
178  Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
179  if(!fIsInitialized) return evt;
180  int NPartBefore = evt->particles_size();
181  Photospp::PhotosHepMCEvent PhotosEvt(evt);
182  PhotosEvt.process();
183  //Fix the vertices and barcodes based on Julia Yarba's solution from TauolaInterface
184  for (HepMC::GenEvent::vertex_const_iterator vtx=evt->vertices_begin(); vtx!=evt->vertices_end(); vtx++ ){
185  std::vector<int> BCodes;
186  BCodes.clear();
187  if(*vtx){
188  for(HepMC::GenVertex::particle_iterator pitr=(*vtx)->particles_begin(HepMC::children);pitr!=(*vtx)->particles_end(HepMC::children);++pitr){
189  if((*pitr)->barcode()>10000){
190  BCodes.push_back((*pitr)->barcode());
191  }
192  }
193  }
194  if(BCodes.size() > 0){
195  for(size_t ibc=0; ibc<BCodes.size(); ibc++){
196  HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
197  int nbc = p1->barcode() - 10000 + NPartBefore;
198  p1->suggest_barcode(nbc);
199  }
200  }
201  }
202  return evt;
203 }
204 
205 double PhotosppInterface::flat(){
206  if ( !fRandomEngine ) {
207  throw cms::Exception("LogicError")
208  << "PhotosppInterface::flat: Attempt to generate random number when engine pointer is null\n"
209  << "This might mean that the code was modified to generate a random number outside the\n"
210  << "event and beginLuminosityBlock methods, which is not allowed.\n";
211  }
212  return fRandomEngine->flat();
213 }
214 
215 void PhotosppInterface::statistics(){Photospp::Photos::iniInfo();}
216 
217 DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosppInterface, "Photospp356");
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
initialize
static AlgebraicMatrix initialize()
Definition: BeamSpotTransientTrackingRecHit.cc:27
gen::PhotosppInterface::setRandomEngine
void setRandomEngine(CLHEP::HepRandomEngine *decayRandomEngine)
Definition: PhotosppInterface.cc:40
i
int i
Definition: DBlmapReader.cc:9
gen::PhotosppInterface::fRandomEngine
static CLHEP::HepRandomEngine * fRandomEngine
Definition: PhotosppInterface.h:44
looper.cfg
tuple cfg
Definition: looper.py:237
decayRandomEngine
CLHEP::HepRandomEngine * decayRandomEngine
Definition: ExternalDecayDriver.cc:23
gen::PhotosppInterface::PhotosppInterface
PhotosppInterface(const edm::ParameterSet &pset)
Definition: PhotosppInterface.cc:22
gen::PhotosppInterface::fPSet
edm::ParameterSet * fPSet
Definition: PhotosppInterface.h:42
gen::v
double v[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:80
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
relval_steps.gen
def gen
Definition: relval_steps.py:243
gen::PhotosppInterface::apply
HepMC::GenEvent * apply(HepMC::GenEvent *)
Definition: PhotosppInterface.cc:177
p1
double p1[4]
Definition: TauolaWrapper.h:89
DEFINE_EDM_PLUGIN
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: PluginFactory.h:101
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
gen::PhotosInterfaceBase::fSpecialSettings
std::vector< std::string > fSpecialSettings
Definition: PhotosInterfaceBase.h:27
Generated for CMSSW Reference Manual by doxygen 1.8.5