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  fSpecialSettings.push_back("QED-brem-off:all");
29  fPSet = new ParameterSet(pset);
30 }
31 
32 void PhotosppInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine){fRandomEngine=decayRandomEngine;}
33 
34 void PhotosppInterface::configureOnlyFor( int ipdg ){
35  fOnlyPDG = ipdg;
36  fSpecialSettings.clear();
37  return;
38 }
39 
40 void PhotosppInterface::init(){
41  if ( fIsInitialized ) return; // do init only once
42  Photospp::Photos::initialize();
43  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
44  for (unsigned int ip=0; ip<par.size(); ++ip ){
45  std::string curSet = par[ip];
46 
47  // Physics settings
48  if(curSet=="maxWtInterference") Photospp::Photos::maxWtInterference(fPSet->getParameter<double>(curSet));
49  if(curSet=="setInfraredCutOff") Photospp::Photos::setInfraredCutOff(fPSet->getParameter<double>(curSet));
50  if(curSet=="setAlphaQED") Photospp::Photos::setAlphaQED(fPSet->getParameter<double>(curSet));
51  if(curSet=="setInterference") Photospp::Photos::setInterference(fPSet->getParameter<bool>(curSet));
52  if(curSet=="setDoubleBrem") Photospp::Photos::setDoubleBrem(fPSet->getParameter<bool>(curSet));
53  if(curSet=="setQuatroBrem") Photospp::Photos::setQuatroBrem(fPSet->getParameter<bool>(curSet));
54  if(curSet=="setExponentiation") Photospp::Photos::setExponentiation(fPSet->getParameter<bool>(curSet));
55  if(curSet=="setCorrectionWtForW") Photospp::Photos::setCorrectionWtForW(fPSet->getParameter<bool>(curSet));
56  if(curSet=="setMeCorrectionWtForScalar") Photospp::Photos::setMeCorrectionWtForScalar(fPSet->getParameter<bool>(curSet));
57  if(curSet=="setMeCorrectionWtForW") Photospp::Photos::setMeCorrectionWtForW(fPSet->getParameter<bool>(curSet));
58  if(curSet=="setMeCorrectionWtForZ") Photospp::Photos::setMeCorrectionWtForZ(fPSet->getParameter<bool>(curSet));
59  if(curSet=="initializeKinematicCorrections") Photospp::Photos::initializeKinematicCorrections(fPSet->getParameter<int>(curSet));
60  if(curSet=="forceMassFrom4Vector") Photospp::Photos::forceMassFrom4Vector(fPSet->getParameter<bool>(curSet));
61  if(curSet=="forceMassFromEventRecord") Photospp::Photos::forceMassFromEventRecord(fPSet->getParameter<int>(curSet));
62  if(curSet=="ignoreParticlesOfStatus") Photospp::Photos::ignoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
63  if(curSet=="deIgnoreParticlesOfStatus") Photospp::Photos::deIgnoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
64  if(curSet=="setMomentumConservationThreshold") Photospp::Photos::setMomentumConservationThreshold(fPSet->getParameter<double>(curSet));
65  if(curSet=="suppressAll") if(fPSet->getParameter<bool>(curSet)==true)Photospp::Photos::suppressAll();
66 
67  // Now setup more complicated radiation/mass supression and forcing.
68  if(curSet=="suppressBremForBranch"){
69  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
70  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
71  for(unsigned int i=0;i<v.size();i++){
72  std::string vs = v[i];
73  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
74  if(vpar.size()==1) Photospp::Photos::suppressBremForBranch(0,vpar[0]);
75  if(vpar.size()==2) Photospp::Photos::suppressBremForBranch(0/*vpar[0]*/,vpar[1]);
76  if(vpar.size()==3) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2]);
77  if(vpar.size()==4) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
78  if(vpar.size()==5) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
79  if(vpar.size()==6) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
80  if(vpar.size()==7) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
81  if(vpar.size()==8) Photospp::Photos::suppressBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
82  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]);
83  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]);
84  }
85  }
86  if(curSet=="suppressBremForDecay"){
87  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
88  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
89  for(unsigned int i=0;i<v.size();i++){
90  std::string vs = v[i];
91  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
92  if(vpar.size()==1) Photospp::Photos::suppressBremForDecay(0,vpar[0]);
93  if(vpar.size()==2) Photospp::Photos::suppressBremForDecay(0/*vpar[0]*/,vpar[1]);
94  if(vpar.size()==3) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2]);
95  if(vpar.size()==4) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
96  if(vpar.size()==5) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
97  if(vpar.size()==6) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
98  if(vpar.size()==7) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
99  if(vpar.size()==8) Photospp::Photos::suppressBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
100  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]);
101  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]);
102  }
103  }
104 
105  if(curSet=="forceBremForBranch"){
106  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
107  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
108  for(unsigned int i=0;i<v.size();i++){
109  std::string vs = v[i];
110  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
111  if(vpar.size()==1) Photospp::Photos::forceBremForBranch(0,vpar[0]);
112  if(vpar.size()==2) Photospp::Photos::forceBremForBranch(0/*vpar[0]*/,vpar[1]);
113  if(vpar.size()==3) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2]);
114  if(vpar.size()==4) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3]);
115  if(vpar.size()==5) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
116  if(vpar.size()==6) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
117  if(vpar.size()==7) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
118  if(vpar.size()==8) Photospp::Photos::forceBremForBranch(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
119  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]);
120  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]);
121  }
122  }
123  if(curSet=="forceBremForDecay"){
124  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
125  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
126  for(unsigned int i=0;i<v.size();i++){
127  std::string vs = v[i];
128  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
129  if(vpar.size()==1) Photospp::Photos::forceBremForDecay(0,vpar[0]);
130  if(vpar.size()==2) Photospp::Photos::forceBremForDecay(0/*vpar[0]*/,vpar[1]);
131  if(vpar.size()==3) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2]);
132  if(vpar.size()==4) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3]);
133  if(vpar.size()==5) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4]);
134  if(vpar.size()==6) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5]);
135  if(vpar.size()==7) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6]);
136  if(vpar.size()==8) Photospp::Photos::forceBremForDecay(vpar[0],vpar[1],vpar[2],vpar[3],vpar[4],vpar[5],vpar[6],vpar[7]);
137  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]);
138  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]);
139  }
140  }
141 
142  if(curSet=="forceMass"){
143  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
144  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
145  for(unsigned int i=0;i<v.size();i++){
146  std::string vs = v[i];
147  std::vector<double> vpar = cfg.getParameter<std::vector<double> >(vs);
148  if(vpar.size()==2) Photospp::Photos::forceMass((int)vpar[0], vpar[1]);
149  }
150  }
151  }
152  // implement options set via c++
153  if(fOnlyPDG>=0){
154  Photospp::Photos::suppressAll();
155  Photospp::Photos::forceBremForBranch(0,fOnlyPDG);
156  Photospp::Photos::forceBremForBranch(0,-1*fOnlyPDG);
157  }
158  if(fAvoidTauLeptonicDecays){
159  Photospp::Photos::suppressBremForDecay(3, 15, 16, 11, -12);
160  Photospp::Photos::suppressBremForDecay(3, -15, -16, -11, 12);
161  Photospp::Photos::suppressBremForDecay(3, 15, 16, 13, -14);
162  Photospp::Photos::suppressBremForDecay(3, -15, -16, -13, -14);
163  }
164  Photospp::Photos::iniInfo();
165  fIsInitialized = true;
166  return;
167 }
168 
169 HepMC::GenEvent* PhotosppInterface::apply( HepMC::GenEvent* evt){
170  Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
171  if(!fIsInitialized) return evt;
172  int NPartBefore = evt->particles_size();
173  Photospp::PhotosHepMCEvent PhotosEvt(evt);
174  PhotosEvt.process();
175  //Fix the vertices and barcodes based on Julia Yarba's solution from TauolaInterface
176  for (HepMC::GenEvent::vertex_const_iterator vtx=evt->vertices_begin(); vtx!=evt->vertices_end(); vtx++ ){
177  std::vector<int> BCodes;
178  BCodes.clear();
179  if(*vtx){
180  for(HepMC::GenVertex::particle_iterator pitr=(*vtx)->particles_begin(HepMC::children);pitr!=(*vtx)->particles_end(HepMC::children);++pitr){
181  if((*pitr)->barcode()>10000){
182  BCodes.push_back((*pitr)->barcode());
183  }
184  }
185  }
186  if(BCodes.size() > 0){
187  for(size_t ibc=0; ibc<BCodes.size(); ibc++){
188  HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
189  int nbc = p1->barcode() - 10000 + NPartBefore;
190  p1->suggest_barcode(nbc);
191  }
192  }
193  }
194  return evt;
195 }
196 
197 double PhotosppInterface::flat(){
198  if ( !fRandomEngine ) {
199  throw cms::Exception("LogicError")
200  << "PhotosppInterface::flat: Attempt to generate random number when engine pointer is null\n"
201  << "This might mean that the code was modified to generate a random number outside the\n"
202  << "event and beginLuminosityBlock methods, which is not allowed.\n";
203  }
204  return fRandomEngine->flat();
205 }
206 
207 void PhotosppInterface::statistics(){Photospp::Photos::iniInfo();}
208 
209 DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosppInterface, "Photospp355");
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:32
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:169
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