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  if(curSet=="setPairEmission") Photospp::Photos::setPairEmission(fPSet->getParameter<bool>(curSet));
75  if(curSet=="setPhotonEmission") Photospp::Photos::setPhotonEmission(fPSet->getParameter<bool>(curSet));
76  if(curSet=="setStopAtCriticalError") Photospp::Photos::setStopAtCriticalError(fPSet->getParameter<bool>(curSet));
77 
78  // Now setup more complicated radiation/mass supression and forcing.
79  if(curSet=="suppressBremForBranch"){
80  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
81  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
82  for(unsigned int i=0;i<v.size();i++){
83  std::string vs = v[i];
84  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
85  if(vpar.size()==1) Photospp::Photos::suppressBremForBranch(0,vpar[0]);
86  if(vpar.size()==2) Photospp::Photos::suppressBremForBranch(0/*vpar[0]*/,vpar[1]);
87  if(vpar.size()==3) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2]);
88  if(vpar.size()==4) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
89  if(vpar.size()==5) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
90  if(vpar.size()==6) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
91  if(vpar.size()==7) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
92  if(vpar.size()==8) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
93  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]);
94  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]);
95  }
96  }
97  if(curSet=="suppressBremForDecay"){
98  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
99  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
100  for(unsigned int i=0;i<v.size();i++){
101  std::string vs = v[i];
102  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
103  if(vpar.size()==1) Photospp::Photos::suppressBremForDecay(0,vpar[0]);
104  if(vpar.size()==2) Photospp::Photos::suppressBremForDecay(0/*vpar[0]*/,vpar[1]);
105  if(vpar.size()==3) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2]);
106  if(vpar.size()==4) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
107  if(vpar.size()==5) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
108  if(vpar.size()==6) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
109  if(vpar.size()==7) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
110  if(vpar.size()==8) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
111  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]);
112  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]);
113  }
114  }
115 
116  if(curSet=="forceBremForBranch"){
117  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
118  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
119  for(unsigned int i=0;i<v.size();i++){
120  std::string vs = v[i];
121  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
122  if(vpar.size()==1) Photospp::Photos::forceBremForBranch(0,vpar[0]);
123  if(vpar.size()==2) Photospp::Photos::forceBremForBranch(0/*vpar[0]*/,vpar[1]);
124  if(vpar.size()==3) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2]);
125  if(vpar.size()==4) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
126  if(vpar.size()==5) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
127  if(vpar.size()==6) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
128  if(vpar.size()==7) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
129  if(vpar.size()==8) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
130  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]);
131  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]);
132  }
133  }
134  if(curSet=="forceBremForDecay"){
135  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
136  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
137  for(unsigned int i=0;i<v.size();i++){
138  std::string vs = v[i];
139  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
140  if(vpar.size()==1) Photospp::Photos::forceBremForDecay(0,vpar[0]);
141  if(vpar.size()==2) Photospp::Photos::forceBremForDecay(0/*vpar[0]*/,vpar[1]);
142  if(vpar.size()==3) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2]);
143  if(vpar.size()==4) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
144  if(vpar.size()==5) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
145  if(vpar.size()==6) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
146  if(vpar.size()==7) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
147  if(vpar.size()==8) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
148  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]);
149  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]);
150  }
151  }
152 
153  if(curSet=="forceMass"){
154  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
155  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
156  for(unsigned int i=0;i<v.size();i++){
157  std::string vs = v[i];
158  std::vector<double> vpar = cfg.getParameter<std::vector<double> >(vs);
159  if(vpar.size()==2) Photospp::Photos::forceMass((int)vpar[0], vpar[1]);
160  }
161  }
162  }
163  // implement options set via c++
164  if(fOnlyPDG>=0){
165  Photospp::Photos::suppressAll();
166  Photospp::Photos::forceBremForBranch(0,fOnlyPDG);
167  Photospp::Photos::forceBremForBranch(0,-1*fOnlyPDG);
168  }
169  if(fAvoidTauLeptonicDecays){
170  Photospp::Photos::suppressBremForDecay(3, 15, 16, 11, -12);
171  Photospp::Photos::suppressBremForDecay(3, -15, -16, -11, 12);
172  Photospp::Photos::suppressBremForDecay(3, 15, 16, 13, -14);
173  Photospp::Photos::suppressBremForDecay(3, -15, -16, -13, -14);
174  }
175  Photospp::Photos::iniInfo();
176  fIsInitialized = true;
177  return;
178 }
179 
180 HepMC::GenEvent* PhotosppInterface::apply( HepMC::GenEvent* evt){
181  Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
182  if(!fIsInitialized) return evt;
183  int NPartBefore = evt->particles_size();
184  Photospp::PhotosHepMCEvent PhotosEvt(evt);
185  PhotosEvt.process();
186  //Fix the vertices and barcodes based on Julia Yarba's solution from TauolaInterface
187  for (HepMC::GenEvent::vertex_const_iterator vtx=evt->vertices_begin(); vtx!=evt->vertices_end(); vtx++ ){
188  std::vector<int> BCodes;
189  BCodes.clear();
190  if(*vtx){
191  for(HepMC::GenVertex::particle_iterator pitr=(*vtx)->particles_begin(HepMC::children);pitr!=(*vtx)->particles_end(HepMC::children);++pitr){
192  if((*pitr)->barcode()>10000){
193  BCodes.push_back((*pitr)->barcode());
194  }
195  }
196  }
197  if(BCodes.size() > 0){
198  for(size_t ibc=0; ibc<BCodes.size(); ibc++){
199  HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
200  int nbc = p1->barcode() - 10000 + NPartBefore;
201  p1->suggest_barcode(nbc);
202  }
203  }
204  }
205  return evt;
206 }
207 
208 double PhotosppInterface::flat(){
209  if ( !fRandomEngine ) {
210  throw cms::Exception("LogicError")
211  << "PhotosppInterface::flat: Attempt to generate random number when engine pointer is null\n"
212  << "This might mean that the code was modified to generate a random number outside the\n"
213  << "event and beginLuminosityBlock methods, which is not allowed.\n";
214  }
215  return fRandomEngine->flat();
216 }
217 
218 void PhotosppInterface::statistics(){Photospp::Photos::iniInfo();}
219 
220 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
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:249
gen::PhotosppInterface::apply
HepMC::GenEvent * apply(HepMC::GenEvent *)
Definition: PhotosppInterface.cc:180
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