File indexing completed on 2024-04-06 12:32:11
0001
0002
0003
0004
0005
0006 #include "Validation/EventGenerator/interface/TauValidation.h"
0007 #include "CLHEP/Units/defs.h"
0008 #include "CLHEP/Units/PhysicalConstants.h"
0009 #include "DataFormats/Math/interface/LorentzVector.h"
0010 #include "Validation/EventGenerator/interface/TauDecay_GenParticle.h"
0011 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
0012 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0013 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0014 #include <iostream>
0015 #include "Validation/EventGenerator/interface/DQMHelper.h"
0016 using namespace edm;
0017
0018 TauValidation::TauValidation(const edm::ParameterSet &iPSet)
0019 :
0020 genparticleCollection_(iPSet.getParameter<edm::InputTag>("genparticleCollection")),
0021 NMODEID(TauDecay::NMODEID - 1),
0022 zsbins(20),
0023 zsmin(-0.5),
0024 zsmax(0.5) {
0025 genparticleCollectionToken_ = consumes<reco::GenParticleCollection>(genparticleCollection_);
0026 fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0027 }
0028
0029 TauValidation::~TauValidation() {}
0030
0031 void TauValidation::dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) {
0032 fPDGTable = c.getHandle(fPDGTableToken);
0033 }
0034
0035 void TauValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) {
0036
0037 DQMHelper dqm(&i);
0038 i.setCurrentFolder("Generator/Tau");
0039
0040 nTaus = dqm.book1dHisto("nTaus", "n analyzed Taus", 1, 0., 1., "bin", "Number of #tau's found");
0041 nPrimeTaus =
0042 dqm.book1dHisto("nPrimeTaus", "n analyzed prime Taus", 1, 0., 1., "bin", "Number of #tau's from Gauge Bosons");
0043
0044
0045 TauPt = dqm.book1dHisto("TauPt", "Tau pT", 100, 0, 100, "P_{T}^{#tau}", "Number of #tau's from Gauge Bosons");
0046 TauEta = dqm.book1dHisto("TauEta", "Tau eta", 100, -2.5, 2.5, "#eta^{#tau}", "Number of #tau's from Gauge Bosons");
0047 TauPhi = dqm.book1dHisto("TauPhi", "Tau phi", 100, -3.14, 3.14, "#phi^{#tau}", "Number of #tau's from Gauge Bosons");
0048 TauProngs = dqm.book1dHisto("TauProngs", "Tau n prongs", 7, 0, 7, "N_{prongs}", "Number of #tau's from Gauge Bosons");
0049 TauDecayChannels = dqm.book1dHisto(
0050 "TauDecayChannels", "Tau decay channels", 13, 0, 13, "Tau POG Decay Mode", "Number of #tau's from Gauge Bosons");
0051 TauDecayChannels->setBinLabel(1 + undetermined, "?");
0052 TauDecayChannels->setBinLabel(1 + electron, "e");
0053 TauDecayChannels->setBinLabel(1 + muon, "mu");
0054 TauDecayChannels->setBinLabel(1 + pi, "#pi^{#pm}");
0055 TauDecayChannels->setBinLabel(1 + rho, "#rho^{#pm}");
0056 TauDecayChannels->setBinLabel(1 + a1, "a_{1}^{#pm}");
0057 TauDecayChannels->setBinLabel(1 + pi1pi0, "#pi^{#pm}#pi^{0}");
0058 TauDecayChannels->setBinLabel(1 + pinpi0, "#pi^{#pm}n#pi^{0}");
0059 TauDecayChannels->setBinLabel(1 + tripi, "3#pi^{#pm}");
0060 TauDecayChannels->setBinLabel(1 + tripinpi0, "3#pi^{#pm}n#pi^{0}");
0061 TauDecayChannels->setBinLabel(1 + K, "K");
0062 TauDecayChannels->setBinLabel(1 + Kstar, "K^{*}");
0063 TauDecayChannels->setBinLabel(1 + stable, "Stable");
0064
0065 TauMothers = dqm.book1dHisto("TauMothers", "Tau mother particles", 10, 0, 10, "Mother of #tau", "Number of #tau's");
0066
0067 TauMothers->setBinLabel(1 + other, "?");
0068 TauMothers->setBinLabel(1 + B, "B Decays");
0069 TauMothers->setBinLabel(1 + D, "D Decays");
0070 TauMothers->setBinLabel(1 + gamma, "#gamma");
0071 TauMothers->setBinLabel(1 + Z, "Z");
0072 TauMothers->setBinLabel(1 + W, "W");
0073 TauMothers->setBinLabel(1 + HSM, "H_{SM}/h^{0}");
0074 TauMothers->setBinLabel(1 + H0, "H^{0}");
0075 TauMothers->setBinLabel(1 + A0, "A^{0}");
0076 TauMothers->setBinLabel(1 + Hpm, "H^{#pm}");
0077
0078 DecayLength = dqm.book1dHisto(
0079 "DecayLength", "#tau Decay Length", 100, -20, 20, "L_{#tau} (cm)", "Number of #tau's from Gauge Bosons");
0080 LifeTime = dqm.book1dHisto(
0081 "LifeTime", "#tau LifeTime ", 500, 0, 10000E-15, "#tau_{#tau} (s)", "Number of #tau's from Gauge Bosons");
0082
0083 TauSpinEffectsW_X = dqm.book1dHisto(
0084 "TauSpinEffectsWX", "X for pion", 50, 0, 1, "X", "Number of #tau#rightarrow#pi#nu from W^{#pm} Bosons");
0085 TauSpinEffectsHpm_X = dqm.book1dHisto(
0086 "TauSpinEffectsHpmX", "X for pion", 50, 0, 1, "X", "Number of #tau#rightarrow#pi#nu from H^{#pm} Bosons");
0087
0088 TauSpinEffectsW_eX = dqm.book1dHisto(
0089 "TauSpinEffectsWeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from W^{#pm} Bosons");
0090 TauSpinEffectsHpm_eX = dqm.book1dHisto(
0091 "TauSpinEffectsHpmeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from H^{#pm} Bosons");
0092
0093 TauSpinEffectsW_muX = dqm.book1dHisto(
0094 "TauSpinEffectsWmuX", "X for mu", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from W^{#pm} Bosons");
0095 TauSpinEffectsHpm_muX = dqm.book1dHisto(
0096 "TauSpinEffectsHpmmuX", "X for mue", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from H^{#pm} Bosons");
0097
0098 TauSpinEffectsW_UpsilonRho = dqm.book1dHisto("TauSpinEffectsWUpsilonRho",
0099 "#Upsilon for #rho",
0100 50,
0101 -1,
0102 1,
0103 "#Upsilon",
0104 "Number of #tau#rightarrow#rho#nu from Gauge Bosons");
0105 TauSpinEffectsHpm_UpsilonRho = dqm.book1dHisto("TauSpinEffectsHpmUpsilonRho",
0106 "#Upsilon for #rho",
0107 50,
0108 -1,
0109 1,
0110 "#Upsilon",
0111 "Number of #tau#rightarrow#rho#nu from Gauge Bosons");
0112
0113 TauSpinEffectsW_UpsilonA1 = dqm.book1dHisto("TauSpinEffectsWUpsilonA1",
0114 "#Upsilon for a1",
0115 50,
0116 -1,
0117 1,
0118 "#Upsilon",
0119 "Number of #tau#rightarrow#pi#pi#pi#nu from Gauge Bosons");
0120 TauSpinEffectsHpm_UpsilonA1 = dqm.book1dHisto("TauSpinEffectsHpmUpsilonA1",
0121 "#Upsilon for a1",
0122 50,
0123 -1,
0124 1,
0125 "#Upsilon",
0126 "Number of #tau#rightarrow#pi#pi#pi#nu from Gauge Bosons");
0127
0128 TauSpinEffectsH_pipiAcoplanarity =
0129 dqm.book1dHisto("TauSpinEffectsH_pipiAcoplanarity",
0130 "H Acoplanarity for #pi^{-}#pi^{+}",
0131 50,
0132 0,
0133 2 * TMath::Pi(),
0134 "Acoplanarity",
0135 "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
0136
0137 TauSpinEffectsH_pipiAcollinearity =
0138 dqm.book1dHisto("TauSpinEffectsH_pipiAcollinearity",
0139 "H Acollinearity for #pi^{-}#pi^{+}",
0140 50,
0141 0,
0142 TMath::Pi(),
0143 "Acollinearity",
0144 "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
0145 TauSpinEffectsH_pipiAcollinearityzoom =
0146 dqm.book1dHisto("TauSpinEffectsH_pipiAcollinearityzoom",
0147 "H Acollinearity for #pi^{-}#pi^{+}",
0148 50,
0149 3,
0150 TMath::Pi(),
0151 "Acollinearity",
0152 "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
0153
0154 TauSpinEffectsZ_MVis =
0155 dqm.book1dHisto("TauSpinEffectsZMVis",
0156 "Mass of pi+ pi-",
0157 25,
0158 0,
0159 1.1,
0160 "M_{#pi^{+}#pi^{-}} (GeV)",
0161 "Number of Z#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
0162 TauSpinEffectsH_MVis =
0163 dqm.book1dHisto("TauSpinEffectsHMVis",
0164 "Mass of pi+ pi-",
0165 25,
0166 0,
0167 1.1,
0168 "M_{#pi^{+}#pi^{-}} (GeV)",
0169 "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
0170
0171 TauSpinEffectsZ_Zs =
0172 dqm.book1dHisto("TauSpinEffectsZZs",
0173 "Z_{s}",
0174 zsbins,
0175 zsmin,
0176 zsmax,
0177 "Z_{s}",
0178 "Number of Z#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu Events");
0179 TauSpinEffectsH_Zs =
0180 dqm.book1dHisto("TauSpinEffectsHZs",
0181 "Z_{s}",
0182 zsbins,
0183 zsmin,
0184 zsmax,
0185 "Z_{s}",
0186 "Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu Events");
0187
0188 TauSpinEffectsZ_X = dqm.book1dHisto(
0189 "TauSpinEffectsZX", "X for pion of #tau^{-}", 25, 0, 1.0, "X", "Number of #tau#rightarrow#pi#nu from Z Bosons");
0190 TauSpinEffectsZ_X50to75 = dqm.book1dHisto("TauSpinEffectsZX50to75",
0191 "X for pion of #tau^{-} (50GeV-75GeV)",
0192 10,
0193 0,
0194 1.0,
0195 "X",
0196 "Number of #tau#rightarrow#pi#nu from Z(50GeV<M<75GeV) Bosons");
0197 TauSpinEffectsZ_X75to88 = dqm.book1dHisto("TauSpinEffectsZX75to88",
0198 "X for pion of #tau^{-} (75GeV-88GeV)",
0199 10,
0200 0,
0201 1.0,
0202 "X",
0203 "Number of #tau#rightarrow#pi#nu from Z(75GeV<M<88GeV) Bosons");
0204 TauSpinEffectsZ_X88to100 = dqm.book1dHisto("TauSpinEffectsZX88to100",
0205 "X for pion of #tau^{-} (88GeV-100GeV)",
0206 10,
0207 0,
0208 1.0,
0209 "X",
0210 "Number of #tau#rightarrow#pi#nu from Z(88GeV<M<100GeV) Bosons");
0211 TauSpinEffectsZ_X100to120 = dqm.book1dHisto("TauSpinEffectsZX100to120",
0212 "X for pion of #tau^{-} (100GeV-120GeV)",
0213 10,
0214 0,
0215 1.0,
0216 "X",
0217 "Number of #tau#rightarrow#pi#nu from Z(100GeV<M<120GeV) Bosons");
0218 TauSpinEffectsZ_X120UP = dqm.book1dHisto("TauSpinEffectsZX120UP",
0219 "X for pion of #tau^{-} (>120GeV)",
0220 10,
0221 0,
0222 1.0,
0223 "X",
0224 "Number of #tau#rightarrow#pi#nu from Z(120GeV<MGeV) Bosons");
0225
0226 TauSpinEffectsH_X = dqm.book1dHisto(
0227 "TauSpinEffectsH_X", "X for pion of #tau^{-}", 25, 0, 1.0, "X", "Number of #tau#rightarrow#pi#nu from H Bosons");
0228
0229 TauSpinEffectsZ_Xf = dqm.book1dHisto("TauSpinEffectsZXf",
0230 "X for pion of forward emitted #tau^{-}",
0231 25,
0232 0,
0233 1.0,
0234 "X_{f}",
0235 "Number of #tau#rightarrow#pi#nu from Z Bosons");
0236 TauSpinEffectsH_Xf = dqm.book1dHisto("TauSpinEffectsHXf",
0237 "X for pion of forward emitted #tau^{-}",
0238 25,
0239 0,
0240 1.0,
0241 "X_{f}",
0242 "Number of #tau#rightarrow#pi#nu from H Bosons");
0243
0244 TauSpinEffectsZ_Xb = dqm.book1dHisto("TauSpinEffectsZXb",
0245 "X for pion of backward emitted #tau^{-}",
0246 25,
0247 0,
0248 1.0,
0249 "X_{b}",
0250 "Number of #tau#rightarrow#pi#nu from Z Bosons");
0251 TauSpinEffectsH_Xb = dqm.book1dHisto("TauSpinEffectsHXb",
0252 "X for pion of backward emitted #tau^{-}",
0253 25,
0254 0,
0255 1.0,
0256 "X_{b}",
0257 "Number of #tau#rightarrow#pi#nu from H Bosons");
0258
0259 TauSpinEffectsZ_eX = dqm.book1dHisto(
0260 "TauSpinEffectsZeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from Gauge Bosons");
0261 TauSpinEffectsH_eX = dqm.book1dHisto(
0262 "TauSpinEffectsHeX", "X for e", 50, 0, 1, "X", "Number of #tau#rightarrowe#nu#nu from Gauge Bosons");
0263
0264 TauSpinEffectsZ_muX = dqm.book1dHisto(
0265 "TauSpinEffectsZmuX", "X for mu", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from Gauge Bosons");
0266 TauSpinEffectsH_muX = dqm.book1dHisto(
0267 "TauSpinEffectsHmuX", "X for mu", 50, 0, 1, "X", "Number of #tau#rightarrow#mu#nu#nu from Gauge Bosons");
0268
0269 TauSpinEffectsH_rhorhoAcoplanarityminus =
0270 dqm.book1dHisto("TauSpinEffectsH_rhorhoAcoplanarityminus",
0271 "#phi^{*-} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}<0)",
0272 32,
0273 0,
0274 2 * TMath::Pi(),
0275 "#phi^{*-} (Acoplanarity)",
0276 "Number of H#rightarrow#tau^{-}(#rightarrow#rho^{-}#nu)#tau^{+}(#rightarrow#rho^{+}#nu) Events");
0277 TauSpinEffectsH_rhorhoAcoplanarityplus =
0278 dqm.book1dHisto("TauSpinEffectsH_rhorhoAcoplanarityplus",
0279 "#phi^{*+} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}>0)",
0280 32,
0281 0,
0282 2 * TMath::Pi(),
0283 "#phi^{*+} (Acoplanarity)",
0284 "Number of H#rightarrow#tau^{-}(#rightarrow#rho^{-}#nu)#tau^{+}(#rightarrow#rho^{+}#nu) Events");
0285
0286 TauFSRPhotonsN = dqm.book1dHisto("TauFSRPhotonsN",
0287 "FSR Photons radiating from/with tau (Gauge Boson)",
0288 5,
0289 -0.5,
0290 4.5,
0291 "N^{FSR Photons radiating from/with #tau}",
0292 "Number of #tau's from Gauge Bosons");
0293 TauFSRPhotonsPt = dqm.book1dHisto("TauFSRPhotonsPt",
0294 "Pt of FSR Photons radiating from/with tau (Gauge Boson)",
0295 100,
0296 0,
0297 100,
0298 "P_{t}^{FSR Photons radiating from/with #tau [per #tau]} (GeV)",
0299 "Number of #tau's from Gauge Bosons");
0300 TauFSRPhotonsPtSum = dqm.book1dHisto("TauFSRPhotonsPtSum",
0301 "Pt of FSR Photons radiating from/with tau (Gauge Boson)",
0302 100,
0303 0,
0304 100,
0305 "P_{t}^{FSR Photons radiating from/with #tau [per #tau]} (GeV)",
0306 "Number of #tau's from Gauge Bosons");
0307
0308 TauBremPhotonsN = dqm.book1dHisto("TauBremPhotonsN",
0309 "Brem. Photons radiating in tau decay",
0310 5,
0311 -0.5,
0312 4.5,
0313 "N FSR Photons radiating from/with tau",
0314 "Number of #tau's from Gauge Bosons");
0315 TauBremPhotonsPt = dqm.book1dHisto("TauBremPhotonsPt",
0316 "Sum Brem Pt ",
0317 100,
0318 0,
0319 100,
0320 "P_{t}^{Brem. Photons radiating in tau decay} (GeV)",
0321 "Number of #tau's from Gauge Bosons");
0322 TauBremPhotonsPtSum = dqm.book1dHisto("TauBremPhotonsPtSum",
0323 "Sum of Brem Pt ",
0324 100,
0325 0,
0326 100,
0327 "Sum P_{t}^{Brem. Photons radiating in tau decay} (GeV)",
0328 "Number of #tau's from Gauge Bosons");
0329
0330 MODEID = dqm.book1dHisto("JAKID", "JAK ID", NMODEID + 1, -0.5, NMODEID + 0.5);
0331 for (unsigned int j = 0; j < NMODEID + 1; j++) {
0332 MODEInvMass.push_back(std::vector<MonitorElement *>());
0333 std::string tmp = "JAKID";
0334 tmp += std::to_string(j);
0335 MODEInvMass.at(j).push_back(dqm.book1dHisto("M" + tmp,
0336 "M_{" + TauDecay::DecayMode(j) + "} (GeV)",
0337 80,
0338 0,
0339 2.0,
0340 "M_{" + TauDecay::DecayMode(j) + "} (GeV)",
0341 "Number of #tau's from Gauge Bosons"));
0342 MODEID->setBinLabel(1 + j, TauDecay::DecayMode(j));
0343 if (j == TauDecay::MODE_3PI || j == TauDecay::MODE_PI2PI0 || j == TauDecay::MODE_KPIK ||
0344 j == TauDecay::MODE_KPIPI) {
0345 MODEInvMass.at(j).push_back(dqm.book1dHisto("M13" + tmp,
0346 "M_{13," + TauDecay::DecayMode(j) + "} (GeV)",
0347 80,
0348 0,
0349 2.0,
0350 "M_{13," + TauDecay::DecayMode(j) + "} (GeV)",
0351 "Number of #tau's from Gauge Bosons"));
0352 MODEInvMass.at(j).push_back(dqm.book1dHisto("M23" + tmp,
0353 "M_{23," + TauDecay::DecayMode(j) + "} (GeV)",
0354 80,
0355 0,
0356 2.0,
0357 "M_{23," + TauDecay::DecayMode(j) + "} (GeV)",
0358 "Number of #tau's from Gauge Bosons"));
0359 MODEInvMass.at(j).push_back(dqm.book1dHisto("M12" + tmp,
0360 "M_{12," + TauDecay::DecayMode(j) + "} (GeV)",
0361 80,
0362 0,
0363 2.0,
0364 "M_{12," + TauDecay::DecayMode(j) + "} (GeV)",
0365 "Number of #tau's from Gauge Bosons"));
0366 }
0367 }
0368 return;
0369 }
0370
0371 void TauValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0372
0373 edm::Handle<reco::GenParticleCollection> genParticles;
0374 iEvent.getByToken(genparticleCollectionToken_, genParticles);
0375
0376 double weight = 1.0;
0377
0378
0379 for (reco::GenParticleCollection::const_iterator iter = genParticles->begin(); iter != genParticles->end(); ++iter) {
0380 if (abs(iter->pdgId()) == PdtPdgMini::Z0 || abs(iter->pdgId()) == PdtPdgMini::Higgs0) {
0381 spinEffectsZH(&(*iter), weight);
0382 }
0383 if (abs(iter->pdgId()) == 15) {
0384 if (isLastTauinChain(&(*iter))) {
0385 nTaus->Fill(0.5, weight);
0386 int mother = tauMother(&(*iter), weight);
0387 if (mother > -1) {
0388 nPrimeTaus->Fill(0.5, weight);
0389 TauPt->Fill(iter->pt(), weight);
0390 TauEta->Fill(iter->eta(), weight);
0391 TauPhi->Fill(iter->phi(), weight);
0392 photons(&(*iter), weight);
0393
0394
0395 TauDecay_GenParticle TD;
0396 unsigned int jak_id, TauBitMask;
0397 if (TD.AnalyzeTau(&(*iter), jak_id, TauBitMask, false, false)) {
0398 MODEID->Fill(jak_id, weight);
0399 TauProngs->Fill(TD.nProng(TauBitMask), weight);
0400 tauDecayChannel(&(*iter), jak_id, TauBitMask, weight);
0401 if (jak_id <= NMODEID) {
0402 int tcharge = iter->pdgId() / abs(iter->pdgId());
0403 std::vector<const reco::GenParticle *> part = TD.Get_TauDecayProducts();
0404 spinEffectsWHpm(&(*iter), mother, jak_id, part, weight);
0405 TLorentzVector LVQ(0, 0, 0, 0);
0406 TLorentzVector LVS12(0, 0, 0, 0);
0407 TLorentzVector LVS13(0, 0, 0, 0);
0408 TLorentzVector LVS23(0, 0, 0, 0);
0409 bool haspart1 = false;
0410 TVector3 PV, SV;
0411 bool hasDL(false);
0412 for (unsigned int i = 0; i < part.size(); i++) {
0413 if (abs(part.at(i)->pdgId()) != PdtPdgMini::nu_tau && TD.isTauFinalStateParticle(part.at(i)->pdgId()) &&
0414 !hasDL) {
0415 hasDL = true;
0416 TLorentzVector tlv(iter->px(), iter->py(), iter->pz(), iter->energy());
0417 PV = TVector3(iter->vx(), iter->vy(), iter->vz());
0418 SV = TVector3(part.at(i)->vx(), part.at(i)->vy(), part.at(i)->vz());
0419 TVector3 DL = SV - PV;
0420 DecayLength->Fill(DL.Dot(tlv.Vect()) / tlv.P(), weight);
0421 double c(2.99792458E8), Ltau(DL.Mag() / 100) , beta(iter->p() / iter->mass());
0422 LifeTime->Fill(Ltau / (c * beta), weight);
0423 }
0424
0425 if (TD.isTauFinalStateParticle(part.at(i)->pdgId()) && abs(part.at(i)->pdgId()) != PdtPdgMini::nu_e &&
0426 abs(part.at(i)->pdgId()) != PdtPdgMini::nu_mu && abs(part.at(i)->pdgId()) != PdtPdgMini::nu_tau) {
0427 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0428 LVQ += LV;
0429 if (jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0 ||
0430 jak_id == TauDecay::MODE_KPIK || jak_id == TauDecay::MODE_KPIPI) {
0431 if ((tcharge == part.at(i)->pdgId() / abs(part.at(i)->pdgId()) && TD.nProng(TauBitMask) == 3) ||
0432 ((jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0) &&
0433 TD.nProng(TauBitMask) == 1 && abs(part.at(i)->pdgId()) == PdtPdgMini::pi_plus)) {
0434 LVS13 += LV;
0435 LVS23 += LV;
0436 } else {
0437 LVS12 += LV;
0438 if (!haspart1 && ((jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0) ||
0439 ((jak_id != TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0) &&
0440 abs(part.at(i)->pdgId()) == PdtPdgMini::K_plus))) {
0441 LVS13 += LV;
0442 haspart1 = true;
0443 } else {
0444 LVS23 += LV;
0445 }
0446 }
0447 }
0448 }
0449 }
0450 part.clear();
0451 MODEInvMass.at(jak_id).at(0)->Fill(LVQ.M(), weight);
0452 if (jak_id == TauDecay::MODE_3PI || jak_id == TauDecay::MODE_PI2PI0 || jak_id == TauDecay::MODE_KPIK ||
0453 jak_id == TauDecay::MODE_KPIPI) {
0454 MODEInvMass.at(jak_id).at(1)->Fill(LVS13.M(), weight);
0455 MODEInvMass.at(jak_id).at(2)->Fill(LVS23.M(), weight);
0456 MODEInvMass.at(jak_id).at(3)->Fill(LVS12.M(), weight);
0457 }
0458 }
0459 } else {
0460 MODEID->Fill(jak_id, weight);
0461 }
0462 }
0463 }
0464 }
0465 }
0466 }
0467
0468 const reco::GenParticle *TauValidation::GetMother(const reco::GenParticle *tau) {
0469 for (unsigned int i = 0; i < tau->numberOfMothers(); i++) {
0470 const reco::GenParticle *mother = static_cast<const reco::GenParticle *>(tau->mother(i));
0471 if (mother->pdgId() == tau->pdgId())
0472 return GetMother(mother);
0473 return mother;
0474 }
0475 return tau;
0476 }
0477
0478 const std::vector<const reco::GenParticle *> TauValidation::GetMothers(const reco::GenParticle *boson) {
0479 std::vector<const reco::GenParticle *> mothers;
0480 for (unsigned int i = 0; i < boson->numberOfMothers(); i++) {
0481 const reco::GenParticle *mother = static_cast<const reco::GenParticle *>(boson->mother(i));
0482 if (mother->pdgId() == boson->pdgId())
0483 return GetMothers(mother);
0484 mothers.push_back(mother);
0485 }
0486 return mothers;
0487 }
0488
0489 int TauValidation::findMother(const reco::GenParticle *tau) { return TauValidation::GetMother(tau)->pdgId(); }
0490
0491 bool TauValidation::isLastTauinChain(const reco::GenParticle *tau) {
0492 for (unsigned int i = 0; i < tau->numberOfDaughters(); i++) {
0493 if (tau->daughter(i)->pdgId() == tau->pdgId())
0494 return false;
0495 }
0496 return true;
0497 }
0498
0499 void TauValidation::findTauList(const reco::GenParticle *tau, std::vector<const reco::GenParticle *> &TauList) {
0500 TauList.insert(TauList.begin(), tau);
0501 for (unsigned int i = 0; i < tau->numberOfMothers(); i++) {
0502 const reco::GenParticle *mother = static_cast<const reco::GenParticle *>(tau->mother(i));
0503 if (mother->pdgId() == tau->pdgId()) {
0504 findTauList(mother, TauList);
0505 }
0506 }
0507 }
0508
0509 void TauValidation::findFSRandBrem(const reco::GenParticle *p,
0510 bool doBrem,
0511 std::vector<const reco::GenParticle *> &ListofFSR,
0512 std::vector<const reco::GenParticle *> &ListofBrem) {
0513
0514 if (abs(p->pdgId()) == 15) {
0515 if (isLastTauinChain(p)) {
0516 doBrem = true;
0517 } else {
0518 doBrem = false;
0519 }
0520 }
0521 int photo_ID = 22;
0522 for (unsigned int i = 0; i < p->numberOfDaughters(); i++) {
0523 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(p->daughter(i));
0524 if (abs((dau)->pdgId()) == abs(photo_ID) && !doBrem) {
0525 ListofFSR.push_back(dau);
0526 }
0527 if (abs((dau)->pdgId()) == abs(photo_ID) && doBrem) {
0528 ListofBrem.push_back(dau);
0529 }
0530 if (abs((dau)->pdgId()) != 111 && abs((dau)->pdgId()) != 221) {
0531 findFSRandBrem(dau, doBrem, ListofFSR, ListofBrem);
0532 }
0533 }
0534 }
0535
0536 void TauValidation::FindPhotosFSR(const reco::GenParticle *p,
0537 std::vector<const reco::GenParticle *> &ListofFSR,
0538 double &BosonScale) {
0539 BosonScale = 0.0;
0540 const reco::GenParticle *m = GetMother(p);
0541 int mother_pid = m->pdgId();
0542 if (m->pdgId() != p->pdgId()) {
0543 for (unsigned int i = 0; i < m->numberOfDaughters(); i++) {
0544 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(m->daughter(i));
0545 if (abs(dau->pdgId()) == 22) {
0546 ListofFSR.push_back(dau);
0547 }
0548 }
0549 }
0550 if (abs(mother_pid) == 24)
0551 BosonScale = 1.0;
0552 if (abs(mother_pid) == 23)
0553 BosonScale = 2.0;
0554 if (abs(mother_pid) == 22)
0555 BosonScale = 2.0;
0556 if (abs(mother_pid) == 25)
0557 BosonScale = 2.0;
0558 if (abs(mother_pid) == 35)
0559 BosonScale = 2.0;
0560 if (abs(mother_pid) == 36)
0561 BosonScale = 2.0;
0562 if (abs(mother_pid) == 37)
0563 BosonScale = 1.0;
0564 }
0565
0566 int TauValidation::tauMother(const reco::GenParticle *tau, double weight) {
0567 if (abs(tau->pdgId()) != 15)
0568 return -3;
0569 int mother_pid = findMother(tau);
0570 if (mother_pid == -2)
0571 return -2;
0572 int label = other;
0573 if (abs(mother_pid) == 24)
0574 label = W;
0575 if (abs(mother_pid) == 23)
0576 label = Z;
0577 if (abs(mother_pid) == 22)
0578 label = gamma;
0579 if (abs(mother_pid) == 25)
0580 label = HSM;
0581 if (abs(mother_pid) == 35)
0582 label = H0;
0583 if (abs(mother_pid) == 36)
0584 label = A0;
0585 if (abs(mother_pid) == 37)
0586 label = Hpm;
0587 int mother_shortpid = (abs(mother_pid) % 10000);
0588 if (mother_shortpid > 500 && mother_shortpid < 600)
0589 label = B;
0590 if (mother_shortpid > 400 && mother_shortpid < 500)
0591 label = D;
0592 TauMothers->Fill(label, weight);
0593 if (label == B || label == D || label == other)
0594 return -1;
0595 return mother_pid;
0596 }
0597
0598 int TauValidation::tauDecayChannel(const reco::GenParticle *tau, int jak_id, unsigned int TauBitMask, double weight) {
0599 int channel = undetermined;
0600 if (tau->status() == 1)
0601 channel = stable;
0602 int allCount = 0, eCount = 0, muCount = 0, pi0Count = 0, piCount = 0, rhoCount = 0, a1Count = 0, KCount = 0,
0603 KstarCount = 0;
0604
0605 countParticles(tau, allCount, eCount, muCount, pi0Count, piCount, rhoCount, a1Count, KCount, KstarCount);
0606
0607
0608 if (KCount >= 1)
0609 channel = K;
0610 if (KstarCount >= 1)
0611 channel = Kstar;
0612 if (a1Count >= 1)
0613 channel = a1;
0614 if (rhoCount >= 1)
0615 channel = rho;
0616 if (channel != undetermined && weight != 0.0)
0617 TauDecayChannels->Fill(channel, weight);
0618
0619
0620 if (piCount == 1 && pi0Count == 0)
0621 channel = pi;
0622 if (piCount == 1 && pi0Count == 1)
0623 channel = pi1pi0;
0624 if (piCount == 1 && pi0Count > 1)
0625 channel = pinpi0;
0626 if (piCount == 3 && pi0Count == 0)
0627 channel = tripi;
0628 if (piCount == 3 && pi0Count > 0)
0629 channel = tripinpi0;
0630 if (eCount == 1)
0631 channel = electron;
0632 if (muCount == 1)
0633 channel = muon;
0634 if (weight != 0.0)
0635 TauDecayChannels->Fill(channel, weight);
0636 return channel;
0637 }
0638
0639 void TauValidation::countParticles(const reco::GenParticle *p,
0640 int &allCount,
0641 int &eCount,
0642 int &muCount,
0643 int &pi0Count,
0644 int &piCount,
0645 int &rhoCount,
0646 int &a1Count,
0647 int &KCount,
0648 int &KstarCount) {
0649 for (unsigned int i = 0; i < p->numberOfDaughters(); i++) {
0650 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(p->daughter(i));
0651 int pid = dau->pdgId();
0652 allCount++;
0653 if (abs(pid) == 11)
0654 eCount++;
0655 else if (abs(pid) == 13)
0656 muCount++;
0657 else if (abs(pid) == 111)
0658 pi0Count++;
0659 else if (abs(pid) == 211)
0660 piCount++;
0661 else if (abs(pid) == 213)
0662 rhoCount++;
0663 else if (abs(pid) == 20213)
0664 a1Count++;
0665 else if (abs(pid) == 321)
0666 KCount++;
0667 else if (abs(pid) == 323)
0668 KstarCount++;
0669 countParticles(dau, allCount, eCount, muCount, pi0Count, piCount, rhoCount, a1Count, KCount, KstarCount);
0670 }
0671 }
0672
0673 void TauValidation::spinEffectsWHpm(
0674 const reco::GenParticle *tau, int mother, int decay, std::vector<const reco::GenParticle *> &part, double weight) {
0675
0676 if (decay == TauDecay::MODE_PION || decay == TauDecay::MODE_MUON || decay == TauDecay::MODE_ELECTRON) {
0677 TLorentzVector momP4 = motherP4(tau);
0678 TLorentzVector pionP4 = leadingPionP4(tau);
0679 pionP4.Boost(-1 * momP4.BoostVector());
0680 double energy = pionP4.E() / (momP4.M() / 2);
0681 if (decay == TauDecay::MODE_PION) {
0682 if (abs(mother) == 24)
0683 TauSpinEffectsW_X->Fill(energy, weight);
0684 else if (abs(mother) == 37)
0685 TauSpinEffectsHpm_X->Fill(energy, weight);
0686 } else if (decay == TauDecay::MODE_MUON) {
0687 if (abs(mother) == 24)
0688 TauSpinEffectsW_muX->Fill(energy, weight);
0689 else if (abs(mother) == 37)
0690 TauSpinEffectsHpm_muX->Fill(energy, weight);
0691 } else if (decay == TauDecay::MODE_ELECTRON) {
0692 if (abs(mother) == 24)
0693 TauSpinEffectsW_eX->Fill(energy, weight);
0694 else if (abs(mother) == 37)
0695 TauSpinEffectsHpm_eX->Fill(energy, weight);
0696 }
0697 } else if (decay == TauDecay::MODE_PIPI0) {
0698 TLorentzVector rho(0, 0, 0, 0), pi(0, 0, 0, 0);
0699 for (unsigned int i = 0; i < part.size(); i++) {
0700 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0701 if (abs(part.at(i)->pdgId()) == PdtPdgMini::pi_plus) {
0702 pi += LV;
0703 rho += LV;
0704 } else if (abs(part.at(i)->pdgId()) == PdtPdgMini::pi0) {
0705 rho += LV;
0706 }
0707 }
0708 if (abs(mother) == 24)
0709 TauSpinEffectsW_UpsilonRho->Fill(2 * pi.P() / rho.P() - 1, weight);
0710 else if (abs(mother) == 37)
0711 TauSpinEffectsHpm_UpsilonRho->Fill(2 * pi.P() / rho.P() - 1, weight);
0712 } else if (decay == TauDecay::MODE_3PI || decay == TauDecay::MODE_PI2PI0) {
0713 TLorentzVector a1(0, 0, 0, 0), pi_p(0, 0, 0, 0), pi_m(0, 0, 0, 0);
0714 int nplus(0), nminus(0);
0715 for (unsigned int i = 0; i < part.size(); i++) {
0716 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0717 if (part.at(i)->pdgId() == PdtPdgMini::pi_plus) {
0718 pi_p += LV;
0719 a1 += LV;
0720 nplus++;
0721 } else if (part.at(i)->pdgId() == PdtPdgMini::pi_minus) {
0722 pi_m += LV;
0723 a1 += LV;
0724 nminus++;
0725 }
0726 }
0727 double gamma = 0;
0728 if (nplus + nminus == 3 && nplus == 1)
0729 gamma = 2 * pi_p.P() / a1.P() - 1;
0730 else if (nplus + nminus == 3 && nminus == 1)
0731 gamma = 2 * pi_m.P() / a1.P() - 1;
0732 else {
0733 pi_p += pi_m;
0734 gamma = 2 * pi_p.P() / a1.P() - 1;
0735 }
0736 if (abs(mother) == 24)
0737 TauSpinEffectsW_UpsilonA1->Fill(gamma, weight);
0738 else if (abs(mother) == 37)
0739 TauSpinEffectsHpm_UpsilonA1->Fill(gamma, weight);
0740 }
0741 }
0742
0743 void TauValidation::spinEffectsZH(const reco::GenParticle *boson, double weight) {
0744 int ntau(0);
0745 for (unsigned int i = 0; i < boson->numberOfDaughters(); i++) {
0746 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(boson->daughter(i));
0747 if (ntau == 1 && dau->pdgId() == 15)
0748 return;
0749 if (boson->pdgId() != 15 && abs(dau->pdgId()) == 15)
0750 ntau++;
0751 }
0752 if (ntau != 2)
0753 return;
0754 if (abs(boson->pdgId()) == PdtPdgMini::Z0 || abs(boson->pdgId()) == PdtPdgMini::Higgs0) {
0755 TLorentzVector tautau(0, 0, 0, 0);
0756 TLorentzVector pipi(0, 0, 0, 0);
0757 TLorentzVector taum(0, 0, 0, 0);
0758 TLorentzVector taup(0, 0, 0, 0);
0759 TLorentzVector rho_plus, rho_minus, pi_rhominus, pi0_rhominus, pi_rhoplus, pi0_rhoplus, pi_plus, pi_minus;
0760 bool hasrho_minus(false), hasrho_plus(false), haspi_minus(false), haspi_plus(false);
0761 int nSinglePionDecays(0);
0762 double x1(0), x2(0);
0763 TLorentzVector Zboson(boson->px(), boson->py(), boson->pz(), boson->energy());
0764 for (unsigned int i = 0; i < boson->numberOfDaughters(); i++) {
0765 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(boson->daughter(i));
0766 int pid = dau->pdgId();
0767 if (abs(findMother(dau)) != 15 && abs(pid) == 15) {
0768 TauDecay_GenParticle TD;
0769 unsigned int jak_id, TauBitMask;
0770 if (TD.AnalyzeTau(dau, jak_id, TauBitMask, false, false)) {
0771 std::vector<const reco::GenParticle *> part = TD.Get_TauDecayProducts();
0772 if (jak_id == TauDecay::MODE_PION || jak_id == TauDecay::MODE_MUON || jak_id == TauDecay::MODE_ELECTRON) {
0773 if (jak_id == TauDecay::MODE_PION)
0774 nSinglePionDecays++;
0775 TLorentzVector LVtau(dau->px(), dau->py(), dau->pz(), dau->energy());
0776 tautau += LVtau;
0777 TLorentzVector LVpi = leadingPionP4(dau);
0778 pipi += LVpi;
0779 const HepPDT::ParticleData *pd = fPDGTable->particle(dau->pdgId());
0780 int charge = (int)pd->charge();
0781 LVtau.Boost(-1 * Zboson.BoostVector());
0782 LVpi.Boost(-1 * Zboson.BoostVector());
0783
0784 if (jak_id == TauDecay::MODE_MUON) {
0785 if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0786 TauSpinEffectsZ_muX->Fill(LVpi.P() / LVtau.E(), weight);
0787 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0788 TauSpinEffectsH_muX->Fill(LVpi.P() / LVtau.E(), weight);
0789 }
0790 if (jak_id == TauDecay::MODE_ELECTRON) {
0791 if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0792 TauSpinEffectsZ_eX->Fill(LVpi.P() / LVtau.E(), weight);
0793 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0794 TauSpinEffectsH_eX->Fill(LVpi.P() / LVtau.E(), weight);
0795 }
0796
0797 if (jak_id == TauDecay::MODE_PION) {
0798 if (abs(boson->pdgId()) == PdtPdgMini::Z0) {
0799 TauSpinEffectsZ_X->Fill(LVpi.P() / LVtau.E(), weight);
0800 if (50.0 < Zboson.M() && Zboson.M() < 75.0)
0801 TauSpinEffectsZ_X50to75->Fill(LVpi.P() / LVtau.E(), weight);
0802 if (75.0 < Zboson.M() && Zboson.M() < 88.0)
0803 TauSpinEffectsZ_X75to88->Fill(LVpi.P() / LVtau.E(), weight);
0804 if (88.0 < Zboson.M() && Zboson.M() < 100.0)
0805 TauSpinEffectsZ_X88to100->Fill(LVpi.P() / LVtau.E(), weight);
0806 if (100.0 < Zboson.M() && Zboson.M() < 120.0)
0807 TauSpinEffectsZ_X100to120->Fill(LVpi.P() / LVtau.E(), weight);
0808 if (120.0 < Zboson.M())
0809 TauSpinEffectsZ_X120UP->Fill(LVpi.P() / LVtau.E(), weight);
0810 }
0811 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0812 TauSpinEffectsH_X->Fill(LVpi.P() / LVtau.E(), weight);
0813 }
0814 if (charge < 0) {
0815 x1 = LVpi.P() / LVtau.E();
0816 taum = LVtau;
0817 } else {
0818 x2 = LVpi.P() / LVtau.E();
0819 }
0820 }
0821 TLorentzVector LVtau(dau->px(), dau->py(), dau->pz(), dau->energy());
0822 if (pid == 15)
0823 taum = LVtau;
0824 if (pid == -15)
0825 taup = LVtau;
0826 if (jak_id == TauDecay::MODE_PIPI0) {
0827 for (unsigned int i = 0; i < part.size(); i++) {
0828 int pid_d = part.at(i)->pdgId();
0829 if (abs(pid_d) == 211 || abs(pid_d) == 111) {
0830 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0831 if (pid == 15) {
0832 hasrho_minus = true;
0833 if (pid_d == -211) {
0834 pi_rhominus = LV;
0835 }
0836 if (abs(pid_d) == 111) {
0837 pi0_rhominus = LV;
0838 }
0839 }
0840 if (pid == -15) {
0841 hasrho_plus = true;
0842 if (pid_d == 211) {
0843 pi_rhoplus = LV;
0844 }
0845 if (abs(pid_d) == 111) {
0846 pi0_rhoplus = LV;
0847 }
0848 }
0849 }
0850 }
0851 }
0852 if (jak_id == TauDecay::MODE_PION) {
0853 for (unsigned int i = 0; i < part.size(); i++) {
0854 int pid_d = part.at(i)->pdgId();
0855 if (abs(pid_d) == 211) {
0856 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0857 if (pid == 15) {
0858 haspi_minus = true;
0859 if (pid_d == -211) {
0860 pi_minus = LV;
0861 }
0862 }
0863 if (pid == -15) {
0864 haspi_plus = true;
0865 if (pid_d == 211) {
0866 pi_plus = LV;
0867 }
0868 }
0869 }
0870 }
0871 }
0872 }
0873 }
0874 }
0875 if (hasrho_minus && hasrho_plus) {
0876
0877 rho_minus = pi_rhominus;
0878 rho_minus += pi0_rhominus;
0879 rho_plus = pi_rhoplus;
0880 rho_plus += pi0_rhoplus;
0881 TLorentzVector rhorho = rho_minus;
0882 rhorho += rho_plus;
0883
0884
0885 TLorentzVector pi_rhoplusb = pi_rhoplus;
0886 pi_rhoplusb.Boost(-1 * rhorho.BoostVector());
0887 TLorentzVector pi0_rhoplusb = pi0_rhoplus;
0888 pi0_rhoplusb.Boost(-1 * rhorho.BoostVector());
0889 TLorentzVector pi_rhominusb = pi_rhominus;
0890 pi_rhominusb.Boost(-1 * rhorho.BoostVector());
0891 TLorentzVector pi0_rhominusb = pi0_rhominus;
0892 pi0_rhominusb.Boost(-1 * rhorho.BoostVector());
0893
0894
0895 TVector3 n_plus = pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
0896 TVector3 n_minus = pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
0897
0898
0899 double Acoplanarity = acos(n_plus.Dot(n_minus) / (n_plus.Mag() * n_minus.Mag()));
0900 if (pi_rhominusb.Vect().Dot(n_plus) > 0) {
0901 Acoplanarity *= -1;
0902 Acoplanarity += 2 * TMath::Pi();
0903 }
0904
0905
0906 pi_rhoplus.Boost(-1 * taup.BoostVector());
0907 pi0_rhoplus.Boost(-1 * taup.BoostVector());
0908 pi_rhominus.Boost(-1 * taum.BoostVector());
0909 pi0_rhominus.Boost(-1 * taum.BoostVector());
0910
0911
0912 double y1 = (pi_rhoplus.E() - pi0_rhoplus.E()) / (pi_rhoplus.E() + pi0_rhoplus.E());
0913 double y2 = (pi_rhominus.E() - pi0_rhominus.E()) / (pi_rhominus.E() + pi0_rhominus.E());
0914
0915
0916 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0 && y1 * y2 < 0)
0917 TauSpinEffectsH_rhorhoAcoplanarityminus->Fill(Acoplanarity, weight);
0918 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0 && y1 * y2 > 0)
0919 TauSpinEffectsH_rhorhoAcoplanarityplus->Fill(Acoplanarity, weight);
0920 }
0921 if (haspi_minus && haspi_plus) {
0922 TLorentzVector tauporig = taup;
0923 TLorentzVector taumorig = taum;
0924
0925
0926 pi_plus.Boost(-1 * Zboson.BoostVector());
0927 pi_minus.Boost(-1 * Zboson.BoostVector());
0928
0929 taup.Boost(-1 * Zboson.BoostVector());
0930 taum.Boost(-1 * Zboson.BoostVector());
0931
0932 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0) {
0933 TauSpinEffectsH_pipiAcollinearity->Fill(
0934 acos(pi_plus.Vect().Dot(pi_minus.Vect()) / (pi_plus.P() * pi_minus.P())));
0935 TauSpinEffectsH_pipiAcollinearityzoom->Fill(
0936 acos(pi_plus.Vect().Dot(pi_minus.Vect()) / (pi_plus.P() * pi_minus.P())));
0937 }
0938
0939 double proj_m = taum.Vect().Dot(pi_minus.Vect()) / (taum.P() * taum.P());
0940 double proj_p = taup.Vect().Dot(pi_plus.Vect()) / (taup.P() * taup.P());
0941 TVector3 Tau_m = taum.Vect();
0942 TVector3 Tau_p = taup.Vect();
0943 Tau_m *= proj_m;
0944 Tau_p *= proj_p;
0945 TVector3 Pit_m = pi_minus.Vect() - Tau_m;
0946 TVector3 Pit_p = pi_plus.Vect() - Tau_p;
0947
0948 double Acoplanarity = acos(Pit_m.Dot(Pit_p) / (Pit_p.Mag() * Pit_m.Mag()));
0949 TVector3 n = Pit_p.Cross(Pit_m);
0950 if (n.Dot(Tau_m) / Tau_m.Mag() > 0) {
0951 Acoplanarity *= -1;
0952 Acoplanarity += 2 * TMath::Pi();
0953 }
0954
0955 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0956 TauSpinEffectsH_pipiAcoplanarity->Fill(Acoplanarity, weight);
0957 taup = tauporig;
0958 taum = taumorig;
0959 }
0960 if (nSinglePionDecays == 2 && tautau.M() != 0) {
0961 for (int i = 0; i < zsbins; i++) {
0962 double zslow = ((double)i) * (zsmax - zsmin) / ((double)zsbins) + zsmin;
0963 double zsup = ((double)i + 1) * (zsmax - zsmin) / ((double)zsbins) + zsmin;
0964 double aup = Zstoa(zsup), alow = Zstoa(zslow);
0965 if (x2 - x1 > alow && x2 - x1 < aup) {
0966 double zs = (zsup + zslow) / 2;
0967 if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0968 TauSpinEffectsZ_Zs->Fill(zs, weight);
0969 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0970 TauSpinEffectsH_Zs->Fill(zs, weight);
0971 break;
0972 }
0973 }
0974 if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0975 TauSpinEffectsZ_MVis->Fill(pipi.M() / tautau.M(), weight);
0976 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0977 TauSpinEffectsH_MVis->Fill(pipi.M() / tautau.M(), weight);
0978
0979 if (x1 != 0) {
0980 const std::vector<const reco::GenParticle *> m = GetMothers(boson);
0981 int q(0), qbar(0);
0982 TLorentzVector Z(0, 0, 0, 0);
0983 for (unsigned int i = 0; i < m.size(); i++) {
0984 if (m.at(i)->pdgId() == PdtPdgMini::d || m.at(i)->pdgId() == PdtPdgMini::u) {
0985 q++;
0986 }
0987 if (m.at(i)->pdgId() == PdtPdgMini::anti_d || m.at(i)->pdgId() == PdtPdgMini::anti_u) {
0988 qbar++;
0989 }
0990 }
0991 if (q == 1 && qbar == 1) {
0992 if (taum.Vect().Dot(Zboson.Vect()) / (Zboson.P() * taum.P()) > 0) {
0993 if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0994 TauSpinEffectsZ_Xf->Fill(x1, weight);
0995 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0996 TauSpinEffectsH_Xf->Fill(x1, weight);
0997 } else {
0998 if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0999 TauSpinEffectsZ_Xb->Fill(x1, weight);
1000 if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
1001 TauSpinEffectsH_Xb->Fill(x1, weight);
1002 }
1003 }
1004 }
1005 }
1006 }
1007 }
1008
1009 double TauValidation::Zstoa(double zs) {
1010 double a = 1 - sqrt(fabs(1.0 - 2 * fabs(zs)));
1011 if (zs < 0) {
1012 a *= -1.0;
1013 }
1014 return a;
1015 }
1016
1017 double TauValidation::leadingPionMomentum(const reco::GenParticle *tau, double weight) {
1018 return leadingPionP4(tau).P();
1019 }
1020
1021 TLorentzVector TauValidation::leadingPionP4(const reco::GenParticle *tau) {
1022 TLorentzVector p4(0, 0, 0, 0);
1023 for (unsigned int i = 0; i < tau->numberOfDaughters(); i++) {
1024 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(tau->daughter(i));
1025 int pid = dau->pdgId();
1026 if (abs(pid) == 15)
1027 return leadingPionP4(dau);
1028 if (!(abs(pid) == 211 || abs(pid) == 13 || abs(pid) == 11))
1029 continue;
1030 if (dau->p() > p4.P())
1031 p4 = TLorentzVector(dau->px(), dau->py(), dau->pz(), dau->energy());
1032 }
1033 return p4;
1034 }
1035
1036 TLorentzVector TauValidation::motherP4(const reco::GenParticle *tau) {
1037 const reco::GenParticle *m = GetMother(tau);
1038 return TLorentzVector(m->px(), m->py(), m->pz(), m->energy());
1039 }
1040
1041 double TauValidation::visibleTauEnergy(const reco::GenParticle *tau) {
1042 TLorentzVector p4(tau->px(), tau->py(), tau->pz(), tau->energy());
1043 for (unsigned int i = 0; i < tau->numberOfDaughters(); i++) {
1044 const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(tau->daughter(i));
1045 int pid = dau->pdgId();
1046 if (abs(pid) == 15)
1047 return visibleTauEnergy(dau);
1048 if (abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
1049 p4 -= TLorentzVector(dau->px(), dau->py(), dau->pz(), dau->energy());
1050 }
1051 }
1052 return p4.E();
1053 }
1054
1055 void TauValidation::photons(const reco::GenParticle *tau, double weight) {
1056
1057 std::vector<const reco::GenParticle *> TauList;
1058 findTauList(tau, TauList);
1059
1060
1061 bool passedW = false;
1062 std::vector<const reco::GenParticle *> ListofFSR;
1063 ListofFSR.clear();
1064 std::vector<const reco::GenParticle *> ListofBrem;
1065 ListofBrem.clear();
1066 std::vector<const reco::GenParticle *> FSR_photos;
1067 FSR_photos.clear();
1068 double BosonScale(1);
1069 if (!TauList.empty()) {
1070 TauValidation::findFSRandBrem(TauList.at(0), passedW, ListofFSR, ListofBrem);
1071 TauValidation::FindPhotosFSR(TauList.at(0), FSR_photos, BosonScale);
1072
1073
1074 TauBremPhotonsN->Fill(ListofBrem.size(), weight);
1075 double photonPtSum = 0;
1076 for (unsigned int i = 0; i < ListofBrem.size(); i++) {
1077 photonPtSum += ListofBrem.at(i)->pt();
1078 TauBremPhotonsPt->Fill(ListofBrem.at(i)->pt(), weight);
1079 }
1080 TauBremPhotonsPtSum->Fill(photonPtSum, weight);
1081
1082
1083 if (BosonScale != 0) {
1084 TauFSRPhotonsN->Fill(ListofFSR.size(), weight);
1085 photonPtSum = 0;
1086 for (unsigned int i = 0; i < ListofFSR.size(); i++) {
1087 photonPtSum += ListofFSR.at(i)->pt();
1088 TauFSRPhotonsPt->Fill(ListofFSR.at(i)->pt(), weight);
1089 }
1090 double FSR_photosSum(0);
1091 for (unsigned int i = 0; i < FSR_photos.size(); i++) {
1092 FSR_photosSum += FSR_photos.at(i)->pt();
1093 TauFSRPhotonsPt->Fill(FSR_photos.at(i)->pt() / BosonScale, weight * BosonScale);
1094 }
1095 TauFSRPhotonsPtSum->Fill(photonPtSum + FSR_photosSum / BosonScale, weight);
1096 }
1097 }
1098 }