CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
TauValidation.cc
Go to the documentation of this file.
1 /*class TauValidation
2  *
3  * Class to fill dqm monitor elements from existing EDM file
4  *
5  */
6 #include "Validation/EventGenerator/interface/TauValidation.h"
7 #include "CLHEP/Units/defs.h"
8 #include "CLHEP/Units/PhysicalConstants.h"
9 #include "DataFormats/Math/interface/LorentzVector.h"
10 #include "Validation/EventGenerator/interface/TauDecay_GenParticle.h"
11 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
12 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
13 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
14 #include <iostream>
15 #include "Validation/EventGenerator/interface/DQMHelper.h"
16 using namespace edm;
17 
18 TauValidation::TauValidation(const edm::ParameterSet &iPSet)
19  : // wmanager_(iPSet,consumesCollector())
20  genparticleCollection_(iPSet.getParameter<edm::InputTag>("genparticleCollection")),
21  NMODEID(TauDecay::NMODEID - 1), // fortran to C++ index
22  zsbins(20),
23  zsmin(-0.5),
24  zsmax(0.5) {
25  genparticleCollectionToken_ = consumes<reco::GenParticleCollection>(genparticleCollection_);
26  fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
27 }
28 
29 TauValidation::~TauValidation() {}
30 
31 void TauValidation::dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) {
32  fPDGTable = c.getHandle(fPDGTableToken);
33 }
34 
35 void TauValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) {
37  DQMHelper dqm(&i);
38  i.setCurrentFolder("Generator/Tau");
39  // Number of analyzed events
40  nTaus = dqm.book1dHisto("nTaus", "n analyzed Taus", 1, 0., 1., "bin", "Number of #tau's found");
41  nPrimeTaus =
42  dqm.book1dHisto("nPrimeTaus", "n analyzed prime Taus", 1, 0., 1., "bin", "Number of #tau's from Gauge Bosons");
43 
44  //Kinematics
45  TauPt = dqm.book1dHisto("TauPt", "Tau pT", 100, 0, 100, "P_{T}^{#tau}", "Number of #tau's from Gauge Bosons");
46  TauEta = dqm.book1dHisto("TauEta", "Tau eta", 100, -2.5, 2.5, "#eta^{#tau}", "Number of #tau's from Gauge Bosons");
47  TauPhi = dqm.book1dHisto("TauPhi", "Tau phi", 100, -3.14, 3.14, "#phi^{#tau}", "Number of #tau's from Gauge Bosons");
48  TauProngs = dqm.book1dHisto("TauProngs", "Tau n prongs", 7, 0, 7, "N_{prongs}", "Number of #tau's from Gauge Bosons");
49  TauDecayChannels = dqm.book1dHisto(
50  "TauDecayChannels", "Tau decay channels", 13, 0, 13, "Tau POG Decay Mode", "Number of #tau's from Gauge Bosons");
51  TauDecayChannels->setBinLabel(1 + undetermined, "?");
52  TauDecayChannels->setBinLabel(1 + electron, "e");
53  TauDecayChannels->setBinLabel(1 + muon, "mu");
54  TauDecayChannels->setBinLabel(1 + pi, "#pi^{#pm}");
55  TauDecayChannels->setBinLabel(1 + rho, "#rho^{#pm}");
56  TauDecayChannels->setBinLabel(1 + a1, "a_{1}^{#pm}");
57  TauDecayChannels->setBinLabel(1 + pi1pi0, "#pi^{#pm}#pi^{0}");
58  TauDecayChannels->setBinLabel(1 + pinpi0, "#pi^{#pm}n#pi^{0}");
59  TauDecayChannels->setBinLabel(1 + tripi, "3#pi^{#pm}");
60  TauDecayChannels->setBinLabel(1 + tripinpi0, "3#pi^{#pm}n#pi^{0}");
61  TauDecayChannels->setBinLabel(1 + K, "K");
62  TauDecayChannels->setBinLabel(1 + Kstar, "K^{*}");
63  TauDecayChannels->setBinLabel(1 + stable, "Stable");
64 
65  TauMothers = dqm.book1dHisto("TauMothers", "Tau mother particles", 10, 0, 10, "Mother of #tau", "Number of #tau's");
66 
67  TauMothers->setBinLabel(1 + other, "?");
68  TauMothers->setBinLabel(1 + B, "B Decays");
69  TauMothers->setBinLabel(1 + D, "D Decays");
70  TauMothers->setBinLabel(1 + gamma, "#gamma");
71  TauMothers->setBinLabel(1 + Z, "Z");
72  TauMothers->setBinLabel(1 + W, "W");
73  TauMothers->setBinLabel(1 + HSM, "H_{SM}/h^{0}");
74  TauMothers->setBinLabel(1 + H0, "H^{0}");
75  TauMothers->setBinLabel(1 + A0, "A^{0}");
76  TauMothers->setBinLabel(1 + Hpm, "H^{#pm}");
77 
78  DecayLength = dqm.book1dHisto(
79  "DecayLength", "#tau Decay Length", 100, -20, 20, "L_{#tau} (cm)", "Number of #tau's from Gauge Bosons");
80  LifeTime = dqm.book1dHisto(
81  "LifeTime", "#tau LifeTime ", 500, 0, 10000E-15, "#tau_{#tau} (s)", "Number of #tau's from Gauge Bosons");
82 
83  TauSpinEffectsW_X = dqm.book1dHisto(
84  "TauSpinEffectsWX", "X for pion", 50, 0, 1, "X", "Number of #tau#rightarrow#pi#nu from W^{#pm} Bosons");
85  TauSpinEffectsHpm_X = dqm.book1dHisto(
86  "TauSpinEffectsHpmX", "X for pion", 50, 0, 1, "X", "Number of #tau#rightarrow#pi#nu from H^{#pm} Bosons");
87 
88  TauSpinEffectsW_eX = dqm.book1dHisto(
89  "TauSpinEffectsWeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from W^{#pm} Bosons");
90  TauSpinEffectsHpm_eX = dqm.book1dHisto(
91  "TauSpinEffectsHpmeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from H^{#pm} Bosons");
92 
93  TauSpinEffectsW_muX = dqm.book1dHisto(
94  "TauSpinEffectsWmuX", "X for mu", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from W^{#pm} Bosons");
95  TauSpinEffectsHpm_muX = dqm.book1dHisto(
96  "TauSpinEffectsHpmmuX", "X for mue", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from H^{#pm} Bosons");
97 
98  TauSpinEffectsW_UpsilonRho = dqm.book1dHisto("TauSpinEffectsWUpsilonRho",
99  "#Upsilon for #rho",
100  50,
101  -1,
102  1,
103  "#Upsilon",
104  "Number of #tau#rightarrow#rho#nu from Gauge Bosons");
105  TauSpinEffectsHpm_UpsilonRho = dqm.book1dHisto("TauSpinEffectsHpmUpsilonRho",
106  "#Upsilon for #rho",
107  50,
108  -1,
109  1,
110  "#Upsilon",
111  "Number of #tau#rightarrow#rho#nu from Gauge Bosons");
112 
113  TauSpinEffectsW_UpsilonA1 = dqm.book1dHisto("TauSpinEffectsWUpsilonA1",
114  "#Upsilon for a1",
115  50,
116  -1,
117  1,
118  "#Upsilon",
119  "Number of #tau#rightarrow#pi#pi#pi#nu from Gauge Bosons");
120  TauSpinEffectsHpm_UpsilonA1 = dqm.book1dHisto("TauSpinEffectsHpmUpsilonA1",
121  "#Upsilon for a1",
122  50,
123  -1,
124  1,
125  "#Upsilon",
126  "Number of #tau#rightarrow#pi#pi#pi#nu from Gauge Bosons");
127 
128  TauSpinEffectsH_pipiAcoplanarity =
129  dqm.book1dHisto("TauSpinEffectsH_pipiAcoplanarity",
130  "H Acoplanarity for #pi^{-}#pi^{+}",
131  50,
132  0,
133  2 * TMath::Pi(),
134  "Acoplanarity",
135  "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
136 
137  TauSpinEffectsH_pipiAcollinearity =
138  dqm.book1dHisto("TauSpinEffectsH_pipiAcollinearity",
139  "H Acollinearity for #pi^{-}#pi^{+}",
140  50,
141  0,
142  TMath::Pi(),
143  "Acollinearity",
144  "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
145  TauSpinEffectsH_pipiAcollinearityzoom =
146  dqm.book1dHisto("TauSpinEffectsH_pipiAcollinearityzoom",
147  "H Acollinearity for #pi^{-}#pi^{+}",
148  50,
149  3,
150  TMath::Pi(),
151  "Acollinearity",
152  "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
153 
154  TauSpinEffectsZ_MVis =
155  dqm.book1dHisto("TauSpinEffectsZMVis",
156  "Mass of pi+ pi-",
157  25,
158  0,
159  1.1,
160  "M_{#pi^{+}#pi^{-}} (GeV)",
161  "Number of Z#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
162  TauSpinEffectsH_MVis =
163  dqm.book1dHisto("TauSpinEffectsHMVis",
164  "Mass of pi+ pi-",
165  25,
166  0,
167  1.1,
168  "M_{#pi^{+}#pi^{-}} (GeV)",
169  "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
170 
171  TauSpinEffectsZ_Zs =
172  dqm.book1dHisto("TauSpinEffectsZZs",
173  "Z_{s}",
174  zsbins,
175  zsmin,
176  zsmax,
177  "Z_{s}",
178  "Number of Z#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu Events");
179  TauSpinEffectsH_Zs =
180  dqm.book1dHisto("TauSpinEffectsHZs",
181  "Z_{s}",
182  zsbins,
183  zsmin,
184  zsmax,
185  "Z_{s}",
186  "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu Events");
187 
188  TauSpinEffectsZ_X = dqm.book1dHisto(
189  "TauSpinEffectsZX", "X for pion of #tau^{-}", 25, 0, 1.0, "X", "Number of #tau#rightarrow#pi#nu from Z Bosons");
190  TauSpinEffectsZ_X50to75 = dqm.book1dHisto("TauSpinEffectsZX50to75",
191  "X for pion of #tau^{-} (50GeV-75GeV)",
192  10,
193  0,
194  1.0,
195  "X",
196  "Number of #tau#rightarrow#pi#nu from Z(50GeV<M<75GeV) Bosons");
197  TauSpinEffectsZ_X75to88 = dqm.book1dHisto("TauSpinEffectsZX75to88",
198  "X for pion of #tau^{-} (75GeV-88GeV)",
199  10,
200  0,
201  1.0,
202  "X",
203  "Number of #tau#rightarrow#pi#nu from Z(75GeV<M<88GeV) Bosons");
204  TauSpinEffectsZ_X88to100 = dqm.book1dHisto("TauSpinEffectsZX88to100",
205  "X for pion of #tau^{-} (88GeV-100GeV)",
206  10,
207  0,
208  1.0,
209  "X",
210  "Number of #tau#rightarrow#pi#nu from Z(88GeV<M<100GeV) Bosons");
211  TauSpinEffectsZ_X100to120 = dqm.book1dHisto("TauSpinEffectsZX100to120",
212  "X for pion of #tau^{-} (100GeV-120GeV)",
213  10,
214  0,
215  1.0,
216  "X",
217  "Number of #tau#rightarrow#pi#nu from Z(100GeV<M<120GeV) Bosons");
218  TauSpinEffectsZ_X120UP = dqm.book1dHisto("TauSpinEffectsZX120UP",
219  "X for pion of #tau^{-} (>120GeV)",
220  10,
221  0,
222  1.0,
223  "X",
224  "Number of #tau#rightarrow#pi#nu from Z(120GeV<MGeV) Bosons");
225 
226  TauSpinEffectsH_X = dqm.book1dHisto(
227  "TauSpinEffectsH_X", "X for pion of #tau^{-}", 25, 0, 1.0, "X", "Number of #tau#rightarrow#pi#nu from H Bosons");
228 
229  TauSpinEffectsZ_Xf = dqm.book1dHisto("TauSpinEffectsZXf",
230  "X for pion of forward emitted #tau^{-}",
231  25,
232  0,
233  1.0,
234  "X_{f}",
235  "Number of #tau#rightarrow#pi#nu from Z Bosons");
236  TauSpinEffectsH_Xf = dqm.book1dHisto("TauSpinEffectsHXf",
237  "X for pion of forward emitted #tau^{-}",
238  25,
239  0,
240  1.0,
241  "X_{f}",
242  "Number of #tau#rightarrow#pi#nu from H Bosons");
243 
244  TauSpinEffectsZ_Xb = dqm.book1dHisto("TauSpinEffectsZXb",
245  "X for pion of backward emitted #tau^{-}",
246  25,
247  0,
248  1.0,
249  "X_{b}",
250  "Number of #tau#rightarrow#pi#nu from Z Bosons");
251  TauSpinEffectsH_Xb = dqm.book1dHisto("TauSpinEffectsHXb",
252  "X for pion of backward emitted #tau^{-}",
253  25,
254  0,
255  1.0,
256  "X_{b}",
257  "Number of #tau#rightarrow#pi#nu from H Bosons");
258 
259  TauSpinEffectsZ_eX = dqm.book1dHisto(
260  "TauSpinEffectsZeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from Gauge Bosons");
261  TauSpinEffectsH_eX = dqm.book1dHisto(
262  "TauSpinEffectsHeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from Gauge Bosons");
263 
264  TauSpinEffectsZ_muX = dqm.book1dHisto(
265  "TauSpinEffectsZmuX", "X for mu", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from Gauge Bosons");
266  TauSpinEffectsH_muX = dqm.book1dHisto(
267  "TauSpinEffectsHmuX", "X for mu", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from Gauge Bosons");
268 
269  TauSpinEffectsH_rhorhoAcoplanarityminus =
270  dqm.book1dHisto("TauSpinEffectsH_rhorhoAcoplanarityminus",
271  "#phi^{*-} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}<0)",
272  32,
273  0,
274  2 * TMath::Pi(),
275  "#phi^{*-} (Acoplanarity)",
276  "Number of H#rightarrow#tau^{-}(#rightarrow#rho^{-}#nu)#tau^{+}(#rightarrow#rho^{+}#nu) Events");
277  TauSpinEffectsH_rhorhoAcoplanarityplus =
278  dqm.book1dHisto("TauSpinEffectsH_rhorhoAcoplanarityplus",
279  "#phi^{*+} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}>0)",
280  32,
281  0,
282  2 * TMath::Pi(),
283  "#phi^{*+} (Acoplanarity)",
284  "Number of H#rightarrow#tau^{-}(#rightarrow#rho^{-}#nu)#tau^{+}(#rightarrow#rho^{+}#nu) Events");
285 
286  TauFSRPhotonsN = dqm.book1dHisto("TauFSRPhotonsN",
287  "FSR Photons radiating from/with tau (Gauge Boson)",
288  5,
289  -0.5,
290  4.5,
291  "N^{FSR Photons radiating from/with #tau}",
292  "Number of #tau's from Gauge Bosons");
293  TauFSRPhotonsPt = dqm.book1dHisto("TauFSRPhotonsPt",
294  "Pt of FSR Photons radiating from/with tau (Gauge Boson)",
295  100,
296  0,
297  100,
298  "P_{t}^{FSR Photons radiating from/with #tau [per #tau]} (GeV)",
299  "Number of #tau's from Gauge Bosons");
300  TauFSRPhotonsPtSum = dqm.book1dHisto("TauFSRPhotonsPtSum",
301  "Pt of FSR Photons radiating from/with tau (Gauge Boson)",
302  100,
303  0,
304  100,
305  "P_{t}^{FSR Photons radiating from/with #tau [per #tau]} (GeV)",
306  "Number of #tau's from Gauge Bosons");
307 
308  TauBremPhotonsN = dqm.book1dHisto("TauBremPhotonsN",
309  "Brem. Photons radiating in tau decay",
310  5,
311  -0.5,
312  4.5,
313  "N FSR Photons radiating from/with tau",
314  "Number of #tau's from Gauge Bosons");
315  TauBremPhotonsPt = dqm.book1dHisto("TauBremPhotonsPt",
316  "Sum Brem Pt ",
317  100,
318  0,
319  100,
320  "P_{t}^{Brem. Photons radiating in tau decay} (GeV)",
321  "Number of #tau's from Gauge Bosons");
322  TauBremPhotonsPtSum = dqm.book1dHisto("TauBremPhotonsPtSum",
323  "Sum of Brem Pt ",
324  100,
325  0,
326  100,
327  "Sum P_{t}^{Brem. Photons radiating in tau decay} (GeV)",
328  "Number of #tau's from Gauge Bosons");
329 
330  MODEID = dqm.book1dHisto("JAKID", "JAK ID", NMODEID + 1, -0.5, NMODEID + 0.5);
331  for (unsigned int j = 0; j < NMODEID + 1; j++) {
332  MODEInvMass.push_back(std::vector<MonitorElement *>());
333  std::string tmp = "JAKID";
334  tmp += std::to_string(j);
335  MODEInvMass.at(j).push_back(dqm.book1dHisto("M" + tmp,
336  "M_{" + TauDecay::DecayMode(j) + "} (GeV)",
337  80,
338  0,
339  2.0,
340  "M_{" + TauDecay::DecayMode(j) + "} (GeV)",
341  "Number of #tau's from Gauge Bosons"));
342  MODEID->setBinLabel(1 + j, TauDecay::DecayMode(j));
343  if (j == TauDecay::MODE_3PI || j == TauDecay::MODE_PI2PI0 || j == TauDecay::MODE_KPIK ||
344  j == TauDecay::MODE_KPIPI) {
345  MODEInvMass.at(j).push_back(dqm.book1dHisto("M13" + tmp,
346  "M_{13," + TauDecay::DecayMode(j) + "} (GeV)",
347  80,
348  0,
349  2.0,
350  "M_{13," + TauDecay::DecayMode(j) + "} (GeV)",
351  "Number of #tau's from Gauge Bosons"));
352  MODEInvMass.at(j).push_back(dqm.book1dHisto("M23" + tmp,
353  "M_{23," + TauDecay::DecayMode(j) + "} (GeV)",
354  80,
355  0,
356  2.0,
357  "M_{23," + TauDecay::DecayMode(j) + "} (GeV)",
358  "Number of #tau's from Gauge Bosons"));
359  MODEInvMass.at(j).push_back(dqm.book1dHisto("M12" + tmp,
360  "M_{12," + TauDecay::DecayMode(j) + "} (GeV)",
361  80,
362  0,
363  2.0,
364  "M_{12," + TauDecay::DecayMode(j) + "} (GeV)",
365  "Number of #tau's from Gauge Bosons"));
366  }
367  }
368  return;
369 }
370 
371 void TauValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
373  edm::Handle<reco::GenParticleCollection> genParticles;
374  iEvent.getByToken(genparticleCollectionToken_, genParticles);
375 
376  double weight = 1.0; //= wmanager_.weight(iEvent);
378  // find taus
379  for (reco::GenParticleCollection::const_iterator iter = genParticles->begin(); iter != genParticles->end(); ++iter) {
380  if (abs(iter->pdgId()) == PdtPdgMini::Z0 || abs(iter->pdgId()) == PdtPdgMini::Higgs0) {
381  spinEffectsZH(&(*iter), weight);
382  }
383  if (abs(iter->pdgId()) == 15) {
384  if (isLastTauinChain(&(*iter))) {
385  nTaus->Fill(0.5, weight);
386  int mother = tauMother(&(*iter), weight);
387  if (mother > -1) { // exclude B, D and other non-signal decay modes
388  nPrimeTaus->Fill(0.5, weight);
389  TauPt->Fill(iter->pt(), weight);
390  TauEta->Fill(iter->eta(), weight);
391  TauPhi->Fill(iter->phi(), weight);
392  photons(&(*iter), weight);
394  // Adding MODEID and Mass information
395  TauDecay_GenParticle TD;
396  unsigned int jak_id, TauBitMask;
397  if (TD.AnalyzeTau(&(*iter), jak_id, TauBitMask, false, false)) {
398  MODEID->Fill(jak_id, weight);
399  TauProngs->Fill(TD.nProng(TauBitMask), weight);
400  tauDecayChannel(&(*iter), jak_id, TauBitMask, weight);
401  if (jak_id <= NMODEID) {
402  int tcharge = iter->pdgId() / abs(iter->pdgId());
403  std::vector<const reco::GenParticle *> part = TD.Get_TauDecayProducts();
404  spinEffectsWHpm(&(*iter), mother, jak_id, part, weight);
405  TLorentzVector LVQ(0, 0, 0, 0);
406  TLorentzVector LVS12(0, 0, 0, 0);
407  TLorentzVector LVS13(0, 0, 0, 0);
408  TLorentzVector LVS23(0, 0, 0, 0);
409  bool haspart1 = false;
410  TVector3 PV, SV;
411  bool hasDL(false);
412  for (unsigned int i = 0; i < part.size(); i++) {
413  if (abs(part.at(i)->pdgId()) != PdtPdgMini::nu_tau && TD.isTauFinalStateParticle(part.at(i)->pdgId()) &&
414  !hasDL) {
415  hasDL = true;
416  TLorentzVector tlv(iter->px(), iter->py(), iter->pz(), iter->energy());
417  PV = TVector3(iter->vx(), iter->vy(), iter->vz());
418  SV = TVector3(part.at(i)->vx(), part.at(i)->vy(), part.at(i)->vz());
419  TVector3 DL = SV - PV;
420  DecayLength->Fill(DL.Dot(tlv.Vect()) / tlv.P(), weight);
421  double c(2.99792458E8), Ltau(DL.Mag() / 100) /*cm->m*/, beta(iter->p() / iter->mass());
422  LifeTime->Fill(Ltau / (c * beta), weight);
423  }
424 
425  if (TD.isTauFinalStateParticle(part.at(i)->pdgId()) && abs(part.at(i)->pdgId()) != PdtPdgMini::nu_e &&
426  abs(part.at(i)->pdgId()) != PdtPdgMini::nu_mu && abs(part.at(i)->pdgId()) != PdtPdgMini::nu_tau) {
427  TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
428  LVQ += LV;
429  if (jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0 ||
430  jak_id == TauDecay::MODE_KPIK || jak_id == TauDecay::MODE_KPIPI) {
431  if ((tcharge == part.at(i)->pdgId() / abs(part.at(i)->pdgId()) && TD.nProng(TauBitMask) == 3) ||
432  ((jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0) &&
433  TD.nProng(TauBitMask) == 1 && abs(part.at(i)->pdgId()) == PdtPdgMini::pi_plus)) {
434  LVS13 += LV;
435  LVS23 += LV;
436  } else {
437  LVS12 += LV;
438  if (!haspart1 && ((jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0) ||
439  ((jak_id != TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0) &&
440  abs(part.at(i)->pdgId()) == PdtPdgMini::K_plus))) {
441  LVS13 += LV;
442  haspart1 = true;
443  } else {
444  LVS23 += LV;
445  }
446  }
447  }
448  }
449  }
450  part.clear();
451  MODEInvMass.at(jak_id).at(0)->Fill(LVQ.M(), weight);
452  if (jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0 || jak_id == TauDecay::MODE_KPIK ||
453  jak_id == TauDecay::MODE_KPIPI) {
454  MODEInvMass.at(jak_id).at(1)->Fill(LVS13.M(), weight);
455  MODEInvMass.at(jak_id).at(2)->Fill(LVS23.M(), weight);
456  MODEInvMass.at(jak_id).at(3)->Fill(LVS12.M(), weight);
457  }
458  }
459  } else {
460  MODEID->Fill(jak_id, weight);
461  }
462  }
463  }
464  }
465  }
466 } //analyze
467 
468 const reco::GenParticle *TauValidation::GetMother(const reco::GenParticle *tau) {
469  for (unsigned int i = 0; i < tau->numberOfMothers(); i++) {
470  const reco::GenParticle *mother = static_cast<const reco::GenParticle *>(tau->mother(i));
471  if (mother->pdgId() == tau->pdgId())
472  return GetMother(mother);
473  return mother;
474  }
475  return tau;
476 }
477 
478 const std::vector<const reco::GenParticle *> TauValidation::GetMothers(const reco::GenParticle *boson) {
479  std::vector<const reco::GenParticle *> mothers;
480  for (unsigned int i = 0; i < boson->numberOfMothers(); i++) {
481  const reco::GenParticle *mother = static_cast<const reco::GenParticle *>(boson->mother(i));
482  if (mother->pdgId() == boson->pdgId())
483  return GetMothers(mother);
484  mothers.push_back(mother);
485  }
486  return mothers;
487 }
488 
489 int TauValidation::findMother(const reco::GenParticle *tau) { return TauValidation::GetMother(tau)->pdgId(); }
490 
491 bool TauValidation::isLastTauinChain(const reco::GenParticle *tau) {
492  for (unsigned int i = 0; i < tau->numberOfDaughters(); i++) {
493  if (tau->daughter(i)->pdgId() == tau->pdgId())
494  return false;
495  }
496  return true;
497 }
498 
499 void TauValidation::findTauList(const reco::GenParticle *tau, std::vector<const reco::GenParticle *> &TauList) {
500  TauList.insert(TauList.begin(), tau);
501  for (unsigned int i = 0; i < tau->numberOfMothers(); i++) {
502  const reco::GenParticle *mother = static_cast<const reco::GenParticle *>(tau->mother(i));
503  if (mother->pdgId() == tau->pdgId()) {
504  findTauList(mother, TauList);
505  }
506  }
507 }
508 
509 void TauValidation::findFSRandBrem(const reco::GenParticle *p,
510  bool doBrem,
511  std::vector<const reco::GenParticle *> &ListofFSR,
512  std::vector<const reco::GenParticle *> &ListofBrem) {
513  // note this code split the FSR and Brem based one if the tau decays into a tau+photon or not with the Fortran Tauola Interface, this is not 100% correct because photos puts the tau with the regular tau decay products.
514  if (abs(p->pdgId()) == 15) {
515  if (isLastTauinChain(p)) {
516  doBrem = true;
517  } else {
518  doBrem = false;
519  }
520  }
521  int photo_ID = 22;
522  for (unsigned int i = 0; i < p->numberOfDaughters(); i++) {
523  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(p->daughter(i));
524  if (abs((dau)->pdgId()) == abs(photo_ID) && !doBrem) {
525  ListofFSR.push_back(dau);
526  }
527  if (abs((dau)->pdgId()) == abs(photo_ID) && doBrem) {
528  ListofBrem.push_back(dau);
529  }
530  if (abs((dau)->pdgId()) != 111 && abs((dau)->pdgId()) != 221) { // remove pi0 and eta decays
531  findFSRandBrem(dau, doBrem, ListofFSR, ListofBrem);
532  }
533  }
534 }
535 
536 void TauValidation::FindPhotosFSR(const reco::GenParticle *p,
537  std::vector<const reco::GenParticle *> &ListofFSR,
538  double &BosonScale) {
539  BosonScale = 0.0;
540  const reco::GenParticle *m = GetMother(p);
541  int mother_pid = m->pdgId();
542  if (m->pdgId() != p->pdgId()) {
543  for (unsigned int i = 0; i < m->numberOfDaughters(); i++) {
544  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(m->daughter(i));
545  if (abs(dau->pdgId()) == 22) {
546  ListofFSR.push_back(dau);
547  }
548  }
549  }
550  if (abs(mother_pid) == 24)
551  BosonScale = 1.0; // W
552  if (abs(mother_pid) == 23)
553  BosonScale = 2.0; // Z;
554  if (abs(mother_pid) == 22)
555  BosonScale = 2.0; // gamma;
556  if (abs(mother_pid) == 25)
557  BosonScale = 2.0; // HSM;
558  if (abs(mother_pid) == 35)
559  BosonScale = 2.0; // H0;
560  if (abs(mother_pid) == 36)
561  BosonScale = 2.0; // A0;
562  if (abs(mother_pid) == 37)
563  BosonScale = 1.0; //Hpm;
564 }
565 
566 int TauValidation::tauMother(const reco::GenParticle *tau, double weight) {
567  if (abs(tau->pdgId()) != 15)
568  return -3;
569  int mother_pid = findMother(tau);
570  if (mother_pid == -2)
571  return -2;
572  int label = other;
573  if (abs(mother_pid) == 24)
574  label = W;
575  if (abs(mother_pid) == 23)
576  label = Z;
577  if (abs(mother_pid) == 22)
578  label = gamma;
579  if (abs(mother_pid) == 25)
580  label = HSM;
581  if (abs(mother_pid) == 35)
582  label = H0;
583  if (abs(mother_pid) == 36)
584  label = A0;
585  if (abs(mother_pid) == 37)
586  label = Hpm;
587  int mother_shortpid = (abs(mother_pid) % 10000);
588  if (mother_shortpid > 500 && mother_shortpid < 600)
589  label = B;
590  if (mother_shortpid > 400 && mother_shortpid < 500)
591  label = D;
592  TauMothers->Fill(label, weight);
593  if (label == B || label == D || label == other)
594  return -1;
595  return mother_pid;
596 }
597 
598 int TauValidation::tauDecayChannel(const reco::GenParticle *tau, int jak_id, unsigned int TauBitMask, double weight) {
599  int channel = undetermined;
600  if (tau->status() == 1)
601  channel = stable;
602  int allCount = 0, eCount = 0, muCount = 0, pi0Count = 0, piCount = 0, rhoCount = 0, a1Count = 0, KCount = 0,
603  KstarCount = 0;
604 
605  countParticles(tau, allCount, eCount, muCount, pi0Count, piCount, rhoCount, a1Count, KCount, KstarCount);
606 
607  // resonances
608  if (KCount >= 1)
609  channel = K;
610  if (KstarCount >= 1)
611  channel = Kstar;
612  if (a1Count >= 1)
613  channel = a1;
614  if (rhoCount >= 1)
615  channel = rho;
616  if (channel != undetermined && weight != 0.0)
617  TauDecayChannels->Fill(channel, weight);
618 
619  // final state products
620  if (piCount == 1 && pi0Count == 0)
621  channel = pi;
622  if (piCount == 1 && pi0Count == 1)
623  channel = pi1pi0;
624  if (piCount == 1 && pi0Count > 1)
625  channel = pinpi0;
626  if (piCount == 3 && pi0Count == 0)
627  channel = tripi;
628  if (piCount == 3 && pi0Count > 0)
629  channel = tripinpi0;
630  if (eCount == 1)
631  channel = electron;
632  if (muCount == 1)
633  channel = muon;
634  if (weight != 0.0)
635  TauDecayChannels->Fill(channel, weight);
636  return channel;
637 }
638 
639 void TauValidation::countParticles(const reco::GenParticle *p,
640  int &allCount,
641  int &eCount,
642  int &muCount,
643  int &pi0Count,
644  int &piCount,
645  int &rhoCount,
646  int &a1Count,
647  int &KCount,
648  int &KstarCount) {
649  for (unsigned int i = 0; i < p->numberOfDaughters(); i++) {
650  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(p->daughter(i));
651  int pid = dau->pdgId();
652  allCount++;
653  if (abs(pid) == 11)
654  eCount++;
655  else if (abs(pid) == 13)
656  muCount++;
657  else if (abs(pid) == 111)
658  pi0Count++;
659  else if (abs(pid) == 211)
660  piCount++;
661  else if (abs(pid) == 213)
662  rhoCount++;
663  else if (abs(pid) == 20213)
664  a1Count++;
665  else if (abs(pid) == 321)
666  KCount++;
667  else if (abs(pid) == 323)
668  KstarCount++;
669  countParticles(dau, allCount, eCount, muCount, pi0Count, piCount, rhoCount, a1Count, KCount, KstarCount);
670  }
671 }
672 
673 void TauValidation::spinEffectsWHpm(
674  const reco::GenParticle *tau, int mother, int decay, std::vector<const reco::GenParticle *> &part, double weight) {
675  // polarization only for 1-prong hadronic taus with no neutral pions
676  if (decay == TauDecay::MODE_PION || decay == TauDecay::MODE_MUON || decay == TauDecay::MODE_ELECTRON) {
677  TLorentzVector momP4 = motherP4(tau);
678  TLorentzVector pionP4 = leadingPionP4(tau);
679  pionP4.Boost(-1 * momP4.BoostVector());
680  double energy = pionP4.E() / (momP4.M() / 2);
681  if (decay == TauDecay::MODE_PION) {
682  if (abs(mother) == 24)
683  TauSpinEffectsW_X->Fill(energy, weight);
684  else if (abs(mother) == 37)
685  TauSpinEffectsHpm_X->Fill(energy, weight);
686  } else if (decay == TauDecay::MODE_MUON) {
687  if (abs(mother) == 24)
688  TauSpinEffectsW_muX->Fill(energy, weight);
689  else if (abs(mother) == 37)
690  TauSpinEffectsHpm_muX->Fill(energy, weight);
691  } else if (decay == TauDecay::MODE_ELECTRON) {
692  if (abs(mother) == 24)
693  TauSpinEffectsW_eX->Fill(energy, weight);
694  else if (abs(mother) == 37)
695  TauSpinEffectsHpm_eX->Fill(energy, weight);
696  }
697  } else if (decay == TauDecay::MODE_PIPI0) {
698  TLorentzVector rho(0, 0, 0, 0), pi(0, 0, 0, 0);
699  for (unsigned int i = 0; i < part.size(); i++) {
700  TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
701  if (abs(part.at(i)->pdgId()) == PdtPdgMini::pi_plus) {
702  pi += LV;
703  rho += LV;
704  } else if (abs(part.at(i)->pdgId()) == PdtPdgMini::pi0) {
705  rho += LV;
706  }
707  }
708  if (abs(mother) == 24)
709  TauSpinEffectsW_UpsilonRho->Fill(2 * pi.P() / rho.P() - 1, weight);
710  else if (abs(mother) == 37)
711  TauSpinEffectsHpm_UpsilonRho->Fill(2 * pi.P() / rho.P() - 1, weight);
712  } else if (decay == TauDecay::MODE_3PI || decay == TauDecay::MODE_PI2PI0) { // only for pi2pi0 for now
713  TLorentzVector a1(0, 0, 0, 0), pi_p(0, 0, 0, 0), pi_m(0, 0, 0, 0);
714  int nplus(0), nminus(0);
715  for (unsigned int i = 0; i < part.size(); i++) {
716  TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
717  if (part.at(i)->pdgId() == PdtPdgMini::pi_plus) {
718  pi_p += LV;
719  a1 += LV;
720  nplus++;
721  } else if (part.at(i)->pdgId() == PdtPdgMini::pi_minus) {
722  pi_m += LV;
723  a1 += LV;
724  nminus++;
725  }
726  }
727  double gamma = 0;
728  if (nplus + nminus == 3 && nplus == 1)
729  gamma = 2 * pi_p.P() / a1.P() - 1;
730  else if (nplus + nminus == 3 && nminus == 1)
731  gamma = 2 * pi_m.P() / a1.P() - 1;
732  else {
733  pi_p += pi_m;
734  gamma = 2 * pi_p.P() / a1.P() - 1;
735  }
736  if (abs(mother) == 24)
737  TauSpinEffectsW_UpsilonA1->Fill(gamma, weight);
738  else if (abs(mother) == 37)
739  TauSpinEffectsHpm_UpsilonA1->Fill(gamma, weight);
740  }
741 }
742 
743 void TauValidation::spinEffectsZH(const reco::GenParticle *boson, double weight) {
744  int ntau(0);
745  for (unsigned int i = 0; i < boson->numberOfDaughters(); i++) {
746  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(boson->daughter(i));
747  if (ntau == 1 && dau->pdgId() == 15)
748  return;
749  if (boson->pdgId() != 15 && abs(dau->pdgId()) == 15)
750  ntau++;
751  }
752  if (ntau != 2)
753  return;
754  if (abs(boson->pdgId()) == PdtPdgMini::Z0 || abs(boson->pdgId()) == PdtPdgMini::Higgs0) {
755  TLorentzVector tautau(0, 0, 0, 0);
756  TLorentzVector pipi(0, 0, 0, 0);
757  TLorentzVector taum(0, 0, 0, 0);
758  TLorentzVector taup(0, 0, 0, 0);
759  TLorentzVector rho_plus, rho_minus, pi_rhominus, pi0_rhominus, pi_rhoplus, pi0_rhoplus, pi_plus, pi_minus;
760  bool hasrho_minus(false), hasrho_plus(false), haspi_minus(false), haspi_plus(false);
761  int nSinglePionDecays(0), nSingleMuonDecays(0), nSingleElectronDecays(0);
762  double x1(0), x2(0);
763  TLorentzVector Zboson(boson->px(), boson->py(), boson->pz(), boson->energy());
764  for (unsigned int i = 0; i < boson->numberOfDaughters(); i++) {
765  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(boson->daughter(i));
766  int pid = dau->pdgId();
767  if (abs(findMother(dau)) != 15 && abs(pid) == 15) {
768  TauDecay_GenParticle TD;
769  unsigned int jak_id, TauBitMask;
770  if (TD.AnalyzeTau(dau, jak_id, TauBitMask, false, false)) {
771  std::vector<const reco::GenParticle *> part = TD.Get_TauDecayProducts();
772  if (jak_id == TauDecay::MODE_PION || jak_id == TauDecay::MODE_MUON || jak_id == TauDecay::MODE_ELECTRON) {
773  if (jak_id == TauDecay::MODE_PION)
774  nSinglePionDecays++;
775  if (jak_id == TauDecay::MODE_MUON)
776  nSingleMuonDecays++;
777  if (jak_id == TauDecay::MODE_ELECTRON)
778  nSingleElectronDecays++;
779  TLorentzVector LVtau(dau->px(), dau->py(), dau->pz(), dau->energy());
780  tautau += LVtau;
781  TLorentzVector LVpi = leadingPionP4(dau);
782  pipi += LVpi;
783  const HepPDT::ParticleData *pd = fPDGTable->particle(dau->pdgId());
784  int charge = (int)pd->charge();
785  LVtau.Boost(-1 * Zboson.BoostVector());
786  LVpi.Boost(-1 * Zboson.BoostVector());
787 
788  if (jak_id == TauDecay::MODE_MUON) {
789  if (abs(boson->pdgId()) == PdtPdgMini::Z0)
790  TauSpinEffectsZ_muX->Fill(LVpi.P() / LVtau.E(), weight);
791  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
792  TauSpinEffectsH_muX->Fill(LVpi.P() / LVtau.E(), weight);
793  }
794  if (jak_id == TauDecay::MODE_ELECTRON) {
795  if (abs(boson->pdgId()) == PdtPdgMini::Z0)
796  TauSpinEffectsZ_eX->Fill(LVpi.P() / LVtau.E(), weight);
797  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
798  TauSpinEffectsH_eX->Fill(LVpi.P() / LVtau.E(), weight);
799  }
800 
801  if (jak_id == TauDecay::MODE_PION) {
802  if (abs(boson->pdgId()) == PdtPdgMini::Z0) {
803  TauSpinEffectsZ_X->Fill(LVpi.P() / LVtau.E(), weight);
804  if (50.0 < Zboson.M() && Zboson.M() < 75.0)
805  TauSpinEffectsZ_X50to75->Fill(LVpi.P() / LVtau.E(), weight);
806  if (75.0 < Zboson.M() && Zboson.M() < 88.0)
807  TauSpinEffectsZ_X75to88->Fill(LVpi.P() / LVtau.E(), weight);
808  if (88.0 < Zboson.M() && Zboson.M() < 100.0)
809  TauSpinEffectsZ_X88to100->Fill(LVpi.P() / LVtau.E(), weight);
810  if (100.0 < Zboson.M() && Zboson.M() < 120.0)
811  TauSpinEffectsZ_X100to120->Fill(LVpi.P() / LVtau.E(), weight);
812  if (120.0 < Zboson.M())
813  TauSpinEffectsZ_X120UP->Fill(LVpi.P() / LVtau.E(), weight);
814  }
815  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
816  TauSpinEffectsH_X->Fill(LVpi.P() / LVtau.E(), weight);
817  }
818  if (charge < 0) {
819  x1 = LVpi.P() / LVtau.E();
820  taum = LVtau;
821  } else {
822  x2 = LVpi.P() / LVtau.E();
823  }
824  }
825  TLorentzVector LVtau(dau->px(), dau->py(), dau->pz(), dau->energy());
826  if (pid == 15)
827  taum = LVtau;
828  if (pid == -15)
829  taup = LVtau;
830  if (jak_id == TauDecay::MODE_PIPI0) {
831  for (unsigned int i = 0; i < part.size(); i++) {
832  int pid_d = part.at(i)->pdgId();
833  if (abs(pid_d) == 211 || abs(pid_d) == 111) {
834  TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
835  if (pid == 15) {
836  hasrho_minus = true;
837  if (pid_d == -211) {
838  pi_rhominus = LV;
839  }
840  if (abs(pid_d) == 111) {
841  pi0_rhominus = LV;
842  }
843  }
844  if (pid == -15) {
845  hasrho_plus = true;
846  if (pid_d == 211) {
847  pi_rhoplus = LV;
848  }
849  if (abs(pid_d) == 111) {
850  pi0_rhoplus = LV;
851  }
852  }
853  }
854  }
855  }
856  if (jak_id == TauDecay::MODE_PION) {
857  for (unsigned int i = 0; i < part.size(); i++) {
858  int pid_d = part.at(i)->pdgId();
859  if (abs(pid_d) == 211) {
860  TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
861  if (pid == 15) {
862  haspi_minus = true;
863  if (pid_d == -211) {
864  pi_minus = LV;
865  }
866  }
867  if (pid == -15) {
868  haspi_plus = true;
869  if (pid_d == 211) {
870  pi_plus = LV;
871  }
872  }
873  }
874  }
875  }
876  }
877  }
878  }
879  if (hasrho_minus && hasrho_plus) {
880  //compute rhorho
881  rho_minus = pi_rhominus;
882  rho_minus += pi0_rhominus;
883  rho_plus = pi_rhoplus;
884  rho_plus += pi0_rhoplus;
885  TLorentzVector rhorho = rho_minus;
886  rhorho += rho_plus;
887 
888  // boost to rhorho cm
889  TLorentzVector pi_rhoplusb = pi_rhoplus;
890  pi_rhoplusb.Boost(-1 * rhorho.BoostVector());
891  TLorentzVector pi0_rhoplusb = pi0_rhoplus;
892  pi0_rhoplusb.Boost(-1 * rhorho.BoostVector());
893  TLorentzVector pi_rhominusb = pi_rhominus;
894  pi_rhominusb.Boost(-1 * rhorho.BoostVector());
895  TLorentzVector pi0_rhominusb = pi0_rhominus;
896  pi0_rhominusb.Boost(-1 * rhorho.BoostVector());
897 
898  // compute n+/-
899  TVector3 n_plus = pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
900  TVector3 n_minus = pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
901 
902  // compute the acoplanarity
903  double Acoplanarity = acos(n_plus.Dot(n_minus) / (n_plus.Mag() * n_minus.Mag()));
904  if (pi_rhominusb.Vect().Dot(n_plus) > 0) {
905  Acoplanarity *= -1;
906  Acoplanarity += 2 * TMath::Pi();
907  }
908 
909  // now boost to tau frame
910  pi_rhoplus.Boost(-1 * taup.BoostVector());
911  pi0_rhoplus.Boost(-1 * taup.BoostVector());
912  pi_rhominus.Boost(-1 * taum.BoostVector());
913  pi0_rhominus.Boost(-1 * taum.BoostVector());
914 
915  // compute y1 and y2
916  double y1 = (pi_rhoplus.E() - pi0_rhoplus.E()) / (pi_rhoplus.E() + pi0_rhoplus.E());
917  double y2 = (pi_rhominus.E() - pi0_rhominus.E()) / (pi_rhominus.E() + pi0_rhominus.E());
918 
919  // fill histograms
920  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0 && y1 * y2 < 0)
921  TauSpinEffectsH_rhorhoAcoplanarityminus->Fill(Acoplanarity, weight);
922  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0 && y1 * y2 > 0)
923  TauSpinEffectsH_rhorhoAcoplanarityplus->Fill(Acoplanarity, weight);
924  }
925  if (haspi_minus && haspi_plus) {
926  TLorentzVector tauporig = taup;
927  TLorentzVector taumorig = taum;
928 
929  // now boost to Higgs frame
930  pi_plus.Boost(-1 * Zboson.BoostVector());
931  pi_minus.Boost(-1 * Zboson.BoostVector());
932 
933  taup.Boost(-1 * Zboson.BoostVector());
934  taum.Boost(-1 * Zboson.BoostVector());
935 
936  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0) {
937  TauSpinEffectsH_pipiAcollinearity->Fill(
938  acos(pi_plus.Vect().Dot(pi_minus.Vect()) / (pi_plus.P() * pi_minus.P())));
939  TauSpinEffectsH_pipiAcollinearityzoom->Fill(
940  acos(pi_plus.Vect().Dot(pi_minus.Vect()) / (pi_plus.P() * pi_minus.P())));
941  }
942 
943  double proj_m = taum.Vect().Dot(pi_minus.Vect()) / (taum.P() * taum.P());
944  double proj_p = taup.Vect().Dot(pi_plus.Vect()) / (taup.P() * taup.P());
945  TVector3 Tau_m = taum.Vect();
946  TVector3 Tau_p = taup.Vect();
947  Tau_m *= proj_m;
948  Tau_p *= proj_p;
949  TVector3 Pit_m = pi_minus.Vect() - Tau_m;
950  TVector3 Pit_p = pi_plus.Vect() - Tau_p;
951 
952  double Acoplanarity = acos(Pit_m.Dot(Pit_p) / (Pit_p.Mag() * Pit_m.Mag()));
953  TVector3 n = Pit_p.Cross(Pit_m);
954  if (n.Dot(Tau_m) / Tau_m.Mag() > 0) {
955  Acoplanarity *= -1;
956  Acoplanarity += 2 * TMath::Pi();
957  }
958  // fill histograms
959  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
960  TauSpinEffectsH_pipiAcoplanarity->Fill(Acoplanarity, weight);
961  taup = tauporig;
962  taum = taumorig;
963  }
964  if (nSinglePionDecays == 2 && tautau.M() != 0) {
965  for (int i = 0; i < zsbins; i++) {
966  double zslow = ((double)i) * (zsmax - zsmin) / ((double)zsbins) + zsmin;
967  double zsup = ((double)i + 1) * (zsmax - zsmin) / ((double)zsbins) + zsmin;
968  double aup = Zstoa(zsup), alow = Zstoa(zslow);
969  if (x2 - x1 > alow && x2 - x1 < aup) {
970  double zs = (zsup + zslow) / 2;
971  if (abs(boson->pdgId()) == PdtPdgMini::Z0)
972  TauSpinEffectsZ_Zs->Fill(zs, weight);
973  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
974  TauSpinEffectsH_Zs->Fill(zs, weight);
975  break;
976  }
977  }
978  if (abs(boson->pdgId()) == PdtPdgMini::Z0)
979  TauSpinEffectsZ_MVis->Fill(pipi.M() / tautau.M(), weight);
980  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
981  TauSpinEffectsH_MVis->Fill(pipi.M() / tautau.M(), weight);
982 
983  if (x1 != 0) {
984  const std::vector<const reco::GenParticle *> m = GetMothers(boson);
985  int q(0), qbar(0);
986  TLorentzVector Z(0, 0, 0, 0);
987  for (unsigned int i = 0; i < m.size(); i++) {
988  if (m.at(i)->pdgId() == PdtPdgMini::d || m.at(i)->pdgId() == PdtPdgMini::u) {
989  q++;
990  }
991  if (m.at(i)->pdgId() == PdtPdgMini::anti_d || m.at(i)->pdgId() == PdtPdgMini::anti_u) {
992  qbar++;
993  }
994  }
995  if (q == 1 && qbar == 1) { // assume q has largest E (valence vs see quarks)
996  if (taum.Vect().Dot(Zboson.Vect()) / (Zboson.P() * taum.P()) > 0) {
997  if (abs(boson->pdgId()) == PdtPdgMini::Z0)
998  TauSpinEffectsZ_Xf->Fill(x1, weight);
999  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
1000  TauSpinEffectsH_Xf->Fill(x1, weight);
1001  } else {
1002  if (abs(boson->pdgId()) == PdtPdgMini::Z0)
1003  TauSpinEffectsZ_Xb->Fill(x1, weight);
1004  if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
1005  TauSpinEffectsH_Xb->Fill(x1, weight);
1006  }
1007  }
1008  }
1009  }
1010  }
1011 }
1012 
1013 double TauValidation::Zstoa(double zs) {
1014  double a = 1 - sqrt(fabs(1.0 - 2 * fabs(zs)));
1015  if (zs < 0) {
1016  a *= -1.0;
1017  }
1018  return a;
1019 }
1020 
1021 double TauValidation::leadingPionMomentum(const reco::GenParticle *tau, double weight) {
1022  return leadingPionP4(tau).P();
1023 }
1024 
1025 TLorentzVector TauValidation::leadingPionP4(const reco::GenParticle *tau) {
1026  TLorentzVector p4(0, 0, 0, 0);
1027  for (unsigned int i = 0; i < tau->numberOfDaughters(); i++) {
1028  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(tau->daughter(i));
1029  int pid = dau->pdgId();
1030  if (abs(pid) == 15)
1031  return leadingPionP4(dau);
1032  if (!(abs(pid) == 211 || abs(pid) == 13 || abs(pid) == 11))
1033  continue;
1034  if (dau->p() > p4.P())
1035  p4 = TLorentzVector(dau->px(), dau->py(), dau->pz(), dau->energy());
1036  }
1037  return p4;
1038 }
1039 
1040 TLorentzVector TauValidation::motherP4(const reco::GenParticle *tau) {
1041  const reco::GenParticle *m = GetMother(tau);
1042  return TLorentzVector(m->px(), m->py(), m->pz(), m->energy());
1043 }
1044 
1045 double TauValidation::visibleTauEnergy(const reco::GenParticle *tau) {
1046  TLorentzVector p4(tau->px(), tau->py(), tau->pz(), tau->energy());
1047  for (unsigned int i = 0; i < tau->numberOfDaughters(); i++) {
1048  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(tau->daughter(i));
1049  int pid = dau->pdgId();
1050  if (abs(pid) == 15)
1051  return visibleTauEnergy(dau);
1052  if (abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
1053  p4 -= TLorentzVector(dau->px(), dau->py(), dau->pz(), dau->energy());
1054  }
1055  }
1056  return p4.E();
1057 }
1058 
1059 void TauValidation::photons(const reco::GenParticle *tau, double weight) {
1060  // Find First tau in chain
1061  std::vector<const reco::GenParticle *> TauList;
1062  findTauList(tau, TauList);
1063 
1064  // Get List of Gauge Boson to tau(s) FSR and Brem
1065  bool passedW = false;
1066  std::vector<const reco::GenParticle *> ListofFSR;
1067  ListofFSR.clear();
1068  std::vector<const reco::GenParticle *> ListofBrem;
1069  ListofBrem.clear();
1070  std::vector<const reco::GenParticle *> FSR_photos;
1071  FSR_photos.clear();
1072  double BosonScale(1);
1073  if (!TauList.empty()) {
1074  TauValidation::findFSRandBrem(TauList.at(0), passedW, ListofFSR, ListofBrem);
1075  TauValidation::FindPhotosFSR(TauList.at(0), FSR_photos, BosonScale);
1076 
1077  // Add the Tau Brem. information
1078  TauBremPhotonsN->Fill(ListofBrem.size(), weight);
1079  double photonPtSum = 0;
1080  for (unsigned int i = 0; i < ListofBrem.size(); i++) {
1081  photonPtSum += ListofBrem.at(i)->pt();
1082  TauBremPhotonsPt->Fill(ListofBrem.at(i)->pt(), weight);
1083  }
1084  TauBremPhotonsPtSum->Fill(photonPtSum, weight);
1085 
1086  // Now add the Gauge Boson FSR information
1087  if (BosonScale != 0) {
1088  TauFSRPhotonsN->Fill(ListofFSR.size(), weight);
1089  photonPtSum = 0;
1090  for (unsigned int i = 0; i < ListofFSR.size(); i++) {
1091  photonPtSum += ListofFSR.at(i)->pt();
1092  TauFSRPhotonsPt->Fill(ListofFSR.at(i)->pt(), weight);
1093  }
1094  double FSR_photosSum(0);
1095  for (unsigned int i = 0; i < FSR_photos.size(); i++) {
1096  FSR_photosSum += FSR_photos.at(i)->pt();
1097  TauFSRPhotonsPt->Fill(FSR_photos.at(i)->pt() / BosonScale, weight * BosonScale);
1098  }
1099  TauFSRPhotonsPtSum->Fill(photonPtSum + FSR_photosSum / BosonScale, weight);
1100  }
1101  }
1102 }
TauValidation::TauFSRPhotonsN
MonitorElement * TauFSRPhotonsN
Definition: TauValidation.h:81
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
TauValidation::findMother
int findMother(const reco::GenParticle *)
Definition: TauValidation.cc:489
TauDecay_GenParticle
Definition: TauDecay_GenParticle.h:23
c
const edm::EventSetup & c
Definition: SiStripLAProfileBooker.cc:66
TauValidation::TauSpinEffectsH_rhorhoAcoplanarityplus
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityplus
Definition: TauValidation.h:81
TauValidation::TauPhi
MonitorElement * TauPhi
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsW_X
MonitorElement * TauSpinEffectsW_X
Definition: TauValidation.h:81
reco::LeafCandidate::pz
double pz() const final
z coordinate of momentum vector
Definition: LeafCandidate.h:144
reco::CompositeRefCandidateT::mother
const Candidate * mother(size_type=0) const override
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
mps_fire.i
i
Definition: mps_fire.py:428
TauValidation::TauSpinEffectsHpm_UpsilonA1
MonitorElement * TauSpinEffectsHpm_UpsilonA1
Definition: TauValidation.h:81
TauValidation::MODEID
MonitorElement * MODEID
Definition: TauValidation.h:94
TauValidation::TauSpinEffectsZ_Xb
MonitorElement * TauSpinEffectsZ_Xb
Definition: TauValidation.h:81
TauValidation::nTaus
MonitorElement * nTaus
Definition: TauValidation.h:80
dqm::implementation::NavigatorBase::setCurrentFolder
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:32
TauValidation::DecayLength
MonitorElement * DecayLength
Definition: TauValidation.h:81
TauValidation::MODEInvMass
std::vector< std::vector< MonitorElement * > > MODEInvMass
Definition: TauValidation.h:95
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
LV
math::XYZTLorentzVectorD LV
Definition: HLTTauDQMPlotter.h:15
TauValidation::TauValidation
TauValidation(const edm::ParameterSet &)
Definition: TauValidation.cc:18
TauValidation::TauPt
MonitorElement * TauPt
Definition: TauValidation.h:81
TauValidation::~TauValidation
~TauValidation() override
Definition: TauValidation.cc:29
TauValidation::TauSpinEffectsH_MVis
MonitorElement * TauSpinEffectsH_MVis
Definition: TauValidation.h:81
TauValidation::TauMothers
MonitorElement * TauMothers
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsHpm_eX
MonitorElement * TauSpinEffectsHpm_eX
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsZ_X100to120
MonitorElement * TauSpinEffectsZ_X100to120
Definition: TauValidation.h:81
reco::CompositeRefCandidateT::numberOfDaughters
size_t numberOfDaughters() const override
number of daughters
TauValidation::FindPhotosFSR
void FindPhotosFSR(const reco::GenParticle *p, std::vector< const reco::GenParticle * > &ListofFSR, double &BosonScale)
Definition: TauValidation.cc:536
TauValidation::TauSpinEffectsW_muX
MonitorElement * TauSpinEffectsW_muX
Definition: TauValidation.h:81
TauValidation::TauBremPhotonsN
MonitorElement * TauBremPhotonsN
Definition: TauValidation.h:81
TauValidation::spinEffectsWHpm
void spinEffectsWHpm(const reco::GenParticle *, int, int, std::vector< const reco::GenParticle * > &part, double weight)
Definition: TauValidation.cc:673
TauValidation::TauSpinEffectsZ_Xf
MonitorElement * TauSpinEffectsZ_Xf
Definition: TauValidation.h:81
TauValidation::isLastTauinChain
bool isLastTauinChain(const reco::GenParticle *tau)
Definition: TauValidation.cc:491
reco::LeafCandidate::status
int status() const final
status word
Definition: LeafCandidate.h:180
TauValidation::fPDGTableToken
edm::ESGetToken< HepPDT::ParticleDataTable, edm::DefaultRecord > fPDGTableToken
Definition: TauValidation.h:78
TauValidation::TauBremPhotonsPtSum
MonitorElement * TauBremPhotonsPtSum
Definition: TauValidation.h:81
genParticleCandidates2GenParticles_cfi.genParticles
tuple genParticles
Definition: genParticleCandidates2GenParticles_cfi.py:3
TauValidation::TauSpinEffectsH_pipiAcoplanarity
MonitorElement * TauSpinEffectsH_pipiAcoplanarity
Definition: TauValidation.h:81
reco::CompositeRefCandidateT::numberOfMothers
size_t numberOfMothers() const override
number of mothers
TauValidation::TauSpinEffectsHpm_UpsilonRho
MonitorElement * TauSpinEffectsHpm_UpsilonRho
Definition: TauValidation.h:81
TauValidation::leadingPionMomentum
double leadingPionMomentum(const reco::GenParticle *, double weight)
Definition: TauValidation.cc:1021
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HLT_FULL_cff.beta
tuple beta
Definition: HLT_FULL_cff.py:8870
TauValidation::TauSpinEffectsZ_eX
MonitorElement * TauSpinEffectsZ_eX
Definition: TauValidation.h:81
TauValidation::TauProngs
MonitorElement * TauProngs
Definition: TauValidation.h:81
TauValidation::motherP4
TLorentzVector motherP4(const reco::GenParticle *)
Definition: TauValidation.cc:1040
reco::LeafCandidate::pdgId
int pdgId() const final
PDG identifier.
Definition: LeafCandidate.h:176
TauDecay::isTauFinalStateParticle
bool isTauFinalStateParticle(int pdgid)
Definition: TauDecay.cc:32
TauValidation::dqmBeginRun
void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override
Definition: TauValidation.cc:31
dqm::impl::MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:290
TauValidation::TauSpinEffectsH_pipiAcollinearityzoom
MonitorElement * TauSpinEffectsH_pipiAcollinearityzoom
Definition: TauValidation.h:81
label
char const * label
Definition: PFTauDecayModeTools.cc:11
visualization-live-secondInstance_cfg.m
tuple m
Definition: visualization-live-secondInstance_cfg.py:79
TauValidation::photons
void photons(const reco::GenParticle *, double weight)
Definition: TauValidation.cc:1059
reco::LeafCandidate::px
double px() const final
x coordinate of momentum vector
Definition: LeafCandidate.h:140
iEvent
int iEvent
Definition: GenABIO.cc:224
TauValidation::TauSpinEffectsW_UpsilonRho
MonitorElement * TauSpinEffectsW_UpsilonRho
Definition: TauValidation.h:81
TauValidation::GetMother
const reco::GenParticle * GetMother(const reco::GenParticle *tau)
Definition: TauValidation.cc:468
TauValidation::TauSpinEffectsZ_X88to100
MonitorElement * TauSpinEffectsZ_X88to100
Definition: TauValidation.h:81
reco::LeafCandidate::p
double p() const final
magnitude of momentum vector
Definition: LeafCandidate.h:123
TauDecay_GenParticle.h
TauValidation::TauSpinEffectsH_rhorhoAcoplanarityminus
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityminus
Definition: TauValidation.h:81
TauValidation::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: TauValidation.cc:371
TauValidation::TauSpinEffectsZ_X50to75
MonitorElement * TauSpinEffectsZ_X50to75
Definition: TauValidation.h:81
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
TauValidation::TauSpinEffectsH_Xf
MonitorElement * TauSpinEffectsH_Xf
Definition: TauValidation.h:81
TauValidation::TauFSRPhotonsPt
MonitorElement * TauFSRPhotonsPt
Definition: TauValidation.h:81
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TauValidation::TauSpinEffectsZ_muX
MonitorElement * TauSpinEffectsZ_muX
Definition: TauValidation.h:81
TauValidation::TauDecayChannels
MonitorElement * TauDecayChannels
Definition: TauValidation.h:81
reco::LeafCandidate::py
double py() const final
y coordinate of momentum vector
Definition: LeafCandidate.h:142
dqm::impl::MonitorElement::setBinLabel
virtual void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
Definition: MonitorElement.cc:762
TauValidation::TauSpinEffectsH_Zs
MonitorElement * TauSpinEffectsH_Zs
Definition: TauValidation.h:81
TauValidation::NMODEID
unsigned int NMODEID
Definition: TauValidation.h:93
TauValidation::TauFSRPhotonsPtSum
MonitorElement * TauFSRPhotonsPtSum
Definition: TauValidation.h:81
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:9
TauValidation::tauMother
int tauMother(const reco::GenParticle *, double weight)
Definition: TauValidation.cc:566
TauDecay_GenParticle::Get_TauDecayProducts
std::vector< const reco::GenParticle * > Get_TauDecayProducts()
Definition: TauDecay_GenParticle.h:31
TauValidation::bookHistograms
void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
Definition: TauValidation.cc:35
TauValidation::GetMothers
const std::vector< const reco::GenParticle * > GetMothers(const reco::GenParticle *boson)
Definition: TauValidation.cc:478
DQMHelper::book1dHisto
MonitorElement * book1dHisto(const std::string &name, const std::string &title, int n, double xmin, double xmax, const std::string &xaxis, const std::string &yaxis)
Definition: DQMHelper.cc:7
TauValidation::TauSpinEffectsZ_MVis
MonitorElement * TauSpinEffectsZ_MVis
Definition: TauValidation.h:81
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
TauValidation::TauSpinEffectsZ_X75to88
MonitorElement * TauSpinEffectsZ_X75to88
Definition: TauValidation.h:81
TauValidation::visibleTauEnergy
double visibleTauEnergy(const reco::GenParticle *)
Definition: TauValidation.cc:1045
TauDecay::DecayMode
static std::string DecayMode(unsigned int &MODE_ID)
Definition: TauDecay.cc:393
TauValidation::tauDecayChannel
int tauDecayChannel(const reco::GenParticle *tau, int jak_id, unsigned int TauBitMask, double weight)
Definition: TauValidation.cc:598
TauValidation::TauSpinEffectsH_pipiAcollinearity
MonitorElement * TauSpinEffectsH_pipiAcollinearity
Definition: TauValidation.h:81
part
part
Definition: HCALResponse.h:20
TauValidation::TauSpinEffectsZ_X
MonitorElement * TauSpinEffectsZ_X
Definition: TauValidation.h:81
TauValidation::fPDGTable
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
Definition: TauValidation.h:77
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
TauValidation::TauSpinEffectsH_eX
MonitorElement * TauSpinEffectsH_eX
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsW_eX
MonitorElement * TauSpinEffectsW_eX
Definition: TauValidation.h:81
TauValidation::genparticleCollectionToken_
edm::EDGetTokenT< reco::GenParticleCollection > genparticleCollectionToken_
Definition: TauValidation.h:100
TauValidation::Zstoa
double Zstoa(double zs)
Definition: TauValidation.cc:1013
TauValidation::findTauList
void findTauList(const reco::GenParticle *tau, std::vector< const reco::GenParticle * > &TauList)
Definition: TauValidation.cc:499
TauValidation::findFSRandBrem
void findFSRandBrem(const reco::GenParticle *p, bool doBrem, std::vector< const reco::GenParticle * > &ListofFSR, std::vector< const reco::GenParticle * > &ListofBrem)
Definition: TauValidation.cc:509
TauValidation::LifeTime
MonitorElement * LifeTime
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsH_muX
MonitorElement * TauSpinEffectsH_muX
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsW_UpsilonA1
MonitorElement * TauSpinEffectsW_UpsilonA1
Definition: TauValidation.h:81
a
double a
Definition: hdecay.h:119
TauValidation::countParticles
void countParticles(const reco::GenParticle *p, int &allCount, int &eCount, int &muCount, int &pi0Count, int &piCount, int &rhoCount, int &a1Count, int &KCount, int &KstarCount)
Definition: TauValidation.cc:639
TauValidation::TauEta
MonitorElement * TauEta
Definition: TauValidation.h:81
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
TauValidation::genparticleCollection_
edm::InputTag genparticleCollection_
Definition: TauValidation.h:74
TauValidation::TauSpinEffectsZ_Zs
MonitorElement * TauSpinEffectsZ_Zs
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsHpm_X
MonitorElement * TauSpinEffectsHpm_X
Definition: TauValidation.h:81
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:157
histoStyle.weight
int weight
Definition: histoStyle.py:51
TauValidation::TauSpinEffectsHpm_muX
MonitorElement * TauSpinEffectsHpm_muX
Definition: TauValidation.h:81
TauDecay_GenParticle::AnalyzeTau
bool AnalyzeTau(const reco::GenParticle *Tau, unsigned int &MODE_ID, unsigned int &TauBitMask, bool dores, bool dopi0)
Definition: TauDecay_GenParticle.cc:11
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
TauValidation::TauSpinEffectsZ_X120UP
MonitorElement * TauSpinEffectsZ_X120UP
Definition: TauValidation.h:81
TauValidation::TauSpinEffectsH_Xb
MonitorElement * TauSpinEffectsH_Xb
Definition: TauValidation.h:81
TauValidation::nPrimeTaus
MonitorElement * nPrimeTaus
Definition: TauValidation.h:80
TauValidation::TauBremPhotonsPt
MonitorElement * TauBremPhotonsPt
Definition: TauValidation.h:81
TauValidation::spinEffectsZH
void spinEffectsZH(const reco::GenParticle *boson, double weight)
Definition: TauValidation.cc:743
reco::CompositeRefCandidateT::daughter
const Candidate * daughter(size_type) const override
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
edm::Run
Definition: Run.h:45
TauDecay::nProng
unsigned int nProng(unsigned int &TauBitMask)
Definition: TauDecay.cc:331
TauValidation::TauSpinEffectsH_X
MonitorElement * TauSpinEffectsH_X
Definition: TauValidation.h:81
TauValidation::leadingPionP4
TLorentzVector leadingPionP4(const reco::GenParticle *)
Definition: TauValidation.cc:1025
reco::LeafCandidate::energy
double energy() const final
energy
Definition: LeafCandidate.h:125
Generated for CMSSW Reference Manual by doxygen 1.8.5