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