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