Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-20 10:42:34

0001 /*class TauValidation
0002  *  
0003  *  Class to fill dqm monitor elements from existing EDM file
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     :  //  wmanager_(iPSet,consumesCollector())
0020       genparticleCollection_(iPSet.getParameter<edm::InputTag>("genparticleCollection")),
0021       NMODEID(TauDecay::NMODEID - 1),  // fortran to C++ index
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   ///Setting the DQM top directories
0037   DQMHelper dqm(&i);
0038   i.setCurrentFolder("Generator/Tau");
0039   // Number of analyzed events
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   //Kinematics
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   ///Gathering the reco::GenParticleCollection information
0373   edm::Handle<reco::GenParticleCollection> genParticles;
0374   iEvent.getByToken(genparticleCollectionToken_, genParticles);
0375 
0376   double weight = 1.0;  //=   wmanager_.weight(iEvent);
0377   //////////////////////////////////////////////
0378   // find taus
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) {  // exclude B, D and other non-signal decay modes
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           // Adding MODEID and Mass information
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) /*cm->m*/, 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 }  //analyze
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   // 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.
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) {  // remove pi0 and eta decays
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;  // W
0552   if (abs(mother_pid) == 23)
0553     BosonScale = 2.0;  // Z;
0554   if (abs(mother_pid) == 22)
0555     BosonScale = 2.0;  // gamma;
0556   if (abs(mother_pid) == 25)
0557     BosonScale = 2.0;  // HSM;
0558   if (abs(mother_pid) == 35)
0559     BosonScale = 2.0;  // H0;
0560   if (abs(mother_pid) == 36)
0561     BosonScale = 2.0;  // A0;
0562   if (abs(mother_pid) == 37)
0563     BosonScale = 1.0;  //Hpm;
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   // resonances
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   // final state products
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   // polarization only for 1-prong hadronic taus with no neutral pions
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) {  // only for pi2pi0 for now
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), nSingleMuonDecays(0), nSingleElectronDecays(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             if (jak_id == TauDecay::MODE_MUON)
0776               nSingleMuonDecays++;
0777             if (jak_id == TauDecay::MODE_ELECTRON)
0778               nSingleElectronDecays++;
0779             TLorentzVector LVtau(dau->px(), dau->py(), dau->pz(), dau->energy());
0780             tautau += LVtau;
0781             TLorentzVector LVpi = leadingPionP4(dau);
0782             pipi += LVpi;
0783             const HepPDT::ParticleData *pd = fPDGTable->particle(dau->pdgId());
0784             int charge = (int)pd->charge();
0785             LVtau.Boost(-1 * Zboson.BoostVector());
0786             LVpi.Boost(-1 * Zboson.BoostVector());
0787 
0788             if (jak_id == TauDecay::MODE_MUON) {
0789               if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0790                 TauSpinEffectsZ_muX->Fill(LVpi.P() / LVtau.E(), weight);
0791               if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0792                 TauSpinEffectsH_muX->Fill(LVpi.P() / LVtau.E(), weight);
0793             }
0794             if (jak_id == TauDecay::MODE_ELECTRON) {
0795               if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0796                 TauSpinEffectsZ_eX->Fill(LVpi.P() / LVtau.E(), weight);
0797               if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0798                 TauSpinEffectsH_eX->Fill(LVpi.P() / LVtau.E(), weight);
0799             }
0800 
0801             if (jak_id == TauDecay::MODE_PION) {
0802               if (abs(boson->pdgId()) == PdtPdgMini::Z0) {
0803                 TauSpinEffectsZ_X->Fill(LVpi.P() / LVtau.E(), weight);
0804                 if (50.0 < Zboson.M() && Zboson.M() < 75.0)
0805                   TauSpinEffectsZ_X50to75->Fill(LVpi.P() / LVtau.E(), weight);
0806                 if (75.0 < Zboson.M() && Zboson.M() < 88.0)
0807                   TauSpinEffectsZ_X75to88->Fill(LVpi.P() / LVtau.E(), weight);
0808                 if (88.0 < Zboson.M() && Zboson.M() < 100.0)
0809                   TauSpinEffectsZ_X88to100->Fill(LVpi.P() / LVtau.E(), weight);
0810                 if (100.0 < Zboson.M() && Zboson.M() < 120.0)
0811                   TauSpinEffectsZ_X100to120->Fill(LVpi.P() / LVtau.E(), weight);
0812                 if (120.0 < Zboson.M())
0813                   TauSpinEffectsZ_X120UP->Fill(LVpi.P() / LVtau.E(), weight);
0814               }
0815               if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0816                 TauSpinEffectsH_X->Fill(LVpi.P() / LVtau.E(), weight);
0817             }
0818             if (charge < 0) {
0819               x1 = LVpi.P() / LVtau.E();
0820               taum = LVtau;
0821             } else {
0822               x2 = LVpi.P() / LVtau.E();
0823             }
0824           }
0825           TLorentzVector LVtau(dau->px(), dau->py(), dau->pz(), dau->energy());
0826           if (pid == 15)
0827             taum = LVtau;
0828           if (pid == -15)
0829             taup = LVtau;
0830           if (jak_id == TauDecay::MODE_PIPI0) {
0831             for (unsigned int i = 0; i < part.size(); i++) {
0832               int pid_d = part.at(i)->pdgId();
0833               if (abs(pid_d) == 211 || abs(pid_d) == 111) {
0834                 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0835                 if (pid == 15) {
0836                   hasrho_minus = true;
0837                   if (pid_d == -211) {
0838                     pi_rhominus = LV;
0839                   }
0840                   if (abs(pid_d) == 111) {
0841                     pi0_rhominus = LV;
0842                   }
0843                 }
0844                 if (pid == -15) {
0845                   hasrho_plus = true;
0846                   if (pid_d == 211) {
0847                     pi_rhoplus = LV;
0848                   }
0849                   if (abs(pid_d) == 111) {
0850                     pi0_rhoplus = LV;
0851                   }
0852                 }
0853               }
0854             }
0855           }
0856           if (jak_id == TauDecay::MODE_PION) {
0857             for (unsigned int i = 0; i < part.size(); i++) {
0858               int pid_d = part.at(i)->pdgId();
0859               if (abs(pid_d) == 211) {
0860                 TLorentzVector LV(part.at(i)->px(), part.at(i)->py(), part.at(i)->pz(), part.at(i)->energy());
0861                 if (pid == 15) {
0862                   haspi_minus = true;
0863                   if (pid_d == -211) {
0864                     pi_minus = LV;
0865                   }
0866                 }
0867                 if (pid == -15) {
0868                   haspi_plus = true;
0869                   if (pid_d == 211) {
0870                     pi_plus = LV;
0871                   }
0872                 }
0873               }
0874             }
0875           }
0876         }
0877       }
0878     }
0879     if (hasrho_minus && hasrho_plus) {
0880       //compute rhorho
0881       rho_minus = pi_rhominus;
0882       rho_minus += pi0_rhominus;
0883       rho_plus = pi_rhoplus;
0884       rho_plus += pi0_rhoplus;
0885       TLorentzVector rhorho = rho_minus;
0886       rhorho += rho_plus;
0887 
0888       // boost to rhorho cm
0889       TLorentzVector pi_rhoplusb = pi_rhoplus;
0890       pi_rhoplusb.Boost(-1 * rhorho.BoostVector());
0891       TLorentzVector pi0_rhoplusb = pi0_rhoplus;
0892       pi0_rhoplusb.Boost(-1 * rhorho.BoostVector());
0893       TLorentzVector pi_rhominusb = pi_rhominus;
0894       pi_rhominusb.Boost(-1 * rhorho.BoostVector());
0895       TLorentzVector pi0_rhominusb = pi0_rhominus;
0896       pi0_rhominusb.Boost(-1 * rhorho.BoostVector());
0897 
0898       // compute n+/-
0899       TVector3 n_plus = pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
0900       TVector3 n_minus = pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
0901 
0902       // compute the acoplanarity
0903       double Acoplanarity = acos(n_plus.Dot(n_minus) / (n_plus.Mag() * n_minus.Mag()));
0904       if (pi_rhominusb.Vect().Dot(n_plus) > 0) {
0905         Acoplanarity *= -1;
0906         Acoplanarity += 2 * TMath::Pi();
0907       }
0908 
0909       // now boost to tau frame
0910       pi_rhoplus.Boost(-1 * taup.BoostVector());
0911       pi0_rhoplus.Boost(-1 * taup.BoostVector());
0912       pi_rhominus.Boost(-1 * taum.BoostVector());
0913       pi0_rhominus.Boost(-1 * taum.BoostVector());
0914 
0915       // compute y1 and y2
0916       double y1 = (pi_rhoplus.E() - pi0_rhoplus.E()) / (pi_rhoplus.E() + pi0_rhoplus.E());
0917       double y2 = (pi_rhominus.E() - pi0_rhominus.E()) / (pi_rhominus.E() + pi0_rhominus.E());
0918 
0919       // fill histograms
0920       if (abs(boson->pdgId()) == PdtPdgMini::Higgs0 && y1 * y2 < 0)
0921         TauSpinEffectsH_rhorhoAcoplanarityminus->Fill(Acoplanarity, weight);
0922       if (abs(boson->pdgId()) == PdtPdgMini::Higgs0 && y1 * y2 > 0)
0923         TauSpinEffectsH_rhorhoAcoplanarityplus->Fill(Acoplanarity, weight);
0924     }
0925     if (haspi_minus && haspi_plus) {
0926       TLorentzVector tauporig = taup;
0927       TLorentzVector taumorig = taum;
0928 
0929       // now boost to Higgs frame
0930       pi_plus.Boost(-1 * Zboson.BoostVector());
0931       pi_minus.Boost(-1 * Zboson.BoostVector());
0932 
0933       taup.Boost(-1 * Zboson.BoostVector());
0934       taum.Boost(-1 * Zboson.BoostVector());
0935 
0936       if (abs(boson->pdgId()) == PdtPdgMini::Higgs0) {
0937         TauSpinEffectsH_pipiAcollinearity->Fill(
0938             acos(pi_plus.Vect().Dot(pi_minus.Vect()) / (pi_plus.P() * pi_minus.P())));
0939         TauSpinEffectsH_pipiAcollinearityzoom->Fill(
0940             acos(pi_plus.Vect().Dot(pi_minus.Vect()) / (pi_plus.P() * pi_minus.P())));
0941       }
0942 
0943       double proj_m = taum.Vect().Dot(pi_minus.Vect()) / (taum.P() * taum.P());
0944       double proj_p = taup.Vect().Dot(pi_plus.Vect()) / (taup.P() * taup.P());
0945       TVector3 Tau_m = taum.Vect();
0946       TVector3 Tau_p = taup.Vect();
0947       Tau_m *= proj_m;
0948       Tau_p *= proj_p;
0949       TVector3 Pit_m = pi_minus.Vect() - Tau_m;
0950       TVector3 Pit_p = pi_plus.Vect() - Tau_p;
0951 
0952       double Acoplanarity = acos(Pit_m.Dot(Pit_p) / (Pit_p.Mag() * Pit_m.Mag()));
0953       TVector3 n = Pit_p.Cross(Pit_m);
0954       if (n.Dot(Tau_m) / Tau_m.Mag() > 0) {
0955         Acoplanarity *= -1;
0956         Acoplanarity += 2 * TMath::Pi();
0957       }
0958       // fill histograms
0959       if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0960         TauSpinEffectsH_pipiAcoplanarity->Fill(Acoplanarity, weight);
0961       taup = tauporig;
0962       taum = taumorig;
0963     }
0964     if (nSinglePionDecays == 2 && tautau.M() != 0) {
0965       for (int i = 0; i < zsbins; i++) {
0966         double zslow = ((double)i) * (zsmax - zsmin) / ((double)zsbins) + zsmin;
0967         double zsup = ((double)i + 1) * (zsmax - zsmin) / ((double)zsbins) + zsmin;
0968         double aup = Zstoa(zsup), alow = Zstoa(zslow);
0969         if (x2 - x1 > alow && x2 - x1 < aup) {
0970           double zs = (zsup + zslow) / 2;
0971           if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0972             TauSpinEffectsZ_Zs->Fill(zs, weight);
0973           if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0974             TauSpinEffectsH_Zs->Fill(zs, weight);
0975           break;
0976         }
0977       }
0978       if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0979         TauSpinEffectsZ_MVis->Fill(pipi.M() / tautau.M(), weight);
0980       if (abs(boson->pdgId()) == PdtPdgMini::Higgs0)
0981         TauSpinEffectsH_MVis->Fill(pipi.M() / tautau.M(), weight);
0982 
0983       if (x1 != 0) {
0984         const std::vector<const reco::GenParticle *> m = GetMothers(boson);
0985         int q(0), qbar(0);
0986         TLorentzVector Z(0, 0, 0, 0);
0987         for (unsigned int i = 0; i < m.size(); i++) {
0988           if (m.at(i)->pdgId() == PdtPdgMini::d || m.at(i)->pdgId() == PdtPdgMini::u) {
0989             q++;
0990           }
0991           if (m.at(i)->pdgId() == PdtPdgMini::anti_d || m.at(i)->pdgId() == PdtPdgMini::anti_u) {
0992             qbar++;
0993           }
0994         }
0995         if (q == 1 && qbar == 1) {  // assume q has largest E (valence vs see quarks)
0996           if (taum.Vect().Dot(Zboson.Vect()) / (Zboson.P() * taum.P()) > 0) {
0997             if (abs(boson->pdgId()) == PdtPdgMini::Z0)
0998               TauSpinEffectsZ_Xf->Fill(x1, weight);
0999             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 }