test/ntuples/HGCalTriggerNtupleGenTau.cc

0001 #include <vector>
0002 #include "DataFormats/Math/interface/LorentzVector.h"
0003 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0004 #include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h"
0005 #include "DataFormats/Math/interface/LorentzVector.h"
0006
0007 typedef math::XYZTLorentzVector LorentzVector;
0008
0009 class HGCalTriggerNtupleGenTau : public HGCalTriggerNtupleBase {
0010 public:
0011   HGCalTriggerNtupleGenTau(const edm::ParameterSet&);
0012
0013   void initialize(TTree&, const edm::ParameterSet&, edm::ConsumesCollector&&) final;
0014   void fill(const edm::Event&, const HGCalTriggerNtupleEventSetup&) final;
0015
0016 private:
0017   void clear() final;
0018
0019   bool isGoodTau(const reco::GenParticle& candidate) const;
0020   bool isStableLepton(const reco::GenParticle& daughter) const;
0021   bool isElectron(const reco::GenParticle& daughter) const;
0022   bool isMuon(const reco::GenParticle& daughter) const;
0023   bool isChargedHadron(const reco::GenParticle& daughter) const;
0024   bool isChargedHadronFromResonance(const reco::GenParticle& daughter) const;
0025   bool isNeutralPion(const reco::GenParticle& daughter) const;
0026   bool isNeutralPionFromResonance(const reco::GenParticle& daughter) const;
0027   bool isIntermediateResonance(const reco::GenParticle& daughter) const;
0028   bool isGamma(const reco::GenParticle& daughter) const;
0029   bool isStableNeutralHadron(const reco::GenParticle& daughter) const;
0030
0031   edm::EDGetToken gen_token_;
0032   bool isPythia8generator_;
0033
0034   std::vector<float> gentau_pt_;
0035   std::vector<float> gentau_eta_;
0036   std::vector<float> gentau_phi_;
0037   std::vector<float> gentau_energy_;
0038   std::vector<float> gentau_mass_;
0039
0040   std::vector<float> gentau_vis_pt_;
0041   std::vector<float> gentau_vis_eta_;
0042   std::vector<float> gentau_vis_phi_;
0043   std::vector<float> gentau_vis_energy_;
0044   std::vector<float> gentau_vis_mass_;
0045   std::vector<int> gentau_decayMode_;
0046   std::vector<int> gentau_totNproducts_;
0047   std::vector<int> gentau_totNgamma_;
0048   std::vector<int> gentau_totNpiZero_;
0049   std::vector<int> gentau_totNcharged_;
0050
0051   std::vector<std::vector<float> > gentau_products_pt_;
0052   std::vector<std::vector<float> > gentau_products_eta_;
0053   std::vector<std::vector<float> > gentau_products_phi_;
0054   std::vector<std::vector<float> > gentau_products_energy_;
0055   std::vector<std::vector<float> > gentau_products_mass_;
0056   std::vector<std::vector<int> > gentau_products_id_;
0057 };
0058
0059 DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleGenTau, "HGCalTriggerNtupleGenTau");
0060
0061 HGCalTriggerNtupleGenTau::HGCalTriggerNtupleGenTau(const edm::ParameterSet& conf) : HGCalTriggerNtupleBase(conf) {
0062   accessEventSetup_ = false;
0063 }
0064
0065 void HGCalTriggerNtupleGenTau::initialize(TTree& tree,
0066                                           const edm::ParameterSet& conf,
0067                                           edm::ConsumesCollector&& collector) {
0068   gen_token_ = collector.consumes<reco::GenParticleCollection>(conf.getParameter<edm::InputTag>("GenParticles"));
0069   isPythia8generator_ = conf.getParameter<bool>("isPythia8");
0070
0071   tree.Branch("gentau_pt", &gentau_pt_);
0072   tree.Branch("gentau_eta", &gentau_eta_);
0073   tree.Branch("gentau_phi", &gentau_phi_);
0074   tree.Branch("gentau_energy", &gentau_energy_);
0075   tree.Branch("gentau_mass", &gentau_mass_);
0076   tree.Branch("gentau_vis_pt", &gentau_vis_pt_);
0077   tree.Branch("gentau_vis_eta", &gentau_vis_eta_);
0078   tree.Branch("gentau_vis_phi", &gentau_vis_phi_);
0079   tree.Branch("gentau_vis_energy", &gentau_vis_energy_);
0080   tree.Branch("gentau_vis_mass", &gentau_vis_mass_);
0081   tree.Branch("gentau_products_pt", &gentau_products_pt_);
0082   tree.Branch("gentau_products_eta", &gentau_products_eta_);
0083   tree.Branch("gentau_products_phi", &gentau_products_phi_);
0084   tree.Branch("gentau_products_energy", &gentau_products_energy_);
0085   tree.Branch("gentau_products_mass", &gentau_products_mass_);
0086   tree.Branch("gentau_products_id", &gentau_products_id_);
0087   tree.Branch("gentau_decayMode", &gentau_decayMode_);
0088   tree.Branch("gentau_totNproducts", &gentau_totNproducts_);
0089   tree.Branch("gentau_totNgamma", &gentau_totNgamma_);
0090   tree.Branch("gentau_totNpiZero", &gentau_totNpiZero_);
0091   tree.Branch("gentau_totNcharged", &gentau_totNcharged_);
0092 }
0093
0094 bool HGCalTriggerNtupleGenTau::isGoodTau(const reco::GenParticle& candidate) const {
0095   return (std::abs(candidate.pdgId()) == 15 && candidate.status() == 2);
0096 }
0097
0098 bool HGCalTriggerNtupleGenTau::isChargedHadron(const reco::GenParticle& candidate) const {
0099   return ((std::abs(candidate.pdgId()) == 211 || std::abs(candidate.pdgId()) == 321) && candidate.status() == 1 &&
0100           candidate.isDirectPromptTauDecayProductFinalState() && candidate.isLastCopy());
0101 }
0102
0103 bool HGCalTriggerNtupleGenTau::isChargedHadronFromResonance(const reco::GenParticle& candidate) const {
0104   return ((std::abs(candidate.pdgId()) == 211 || std::abs(candidate.pdgId()) == 321) && candidate.status() == 1 &&
0105           candidate.isLastCopy());
0106 }
0107
0108 bool HGCalTriggerNtupleGenTau::isStableLepton(const reco::GenParticle& candidate) const {
0109   return ((std::abs(candidate.pdgId()) == 11 || std::abs(candidate.pdgId()) == 13) && candidate.status() == 1 &&
0110           candidate.isDirectPromptTauDecayProductFinalState() && candidate.isLastCopy());
0111 }
0112
0113 bool HGCalTriggerNtupleGenTau::isElectron(const reco::GenParticle& candidate) const {
0114   return (std::abs(candidate.pdgId()) == 11 && candidate.isDirectPromptTauDecayProductFinalState() &&
0115           candidate.isLastCopy());
0116 }
0117
0118 bool HGCalTriggerNtupleGenTau::isMuon(const reco::GenParticle& candidate) const {
0119   return (std::abs(candidate.pdgId()) == 13 && candidate.isDirectPromptTauDecayProductFinalState() &&
0120           candidate.isLastCopy());
0121 }
0122
0123 bool HGCalTriggerNtupleGenTau::isNeutralPion(const reco::GenParticle& candidate) const {
0124   return (std::abs(candidate.pdgId()) == 111 && candidate.status() == 2 &&
0125           candidate.statusFlags().isTauDecayProduct() && !candidate.isDirectPromptTauDecayProductFinalState());
0126 }
0127
0128 bool HGCalTriggerNtupleGenTau::isNeutralPionFromResonance(const reco::GenParticle& candidate) const {
0129   return (std::abs(candidate.pdgId()) == 111 && candidate.status() == 2 && candidate.statusFlags().isTauDecayProduct());
0130 }
0131
0132 bool HGCalTriggerNtupleGenTau::isGamma(const reco::GenParticle& candidate) const {
0133   return (std::abs(candidate.pdgId()) == 22 && candidate.status() == 1 && candidate.statusFlags().isTauDecayProduct() &&
0134           !candidate.isDirectPromptTauDecayProductFinalState() && candidate.isLastCopy());
0135 }
0136
0137 bool HGCalTriggerNtupleGenTau::isIntermediateResonance(const reco::GenParticle& candidate) const {
0138   return ((std::abs(candidate.pdgId()) == 213 || std::abs(candidate.pdgId()) == 20213 ||
0139            std::abs(candidate.pdgId()) == 24) &&
0140           candidate.status() == 2);
0141 }
0142
0143 bool HGCalTriggerNtupleGenTau::isStableNeutralHadron(const reco::GenParticle& candidate) const {
0144   return (!(std::abs(candidate.pdgId()) > 10 && std::abs(candidate.pdgId()) < 17) && !isChargedHadron(candidate) &&
0145           candidate.status() == 1);
0146 }
0147
0148 void HGCalTriggerNtupleGenTau::fill(const edm::Event& e, const HGCalTriggerNtupleEventSetup& es) {
0149   edm::Handle<reco::GenParticleCollection> gen_particles_h;
0150   e.getByToken(gen_token_, gen_particles_h);
0151   const reco::GenParticleCollection& gen_particles = *gen_particles_h;
0152
0153   clear();
0154
0155   for (const auto& particle : gen_particles) {
0156     /* select good taus */
0157     if (isGoodTau(particle)) {
0158       LorentzVector tau_p4vis(0., 0., 0., 0.);
0159       gentau_pt_.emplace_back(particle.pt());
0160       gentau_eta_.emplace_back(particle.eta());
0161       gentau_phi_.emplace_back(particle.phi());
0162       gentau_energy_.emplace_back(particle.energy());
0163       gentau_mass_.emplace_back(particle.mass());
0164
0165       int n_pi = 0;
0166       int n_piZero = 0;
0167       int n_gamma = 0;
0168       int n_ele = 0;
0169       int n_mu = 0;
0170
0171       std::vector<float> tau_products_pt;
0172       std::vector<float> tau_products_eta;
0173       std::vector<float> tau_products_phi;
0174       std::vector<float> tau_products_energy;
0175       std::vector<float> tau_products_mass;
0176       std::vector<int> tau_products_id;
0177
0178       /* loop over tau daughters */
0179       const reco::GenParticleRefVector& daughters = particle.daughterRefVector();
0180
0181       for (const auto& daughter : daughters) {
0182         reco::GenParticleRefVector finalProds;
0183
0184         if (isStableLepton(*daughter)) {
0185           if (isElectron(*daughter)) {
0186             n_ele++;
0187           } else if (isMuon(*daughter)) {
0188             n_mu++;
0189           }
0190           tau_p4vis += (daughter->p4());
0191           finalProds.push_back(daughter);
0192         }
0193
0194         else if (isChargedHadron(*daughter)) {
0195           n_pi++;
0196           tau_p4vis += (daughter->p4());
0197           finalProds.push_back(daughter);
0198         }
0199
0200         else if (isNeutralPion(*daughter)) {
0201           n_piZero++;
0202           const reco::GenParticleRefVector& granddaughters = daughter->daughterRefVector();
0203           for (const auto& granddaughter : granddaughters) {
0204             if (isGamma(*granddaughter)) {
0205               n_gamma++;
0206               tau_p4vis += (granddaughter->p4());
0207               finalProds.push_back(granddaughter);
0208             }
0209           }
0210         }
0211
0212         else if (isStableNeutralHadron(*daughter)) {
0213           tau_p4vis += (daughter->p4());
0214           finalProds.push_back(daughter);
0215         }
0216
0217         else {
0218           const reco::GenParticleRefVector& granddaughters = daughter->daughterRefVector();
0219
0220           for (const auto& granddaughter : granddaughters) {
0221             if (isStableNeutralHadron(*granddaughter)) {
0222               tau_p4vis += (granddaughter->p4());
0223               finalProds.push_back(granddaughter);
0224             }
0225           }
0226         }
0227
0228         /* Here the selection of the decay product according to the Pythia6 decayTree */
0229         if (!isPythia8generator_) {
0230           if (isIntermediateResonance(*daughter)) {
0231             const reco::GenParticleRefVector& grandaughters = daughter->daughterRefVector();
0232             for (const auto& grandaughter : grandaughters) {
0233               if (isChargedHadron(*grandaughter) || isChargedHadronFromResonance(*grandaughter)) {
0234                 n_pi++;
0235                 tau_p4vis += (grandaughter->p4());
0236                 finalProds.push_back(daughter);
0237               } else if (isNeutralPion(*grandaughter) || isNeutralPionFromResonance(*grandaughter)) {
0238                 n_piZero++;
0239                 const reco::GenParticleRefVector& descendants = grandaughter->daughterRefVector();
0240                 for (const auto& descendant : descendants) {
0241                   if (isGamma(*descendant)) {
0242                     n_gamma++;
0243                     tau_p4vis += (descendant->p4());
0244                     finalProds.push_back(daughter);
0245                   }
0246                 }
0247               }
0248             }
0249           }
0250         }
0251
0252         /* Fill daughter informations */
0253         for (const auto& prod : finalProds) {
0254           tau_products_pt.emplace_back(prod->pt());
0255           tau_products_eta.emplace_back(prod->eta());
0256           tau_products_phi.emplace_back(prod->phi());
0257           tau_products_energy.emplace_back(prod->energy());
0258           tau_products_mass.emplace_back(prod->mass());
0259           tau_products_id.emplace_back(prod->pdgId());
0260         }
0261       }
0262
0263       /* assign the tau-variables */
0264       gentau_vis_pt_.emplace_back(tau_p4vis.Pt());
0265       gentau_vis_eta_.emplace_back(tau_p4vis.Eta());
0266       gentau_vis_phi_.emplace_back(tau_p4vis.Phi());
0267       gentau_vis_energy_.emplace_back(tau_p4vis.E());
0268       gentau_vis_mass_.emplace_back(tau_p4vis.M());
0269       gentau_totNproducts_.emplace_back(n_pi + n_gamma);
0270       gentau_totNgamma_.emplace_back(n_gamma);
0271       gentau_totNpiZero_.emplace_back(n_piZero);
0272       gentau_totNcharged_.emplace_back(n_pi);
0273
0274       gentau_products_pt_.emplace_back(tau_products_pt);
0275       gentau_products_eta_.emplace_back(tau_products_eta);
0276       gentau_products_phi_.emplace_back(tau_products_phi);
0277       gentau_products_energy_.emplace_back(tau_products_energy);
0278       gentau_products_mass_.emplace_back(tau_products_mass);
0279       gentau_products_id_.emplace_back(tau_products_id);
0280
0281       /* leptonic tau decays */
0282       if (n_pi == 0 && n_piZero == 0 && n_ele == 1) {
0283         gentau_decayMode_.emplace_back(11);
0284       } else if (n_pi == 0 && n_piZero == 0 && n_mu == 1) {
0285         gentau_decayMode_.emplace_back(13);
0286       }
0287       /* 1-prong */
0288       else if (n_pi == 1 && n_piZero == 0) {
0289         gentau_decayMode_.emplace_back(0);
0290       }
0291       /* 1-prong + pi0s */
0292       else if (n_pi == 1 && n_piZero >= 1) {
0293         gentau_decayMode_.emplace_back(1);
0294       }
0295       /* 3-prongs */
0296       else if (n_pi == 3 && n_piZero == 0) {
0297         gentau_decayMode_.emplace_back(4);
0298       }
0299       /* 3-prongs + pi0s */
0300       else if (n_pi == 3 && n_piZero >= 1) {
0301         gentau_decayMode_.emplace_back(5);
0302       }
0303       /* other decays */
0304       else {
0305         gentau_decayMode_.emplace_back(-1);
0306       }
0307     }
0308   }
0309 }
0310
0311 void HGCalTriggerNtupleGenTau::clear() {
0312   gentau_pt_.clear();
0313   gentau_eta_.clear();
0314   gentau_phi_.clear();
0315   gentau_energy_.clear();
0316   gentau_mass_.clear();
0317   gentau_decayMode_.clear();
0318   gentau_vis_pt_.clear();
0319   gentau_vis_eta_.clear();
0320   gentau_vis_phi_.clear();
0321   gentau_vis_energy_.clear();
0322   gentau_vis_mass_.clear();
0323   gentau_totNproducts_.clear();
0324   gentau_totNgamma_.clear();
0325   gentau_totNcharged_.clear();
0326   gentau_products_pt_.clear();
0327   gentau_products_eta_.clear();
0328   gentau_products_phi_.clear();
0329   gentau_products_energy_.clear();
0330   gentau_products_mass_.clear();
0331   gentau_products_id_.clear();
0332 }