CMS 3D CMS Logo
TestPythiaDecays.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Lukas/TestPythiaDecays
4 // Class: TestPythiaDecays
5 //
13 //
14 // Original Author: Lukas Vanelderen
15 // Created: Tue, 13 May 2014 09:50:05 GMT
16 //
17 //
18 
19 // system include files
20 #include <memory>
21 
22 // user include files
23 #include "FWCore/Framework/interface/Frameworkfwd.h"
24 #include "FWCore/Framework/interface/stream/EDAnalyzer.h"
25 
26 #include "FWCore/Framework/interface/Event.h"
27 #include "FWCore/Framework/interface/MakerMacros.h"
28 
29 #include "FWCore/ParameterSet/interface/ParameterSet.h"
30 
31 #include "SimDataFormats/Track/interface/SimTrack.h"
32 #include "SimDataFormats/Vertex/interface/SimVertex.h"
33 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
34 
35 // pythia
36 #include <Pythia8/Pythia.h>
37 
38 // root
39 #include "TH1D.h"
40 #include "TFile.h"
41 #include "TLorentzVector.h"
42 
43 #if PYTHIA_VERSION_INTEGER < 8304
44 typedef Pythia8::ParticleDataEntry* ParticleDataEntryPtr;
45 #else
46 typedef Pythia8::ParticleDataEntryPtr ParticleDataEntryPtr;
47 #endif
48 //
49 // class declaration
50 //
51 
52 class TestPythiaDecays : public edm::stream::EDAnalyzer<> {
53 public:
54  explicit TestPythiaDecays(const edm::ParameterSet&);
55  ~TestPythiaDecays() override;
56 
57  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
58 
59 private:
60  void analyze(const edm::Event&, const edm::EventSetup&) override;
61 
62  //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
63  //virtual void endRun(edm::Run const&, edm::EventSetup const&) override;
64  //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
65  //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
66 
67  // ----------member data ---------------------------
68  std::vector<int> pids;
69  std::map<int, TH1D*> h_mass;
70  std::map<int, TH1D*> h_p;
71  std::map<int, TH1D*> h_v;
72  std::map<int, TH1D*> h_mass_ref;
73  std::map<int, TH1D*> h_plt; // plt: proper life time
74  std::map<int, TH1D*> h_originVertexRho;
75  std::map<int, TH1D*> h_originVertexZ;
76  std::map<int, TH1D*> h_decayVertexRho;
77  std::map<int, TH1D*> h_decayVertexZ;
78  std::map<int, TH1D*> h_plt_ref; // plt: proper life time
79  std::map<int, TH1D*> h_br;
80  std::map<int, TH1D*> h_br_ref;
81 
82  std::map<int, std::vector<std::string> > knownDecayModes;
83 
84  Pythia8::Pythia* pythia;
85  std::string outputFile;
86 };
87 
88 //
89 // constants, enums and typedefs
90 //
91 
92 //
93 // static data member definitions
94 //
95 
96 //
97 // constructors and destructor
98 //
99 TestPythiaDecays::TestPythiaDecays(const edm::ParameterSet& iConfig) {
100  // output file
101  outputFile = iConfig.getParameter<std::string>("outputFile");
102 
103  // create pythia8 instance to access particle data
104  pythia = new Pythia8::Pythia();
105  pythia->init();
106  Pythia8::ParticleData pdt = pythia->particleData;
107 
108  // which particles will we study?
109  pids.push_back(15); // tau
110  pids.push_back(211); // pi+
111  pids.push_back(111); // pi0
112  pids.push_back(130); // K0L
113  pids.push_back(321); // K+
114  pids.push_back(323); // K*(392)
115  pids.push_back(411); // D+
116  pids.push_back(521); // B+
117 
118  // define histograms
119  for (size_t i = 0; i < pids.size(); ++i) {
120  int pid = abs(pids[i]);
121 
122  // get particle data
123  if (!pdt.isParticle(pid)) {
124  std::cout << "ERROR: BAD PARTICLE, pythia is not aware of pid " << pid << std::endl;
125  std::exit(1);
126  }
127  ParticleDataEntryPtr pd = pdt.particleDataEntryPtr(pid);
128 
129  // mass histograms
130  double m0 = pd->m0();
131  double w = pd->mWidth();
132  double mmin, mmax;
133  if (w == 0) {
134  mmin = m0 - m0 / 1000.;
135  mmax = m0 + m0 / 1000.;
136  } else {
137  mmin = m0 - 2 * w;
138  mmax = m0 + 2 * w;
139  }
140  std::stringstream strstr;
141  strstr << "mass_" << pid;
142  h_mass[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, mmin, mmax);
143  h_mass_ref[pid] = (TH1D*)(h_mass[pid]->Clone(strstr.str().c_str()));
144  h_mass_ref[pid]->SetTitle(h_mass_ref[pid]->GetName());
145  if (w == 0)
146  h_mass_ref[pid]->Fill(m0);
147  else {
148  for (int b = 1; b <= h_mass_ref[pid]->GetNbinsX(); ++b) {
149  double _val = h_mass_ref[pid]->GetBinCenter(b);
150  h_mass_ref[pid]->SetBinContent(b, TMath::BreitWigner(_val, m0, w));
151  }
152  }
153 
154  // p histogram
155  strstr.str("");
156  strstr << "p_" << pid;
157  h_p[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 20);
158 
159  // v histogram
160  strstr.str("");
161  strstr << "v_" << pid;
162  h_v[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 1.);
163 
164  // ctau histograms
165  double ctau0 = pd->tau0() / 10.;
166  strstr.str("");
167  strstr << "plt_" << pid;
168  h_plt[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, std::min(5. * ctau0, 500.));
169  h_plt_ref[pid] = (TH1D*)(h_plt[pid]->Clone(strstr.str().c_str()));
170  h_plt_ref[pid]->SetTitle(h_plt_ref[pid]->GetName());
171  for (int b = 1; b <= h_plt_ref[pid]->GetNbinsX(); ++b) {
172  double _val = h_plt_ref[pid]->GetBinCenter(b);
173  h_plt_ref[pid]->SetBinContent(b, TMath::Exp(-_val / ctau0)); //convert mm to cm
174  }
175 
176  // br histograms
177  strstr.str("");
178  strstr << "br_" << pid;
179  h_br[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 0, 0, 0);
180  h_br[pid]->SetCanExtend(TH1::kAllAxes);
181  h_br_ref[pid] = (TH1D*)(h_br[pid]->Clone(strstr.str().c_str()));
182  h_br_ref[pid]->SetTitle(h_br_ref[pid]->GetName());
183  knownDecayModes[pid] = std::vector<std::string>();
184  for (int d = 0; d < pd->sizeChannels(); ++d) {
185  Pythia8::DecayChannel& channel = pd->channel(d);
186  std::vector<int> prod;
187  for (int p = 0; p < channel.multiplicity(); ++p) {
188  int pId = abs(channel.product(p));
189  // from FastSimulation/Event/src/KineParticleFilter.cc
190  bool particleCut =
191  (pId > 10 && pId != 12 && pId != 14 && pId != 16 && pId != 18 && pId != 21 && (pId < 23 || pId > 40) &&
192  (pId < 81 || pId > 100) && pId != 2101 && pId != 3101 && pId != 3201 && pId != 1103 && pId != 2103 &&
193  pId != 2203 && pId != 3103 && pId != 3203 && pId != 3303);
194  if (particleCut)
195  prod.push_back(abs(channel.product(p)));
196  }
197  std::sort(prod.begin(), prod.end());
198  strstr.str("");
199  for (size_t p = 0; p < prod.size(); ++p) {
200  strstr << "_" << prod[p];
201  }
202  std::string label = strstr.str();
203  h_br[pid]->Fill(label.c_str(), 0.);
204  h_br_ref[pid]->Fill(label.c_str(), channel.bRatio());
205  h_br[pid]->SetEntries(0);
206  knownDecayModes[pid].push_back(label);
207  }
208 
209  // vertex plots
210  strstr.str("");
211  strstr << "originVertexRho_" << pid;
212  h_originVertexRho[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 200);
213  strstr.str("");
214  strstr << "originVertexZ_" << pid;
215  h_originVertexZ[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 400);
216  strstr.str("");
217  strstr << "decayVertexRho_" << pid;
218  h_decayVertexRho[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 200);
219  strstr.str("");
220  strstr << "decayVertexZ_" << pid;
221  h_decayVertexZ[pid] = new TH1D(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 400);
222  }
223 }
224 
225 TestPythiaDecays::~TestPythiaDecays() {
226  // do anything here that needs to be done at desctruction time
227  // (e.g. close files, deallocate resources etc.)
228  TFile* f = TFile::Open(outputFile.c_str(), "RECREATE");
229  f->cd();
230  f->mkdir("observed");
231  f->mkdir("prediction");
232  for (size_t i = 0; i < pids.size(); ++i) {
233  int pid = pids[i];
234  f->cd("observed");
235  h_mass[pid]->Write();
236  h_plt[pid]->Write();
237  h_br[pid]->Write();
238  h_originVertexZ[pid]->Write();
239  h_originVertexRho[pid]->Write();
240  h_decayVertexZ[pid]->Write();
241  h_decayVertexRho[pid]->Write();
242  h_p[pid]->Write();
243  h_v[pid]->Write();
244  f->cd("prediction");
245  h_mass_ref[pid]->Write();
246  h_plt_ref[pid]->Write();
247  h_br_ref[pid]->Write();
248  }
249  f->Close();
250 }
251 
252 //
253 // member functions
254 //
255 
256 // ------------ method called for each event ------------
257 void TestPythiaDecays::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
258  using namespace edm;
259 
260  Handle<std::vector<SimTrack> > simtracks;
261  iEvent.getByLabel("famosSimHits", simtracks);
262 
263  Handle<std::vector<SimVertex> > simvertices;
264  iEvent.getByLabel("famosSimHits", simvertices);
265 
266  // create maps
267 
268  // initialize
269  std::map<size_t, std::vector<size_t> > childMap; // child indices vs parent index
270  std::map<size_t, int> parentMap; // parent index vs child index
271  for (size_t j = 0; j < simtracks->size(); j++) {
272  childMap[j] = std::vector<size_t>();
273  parentMap[j] = -1;
274  }
275 
276  // do the mapping
277  for (size_t j = 0; j < simtracks->size(); j++) {
278  size_t childIndex = j;
279  const SimTrack& child = simtracks->at(childIndex);
280  if (child.noVertex())
281  continue;
282  const SimVertex& vertex = simvertices->at(child.vertIndex());
283  if (vertex.noParent())
284  continue;
285  size_t parentIndex = vertex.parentIndex();
286  childMap[parentIndex].push_back(childIndex);
287  parentMap[childIndex] = int(parentIndex);
288  }
289 
290  for (size_t j = 0; j < simtracks->size(); j++) {
291  const SimTrack& parent = simtracks->at(j);
292  int pid = abs(parent.type());
293  if (std::find(pids.begin(), pids.end(), pid) == pids.end())
294  continue;
295 
296  // fill mass hist
297  double mass = parent.momentum().M();
298  h_mass[pid]->Fill(mass);
299 
300  // fill p hist
301  h_p[pid]->Fill(parent.momentum().P());
302 
303  // fill vertex position hist
304  if (!parent.noVertex()) {
305  const SimVertex& originVertex = simvertices->at(parent.vertIndex());
306  h_originVertexRho[pid]->Fill(originVertex.position().Rho());
307  h_originVertexZ[pid]->Fill(std::fabs(originVertex.position().Z()));
308  }
309  if (!childMap[j].empty()) {
310  const SimTrack& child = simtracks->at(childMap[j][0]);
311  const SimVertex& decayVertex = simvertices->at(child.vertIndex());
312  h_decayVertexRho[pid]->Fill(decayVertex.position().Rho());
313  h_decayVertexZ[pid]->Fill(std::fabs(decayVertex.position().Z()));
314  }
315  }
316 
317  for (std::map<size_t, std::vector<size_t> >::iterator it = childMap.begin(); it != childMap.end(); ++it) {
318  // fill ctau hist
319  size_t parentIndex = it->first;
320  const SimTrack& parent = simtracks->at(parentIndex);
321  int pid = abs(parent.type());
322  std::vector<size_t>& childIndices = it->second;
323  if (childIndices.empty())
324  continue;
325 
326  if (std::find(pids.begin(), pids.end(), pid) == pids.end())
327  continue;
328 
329  const SimVertex& origin_vertex = simvertices->at(parent.vertIndex());
330  const SimTrack& child0 = simtracks->at(childIndices[0]);
331  const SimVertex& decay_vertex = simvertices->at(child0.vertIndex());
332 
333  TLorentzVector lv_origin_vertex(origin_vertex.position().X(),
334  origin_vertex.position().Y(),
335  origin_vertex.position().Z(),
336  origin_vertex.position().T());
337  TLorentzVector lv_decay_vertex(decay_vertex.position().X(),
338  decay_vertex.position().Y(),
339  decay_vertex.position().Z(),
340  decay_vertex.position().T());
341  TLorentzVector lv_dist = lv_decay_vertex - lv_origin_vertex;
342  TLorentzVector lv_parent(
343  parent.momentum().Px(), parent.momentum().Py(), parent.momentum().Pz(), parent.momentum().E());
344  TVector3 boost = lv_parent.BoostVector();
345  lv_dist.Boost(-boost);
346  h_v[pid]->Fill(boost.Mag());
347  double plt = lv_dist.T();
348  h_plt[pid]->Fill(plt);
349 
350  // fill br hist
351  std::vector<int> prod;
352  for (size_t d = 0; d < childIndices.size(); ++d) {
353  prod.push_back(abs(simtracks->at(childIndices[d]).type()));
354  }
355  std::sort(prod.begin(), prod.end());
356  std::stringstream strstr;
357  for (size_t p = 0; p < prod.size(); ++p) {
358  strstr << "_" << prod[p];
359  }
360  std::string label = strstr.str();
361  if (std::find(knownDecayModes[pid].begin(), knownDecayModes[pid].end(), label) == knownDecayModes[pid].end())
362  label = "u" + label;
363  h_br[pid]->Fill(label.c_str(), 1.);
364  h_br_ref[pid]->Fill(label.c_str(), 0.); // keep h_br and h_br_ref in sync
365  }
366 }
367 
368 // ------------ method called when starting to processes a run ------------
369 /*
370 void
371 TestPythiaDecays::beginRun(edm::Run const&, edm::EventSetup const&)
372 {
373 }
374 */
375 
376 // ------------ method called when ending the processing of a run ------------
377 /*
378 void
379 TestPythiaDecays::endRun(edm::Run const&, edm::EventSetup const&)
380 {
381 }
382 */
383 
384 // ------------ method called when starting to processes a luminosity block ------------
385 /*
386 void
387 TestPythiaDecays::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
388 {
389 }
390 */
391 
392 // ------------ method called when ending the processing of a luminosity block ------------
393 /*
394 void
395 TestPythiaDecays::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
396 {
397 }
398 */
399 
400 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
401 void TestPythiaDecays::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
402  //The following says we do not know what parameters are allowed so do no validation
403  // Please change this to state exactly what you do use, even if it is no parameters
404  edm::ParameterSetDescription desc;
405  desc.setUnknown();
406  descriptions.addDefault(desc);
407 }
408 
409 //define this as a plug-in
410 DEFINE_FWK_MODULE(TestPythiaDecays);
TestPythiaDecays::h_mass
std::map< int, TH1D * > h_mass
Definition: TestPythiaDecays.cc:69
mps_fire.i
i
Definition: mps_fire.py:428
SimVertex
Definition: SimVertex.h:5
TestPythiaDecays::pythia
Pythia8::Pythia * pythia
Definition: TestPythiaDecays.cc:84
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
EDAnalyzer.h
TestPythiaDecays::knownDecayModes
std::map< int, std::vector< std::string > > knownDecayModes
Definition: TestPythiaDecays.cc:82
min
T min(T a, T b)
Definition: MathUtil.h:58
edm
HLT enums.
Definition: AlignableModifier.h:19
gather_cfg.cout
cout
Definition: gather_cfg.py:144
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
TestPythiaDecays::h_decayVertexZ
std::map< int, TH1D * > h_decayVertexZ
Definition: TestPythiaDecays.cc:77
TestPythiaDecays::~TestPythiaDecays
~TestPythiaDecays() override
Definition: TestPythiaDecays.cc:225
TestPythiaDecays::h_v
std::map< int, TH1D * > h_v
Definition: TestPythiaDecays.cc:71
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
boost
Definition: CLHEP.h:16
edm::Handle
Definition: AssociativeIterator.h:50
TestPythiaDecays::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: TestPythiaDecays.cc:401
TestPythiaDecays::TestPythiaDecays
TestPythiaDecays(const edm::ParameterSet &)
Definition: TestPythiaDecays.cc:99
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:9
TestPythiaDecays::h_decayVertexRho
std::map< int, TH1D * > h_decayVertexRho
Definition: TestPythiaDecays.cc:76
TestPythiaDecays::outputFile
std::string outputFile
Definition: TestPythiaDecays.cc:85
MakerMacros.h
PSimHit.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
TestPythiaDecays::h_mass_ref
std::map< int, TH1D * > h_mass_ref
Definition: TestPythiaDecays.cc:72
w
const double w
Definition: UKUtility.cc:23
SimVertex.h
mps_fire.end
end
Definition: mps_fire.py:242
dumpMFGeometry_cfg.prod
prod
Definition: dumpMFGeometry_cfg.py:24
b
double b
Definition: hdecay.h:118
TestPythiaDecays::h_plt
std::map< int, TH1D * > h_plt
Definition: TestPythiaDecays.cc:73
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
bphysicsOniaDQM_cfi.vertex
vertex
Definition: bphysicsOniaDQM_cfi.py:7
edm::ParameterSet
Definition: ParameterSet.h:47
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
edm::stream::EDAnalyzer
Definition: EDAnalyzer.h:31
Event.h
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
createfilelist.int
int
Definition: createfilelist.py:10
iEvent
int iEvent
Definition: GenABIO.cc:224
TestPythiaDecays::h_originVertexRho
std::map< int, TH1D * > h_originVertexRho
Definition: TestPythiaDecays.cc:74
edm::EventSetup
Definition: EventSetup.h:58
ParticleDataEntryPtr
Pythia8::ParticleDataEntry * ParticleDataEntryPtr
Definition: TestPythiaDecays.cc:44
CoreSimVertex::position
const math::XYZTLorentzVectorD & position() const
Definition: CoreSimVertex.h:21
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TestPythiaDecays::h_br
std::map< int, TH1D * > h_br
Definition: TestPythiaDecays.cc:79
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
SimTrack
Definition: SimTrack.h:9
Frameworkfwd.h
relativeConstraints.empty
bool empty
Definition: relativeConstraints.py:46
TestPythiaDecays
Definition: TestPythiaDecays.cc:52
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
SimTrack::vertIndex
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
Definition: SimTrack.h:33
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
TestPythiaDecays::h_p
std::map< int, TH1D * > h_p
Definition: TestPythiaDecays.cc:70
ztail.d
d
Definition: ztail.py:151
TestPythiaDecays::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: TestPythiaDecays.cc:257
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
genParticles_cff.map
map
Definition: genParticles_cff.py:11
SimTrack.h
beamvalidation.exit
def exit(msg="")
Definition: beamvalidation.py:52
ParameterSet.h
TestPythiaDecays::h_originVertexZ
std::map< int, TH1D * > h_originVertexZ
Definition: TestPythiaDecays.cc:75
TestPythiaDecays::pids
std::vector< int > pids
Definition: TestPythiaDecays.cc:68
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
edm::Event
Definition: Event.h:73
child
Definition: simpleInheritance.h:11
edm::ConfigurationDescriptions::addDefault
void addDefault(ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:99
class-composition.parent
parent
Definition: class-composition.py:98
TestPythiaDecays::h_br_ref
std::map< int, TH1D * > h_br_ref
Definition: TestPythiaDecays.cc:80
label
const char * label
Definition: PFTauDecayModeTools.cc:11
TestPythiaDecays::h_plt_ref
std::map< int, TH1D * > h_plt_ref
Definition: TestPythiaDecays.cc:78
Generated for CMSSW Reference Manual by doxygen 1.8.16