CMS 3D CMS Logo
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), fAvoidTauLeptonicDecays(false), fIsInitialized(false), fPSet(nullptr) {
24  // add ability to keep brem from hadronizer and only modify specific channels 10/27/2014
25  bool UseHadronizerQEDBrem = false;
26  fPSet = new ParameterSet(pset);
27  std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
28  for (unsigned int ip = 0; ip < par.size(); ++ip) {
29  std::string curSet = par[ip];
30  // Physics settings
31  if (curSet == "UseHadronizerQEDBrem")
32  UseHadronizerQEDBrem = true;
33  }
34  if (!UseHadronizerQEDBrem)
35  fSpecialSettings.push_back("QED-brem-off:all");
36 }
37 
38 void PhotosppInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine) {
39  fRandomEngine = decayRandomEngine;
40 }
41 
42 void PhotosppInterface::configureOnlyFor(int ipdg) {
43  fOnlyPDG = ipdg;
44  fSpecialSettings.clear();
45  return;
46 }
47 
48 void PhotosppInterface::init() {
49  if (fIsInitialized)
50  return; // do init only once
51  Photospp::Photos::initialize();
52  std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
53  for (unsigned int ip = 0; ip < par.size(); ++ip) {
54  std::string curSet = par[ip];
55 
56  // Physics settings
57  if (curSet == "maxWtInterference")
58  Photospp::Photos::maxWtInterference(fPSet->getParameter<double>(curSet));
59  if (curSet == "setInfraredCutOff")
60  Photospp::Photos::setInfraredCutOff(fPSet->getParameter<double>(curSet));
61  if (curSet == "setAlphaQED")
62  Photospp::Photos::setAlphaQED(fPSet->getParameter<double>(curSet));
63  if (curSet == "setInterference")
64  Photospp::Photos::setInterference(fPSet->getParameter<bool>(curSet));
65  if (curSet == "setDoubleBrem")
66  Photospp::Photos::setDoubleBrem(fPSet->getParameter<bool>(curSet));
67  if (curSet == "setQuatroBrem")
68  Photospp::Photos::setQuatroBrem(fPSet->getParameter<bool>(curSet));
69  if (curSet == "setExponentiation")
70  Photospp::Photos::setExponentiation(fPSet->getParameter<bool>(curSet));
71  if (curSet == "setCorrectionWtForW")
72  Photospp::Photos::setCorrectionWtForW(fPSet->getParameter<bool>(curSet));
73  if (curSet == "setMeCorrectionWtForScalar")
74  Photospp::Photos::setMeCorrectionWtForScalar(fPSet->getParameter<bool>(curSet));
75  if (curSet == "setMeCorrectionWtForW")
76  Photospp::Photos::setMeCorrectionWtForW(fPSet->getParameter<bool>(curSet));
77  if (curSet == "setMeCorrectionWtForZ")
78  Photospp::Photos::setMeCorrectionWtForZ(fPSet->getParameter<bool>(curSet));
79  if (curSet == "initializeKinematicCorrections")
80  Photospp::Photos::initializeKinematicCorrections(fPSet->getParameter<int>(curSet));
81  if (curSet == "forceMassFrom4Vector")
82  Photospp::Photos::forceMassFrom4Vector(fPSet->getParameter<bool>(curSet));
83  if (curSet == "forceMassFromEventRecord")
84  Photospp::Photos::forceMassFromEventRecord(fPSet->getParameter<int>(curSet));
85  if (curSet == "ignoreParticlesOfStatus")
86  Photospp::Photos::ignoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
87  if (curSet == "deIgnoreParticlesOfStatus")
88  Photospp::Photos::deIgnoreParticlesOfStatus(fPSet->getParameter<int>(curSet));
89  if (curSet == "setMomentumConservationThreshold")
90  Photospp::Photos::setMomentumConservationThreshold(fPSet->getParameter<double>(curSet));
91  if (curSet == "suppressAll")
92  if (fPSet->getParameter<bool>(curSet) == true)
93  Photospp::Photos::suppressAll();
94  if (curSet == "setPairEmission")
95  Photospp::Photos::setPairEmission(fPSet->getParameter<bool>(curSet));
96  if (curSet == "setPhotonEmission")
97  Photospp::Photos::setPhotonEmission(fPSet->getParameter<bool>(curSet));
98  if (curSet == "setStopAtCriticalError")
99  Photospp::Photos::setStopAtCriticalError(fPSet->getParameter<bool>(curSet));
100 
101  // Now setup more complicated radiation/mass supression and forcing.
102  if (curSet == "suppressBremForBranch") {
103  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
104  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
105  for (unsigned int i = 0; i < v.size(); i++) {
106  std::string vs = v[i];
107  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
108  if (vpar.size() == 1)
109  Photospp::Photos::suppressBremForBranch(0, vpar[0]);
110  if (vpar.size() == 2)
111  Photospp::Photos::suppressBremForBranch(0 /*vpar[0]*/, vpar[1]);
112  if (vpar.size() == 3)
113  Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2]);
114  if (vpar.size() == 4)
115  Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3]);
116  if (vpar.size() == 5)
117  Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
118  if (vpar.size() == 6)
119  Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
120  if (vpar.size() == 7)
121  Photospp::Photos::suppressBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
122  if (vpar.size() == 8)
123  Photospp::Photos::suppressBremForBranch(
124  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
125  if (vpar.size() == 9)
126  Photospp::Photos::suppressBremForBranch(
127  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
128  if (vpar.size() == 10)
129  Photospp::Photos::suppressBremForBranch(
130  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
131  }
132  }
133  if (curSet == "suppressBremForDecay") {
134  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
135  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
136  for (unsigned int i = 0; i < v.size(); i++) {
137  std::string vs = v[i];
138  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
139  if (vpar.size() == 1)
140  Photospp::Photos::suppressBremForDecay(0, vpar[0]);
141  if (vpar.size() == 2)
142  Photospp::Photos::suppressBremForDecay(0 /*vpar[0]*/, vpar[1]);
143  if (vpar.size() == 3)
144  Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2]);
145  if (vpar.size() == 4)
146  Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3]);
147  if (vpar.size() == 5)
148  Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
149  if (vpar.size() == 6)
150  Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
151  if (vpar.size() == 7)
152  Photospp::Photos::suppressBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
153  if (vpar.size() == 8)
154  Photospp::Photos::suppressBremForDecay(
155  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
156  if (vpar.size() == 9)
157  Photospp::Photos::suppressBremForDecay(
158  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
159  if (vpar.size() == 10)
160  Photospp::Photos::suppressBremForDecay(
161  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
162  }
163  }
164 
165  if (curSet == "forceBremForBranch") {
166  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
167  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
168  for (unsigned int i = 0; i < v.size(); i++) {
169  std::string vs = v[i];
170  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
171  if (vpar.size() == 1)
172  Photospp::Photos::forceBremForBranch(0, vpar[0]);
173  if (vpar.size() == 2)
174  Photospp::Photos::forceBremForBranch(0 /*vpar[0]*/, vpar[1]);
175  if (vpar.size() == 3)
176  Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2]);
177  if (vpar.size() == 4)
178  Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3]);
179  if (vpar.size() == 5)
180  Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
181  if (vpar.size() == 6)
182  Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
183  if (vpar.size() == 7)
184  Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
185  if (vpar.size() == 8)
186  Photospp::Photos::forceBremForBranch(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
187  if (vpar.size() == 9)
188  Photospp::Photos::forceBremForBranch(
189  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
190  if (vpar.size() == 10)
191  Photospp::Photos::forceBremForBranch(
192  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
193  }
194  }
195  if (curSet == "forceBremForDecay") {
196  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
197  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
198  for (unsigned int i = 0; i < v.size(); i++) {
199  std::string vs = v[i];
200  std::vector<int> vpar = cfg.getParameter<std::vector<int> >(vs);
201  if (vpar.size() == 1)
202  Photospp::Photos::forceBremForDecay(0, vpar[0]);
203  if (vpar.size() == 2)
204  Photospp::Photos::forceBremForDecay(0 /*vpar[0]*/, vpar[1]);
205  if (vpar.size() == 3)
206  Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2]);
207  if (vpar.size() == 4)
208  Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3]);
209  if (vpar.size() == 5)
210  Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4]);
211  if (vpar.size() == 6)
212  Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5]);
213  if (vpar.size() == 7)
214  Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6]);
215  if (vpar.size() == 8)
216  Photospp::Photos::forceBremForDecay(vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7]);
217  if (vpar.size() == 9)
218  Photospp::Photos::forceBremForDecay(
219  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8]);
220  if (vpar.size() == 10)
221  Photospp::Photos::forceBremForDecay(
222  vpar[0], vpar[1], vpar[2], vpar[3], vpar[4], vpar[5], vpar[6], vpar[7], vpar[8], vpar[9]);
223  }
224  }
225 
226  if (curSet == "forceMass") {
227  edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
228  std::vector<std::string> v = cfg.getParameter<std::vector<std::string> >("parameterSets");
229  for (unsigned int i = 0; i < v.size(); i++) {
230  std::string vs = v[i];
231  std::vector<double> vpar = cfg.getParameter<std::vector<double> >(vs);
232  if (vpar.size() == 2)
233  Photospp::Photos::forceMass((int)vpar[0], vpar[1]);
234  }
235  }
236  }
237  // implement options set via c++
238  if (fOnlyPDG >= 0) {
239  Photospp::Photos::suppressAll();
240  Photospp::Photos::forceBremForBranch(0, fOnlyPDG);
241  Photospp::Photos::forceBremForBranch(0, -1 * fOnlyPDG);
242  }
243  if (fAvoidTauLeptonicDecays) {
244  Photospp::Photos::suppressBremForDecay(3, 15, 16, 11, -12);
245  Photospp::Photos::suppressBremForDecay(3, -15, -16, -11, 12);
246  Photospp::Photos::suppressBremForDecay(3, 15, 16, 13, -14);
247  Photospp::Photos::suppressBremForDecay(3, -15, -16, -13, -14);
248  }
249  Photospp::Photos::iniInfo();
250  fIsInitialized = true;
251  return;
252 }
253 
254 HepMC::GenEvent* PhotosppInterface::apply(HepMC::GenEvent* evt) {
255  Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
256  if (!fIsInitialized)
257  return evt;
258  int NPartBefore = evt->particles_size();
259  Photospp::PhotosHepMCEvent PhotosEvt(evt);
260  PhotosEvt.process();
261  //Fix the vertices and barcodes based on Julia Yarba's solution from TauolaInterface
262  for (HepMC::GenEvent::vertex_const_iterator vtx = evt->vertices_begin(); vtx != evt->vertices_end(); vtx++) {
263  std::vector<int> BCodes;
264  BCodes.clear();
265  if (*vtx) {
266  for (HepMC::GenVertex::particle_iterator pitr = (*vtx)->particles_begin(HepMC::children);
267  pitr != (*vtx)->particles_end(HepMC::children);
268  ++pitr) {
269  if ((*pitr)->barcode() > 10000) {
270  BCodes.push_back((*pitr)->barcode());
271  }
272  }
273  }
274  if (!BCodes.empty()) {
275  for (size_t ibc = 0; ibc < BCodes.size(); ibc++) {
276  HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
277  int nbc = p1->barcode() - 10000 + NPartBefore;
278  p1->suggest_barcode(nbc);
279  }
280  }
281  }
282  return evt;
283 }
284 
285 double PhotosppInterface::flat() {
286  if (!fRandomEngine) {
287  throw cms::Exception("LogicError")
288  << "PhotosppInterface::flat: Attempt to generate random number when engine pointer is null\n"
289  << "This might mean that the code was modified to generate a random number outside the\n"
290  << "event and beginLuminosityBlock methods, which is not allowed.\n";
291  }
292  return fRandomEngine->flat();
293 }
294 
295 void PhotosppInterface::statistics() { Photospp::Photos::iniInfo(); }
296 
297 DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosppInterface, "Photospp356");
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
initialize
static AlgebraicMatrix initialize()
Definition: BeamSpotTransientTrackingRecHit.cc:24
gen::PhotosppInterface::fRandomEngine
static CLHEP::HepRandomEngine * fRandomEngine
Definition: PhotosppInterface.h:44
mps_fire.i
i
Definition: mps_fire.py:341
gen::PhotosppInterface::setRandomEngine
void setRandomEngine(CLHEP::HepRandomEngine *decayRandomEngine) override
Definition: PhotosppInterface.cc:38
decayRandomEngine
CLHEP::HepRandomEngine * decayRandomEngine
Definition: ExternalDecayDriver.cc:23
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
gen::PhotosppInterface::PhotosppInterface
PhotosppInterface(const edm::ParameterSet &pset)
Definition: PhotosppInterface.cc:22
nullptr
#define nullptr
Definition: GCC11Compatibility.h:37
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
gen::PhotosppInterface::configureOnlyFor
void configureOnlyFor(int) override
Definition: PhotosppInterface.cc:42
gen::PhotosppInterface::fPSet
edm::ParameterSet * fPSet
Definition: PhotosppInterface.h:42
gen::v
double v[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:74
PhotosppPythia8TauSettings_cfi.suppressAll
suppressAll
Definition: PhotosppPythia8TauSettings_cfi.py:7
PhotosppPythia8TauSettings_cfi.setMomentumConservationThreshold
setMomentumConservationThreshold
Definition: PhotosppPythia8TauSettings_cfi.py:6
edm
HLT enums.
Definition: AlignableModifier.h:19
p1
double p1[4]
Definition: TauolaWrapper.h:89
DEFINE_EDM_PLUGIN
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: PluginFactory.h:124
badGlobalMuonTaggersAOD_cff.vtx
vtx
Definition: badGlobalMuonTaggersAOD_cff.py:5
PhotosppPythia8TauSettings_cfi.forceBremForDecay
forceBremForDecay
Definition: PhotosppPythia8TauSettings_cfi.py:8
gen::PhotosppInterface::apply
HepMC::GenEvent * apply(HepMC::GenEvent *) override
Definition: PhotosppInterface.cc:254
gen::PhotosInterfaceBase::fSpecialSettings
std::vector< std::string > fSpecialSettings
Definition: PhotosInterfaceBase.h:28
Generated for CMSSW Reference Manual by doxygen 1.8.11