NanoAOD/plugins/TauSpinnerTableProducer.cc

0001 /**\class TauSpinnerTableProducer
0002
0003  Description: Produces FlatTable with TauSpinner weights for H->tau,tau events
0004
0005  Original Author: D. Winterbottom (IC)
0006  Update and adaptation to NanoAOD: M. Bluj (NCBJ)
0007
0008 */
0009
0010 #include "FWCore/Framework/interface/one/EDProducer.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ParameterSet/interface/allowedValues.h"
0019 #include "DataFormats/Common/interface/View.h"
0020 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0021 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0022
0023 #include "Tauola/Tauola.h"
0024 #include "TauSpinner/SimpleParticle.h"
0025 #include "TauSpinner/tau_reweight_lib.h"
0026
0027 class TauSpinnerTableProducer : public edm::one::EDProducer<edm::one::SharedResources> {
0028 public:
0029   explicit TauSpinnerTableProducer(const edm::ParameterSet &);
0030
0031   void produce(edm::Event &, const edm::EventSetup &) final;
0032   void beginJob() final;
0033   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0034     edm::ParameterSetDescription desc;
0035     desc.add<edm::InputTag>("src")->setComment("input genParticle collection");
0036     desc.add<std::string>("name")->setComment("name of the TauSpinner weights table");
0037     desc.add<std::vector<double>>("theta")->setComment("values of CP-even and CP-odd tau Yukawa mixing angle");
0038     desc.ifValue(edm::ParameterDescription<int>("bosonPdgId", 25, true), edm::allowedValues<int>(25, 35, 36))
0039         ->setComment("boson pdgId, default: 25");  // Allow only neutral Higgs bosons
0040     desc.add<double>("defaultWeight", 1)
0041         ->setComment("default weight stored in case of presence of a tau decay unsupported by TauSpinner");
0042     descriptions.addWithDefaultLabel(desc);
0043   }
0044
0045 private:
0046   void getBosons(edm::RefVector<edm::View<reco::GenParticle>> &bosons, const edm::View<reco::GenParticle> &parts) const;
0047   static reco::GenParticleRef getLastCopy(const reco::GenParticleRef &part);
0048   static void getTaus(reco::GenParticleRefVector &taus, const reco::GenParticle &boson);
0049   static bool getTauDaughters(reco::GenParticleRefVector &tau_daughters, const reco::GenParticle &tau);
0050   TauSpinner::SimpleParticle convertToSimplePart(const reco::GenParticle &input_part) const {
0051     return TauSpinner::SimpleParticle(
0052         input_part.px(), input_part.py(), input_part.pz(), input_part.energy(), input_part.pdgId());
0053   }
0054
0055   static std::vector<std::pair<std::string, double>> nameAndValue(const std::vector<double> &val_vec) {
0056     std::vector<std::pair<std::string, double>> out;
0057     for (auto val : val_vec) {
0058       std::string name = std::to_string(val);
0059       name.erase(name.find_last_not_of('0') + 1, std::string::npos);
0060       name.erase(name.find_last_not_of('.') + 1, std::string::npos);
0061       size_t pos = name.find(".");
0062       if (pos != std::string::npos)
0063         name.replace(pos, 1, "p");
0064       pos = name.find("-");
0065       if (pos != std::string::npos)
0066         name.replace(pos, 1, "minus");
0067       out.push_back(std::make_pair(name, val));
0068     }
0069     return out;
0070   }
0071
0072   void printModuleInfo(edm::ParameterSet const &config) const {
0073     std::cout << std::string(78, '-') << "\n";
0074     std::cout << config.getParameter<std::string>("@module_type") << '/'
0075               << config.getParameter<std::string>("@module_label") << "\n";
0076     std::cout << "Input collection: " << config.getParameter<edm::InputTag>("src").encode() << '\n';
0077     std::cout << "Table name: " << config.getParameter<std::string>("name") << '\n';
0078     std::string thetaStr;
0079     for (const auto &theta : theta_vec_)
0080       thetaStr += theta.first + ",";
0081     std::cout << "Theta: " << thetaStr << std::endl;
0082   }
0083
0084   const edm::EDGetTokenT<edm::View<reco::GenParticle>> genPartsToken_;
0085   const std::string name_;
0086   const std::vector<std::pair<std::string, double>> theta_vec_;
0087   const int bosonPdgId_;
0088   const std::string tauSpinnerPDF_;
0089   const bool ipp_;
0090   const int ipol_;
0091   const int nonSM2_;
0092   const int nonSMN_;
0093   const double cmsE_;
0094   const double default_weight_;
0095 };
0096
0097 TauSpinnerTableProducer::TauSpinnerTableProducer(const edm::ParameterSet &config)
0098     : genPartsToken_(consumes(config.getParameter<edm::InputTag>("src"))),
0099       name_(config.getParameter<std::string>("name")),
0100       theta_vec_(nameAndValue(config.getParameter<std::vector<double>>("theta"))),
0101       bosonPdgId_(config.getParameter<int>("bosonPdgId")),
0102       tauSpinnerPDF_(
0103           "NNPDF31_nnlo_hessian_pdfas"),  // PDF set for TauSpinner, relevant only in case of Z/gamma* polarization weights (set "sensible" default)
0104       ipp_(true),  // pp collisions
0105       ipol_(0),
0106       nonSM2_(0),
0107       nonSMN_(0),
0108       cmsE_(
0109           13600),  // collision energy in GeV, relevant only in case of Z/gamma* polarization weights (set "sensible" default)
0110       default_weight_(config.getParameter<double>(
0111           "defaultWeight"))  // default weight stored in case of presence of a tau decay unsupported by TauSpinner
0112 {
0113   printModuleInfo(config);
0114
0115   // State that we use tauola/tauspinner resource
0116   usesResource(edm::SharedResourceNames::kTauola);
0117
0118   produces<nanoaod::FlatTable>();
0119 }
0120
0121 void TauSpinnerTableProducer::getBosons(edm::RefVector<edm::View<reco::GenParticle>> &bosons,
0122                                         const edm::View<reco::GenParticle> &parts) const {
0123   unsigned idx = 0;
0124   for (const auto &part : parts) {
0125     if (std::abs(part.pdgId()) == bosonPdgId_ && part.isLastCopy()) {
0126       edm::Ref<edm::View<reco::GenParticle>> partRef(&parts, idx);
0127       bosons.push_back(partRef);
0128     }
0129     ++idx;
0130   }
0131 }
0132
0133 reco::GenParticleRef TauSpinnerTableProducer::getLastCopy(const reco::GenParticleRef &part) {
0134   if (part->statusFlags().isLastCopy())
0135     return part;
0136   for (const auto &daughter : part->daughterRefVector()) {
0137     if (daughter->pdgId() == part->pdgId() && daughter->statusFlags().fromHardProcess()) {
0138       return getLastCopy(daughter);
0139     }
0140   }
0141   throw std::runtime_error("getLastCopy: no last copy found");
0142 }
0143
0144 void TauSpinnerTableProducer::getTaus(reco::GenParticleRefVector &taus, const reco::GenParticle &boson) {
0145   for (const auto &daughterRef : boson.daughterRefVector()) {
0146     if (std::abs(daughterRef->pdgId()) == 15)
0147       taus.push_back(getLastCopy(daughterRef));
0148   }
0149 }
0150
0151 bool TauSpinnerTableProducer::getTauDaughters(reco::GenParticleRefVector &tau_daughters, const reco::GenParticle &tau) {
0152   static const std::set<int> directTauProducts = {11, 12, 13, 14, 16, 22};
0153   static const std::set<int> finalHadrons = {111, 130, 211, 310, 311, 321};
0154   static const std::set<int> intermediateHadrons = {221, 223, 323};
0155   for (auto daughterRef : tau.daughterRefVector()) {
0156     const int daughter_pdgId = std::abs(daughterRef->pdgId());
0157     if ((std::abs(tau.pdgId()) == 15 && directTauProducts.count(daughter_pdgId)) ||
0158         finalHadrons.count(daughter_pdgId)) {
0159       tau_daughters.push_back(daughterRef);
0160     } else if (intermediateHadrons.count(daughter_pdgId)) {
0161       if (!getTauDaughters(tau_daughters, *daughterRef))
0162         return false;
0163     } else {
0164       edm::LogWarning("TauSpinnerTableProducer::getTauDaughters")
0165           << "Unsupported decay with " << daughter_pdgId << " being daughter of " << std::abs(tau.pdgId()) << "\n";
0166       return false;
0167     }
0168   }
0169   return true;
0170 }
0171
0172 void TauSpinnerTableProducer::beginJob() {
0173   // Initialize TauSpinner
0174   Tauolapp::Tauola::setNewCurrents(0);
0175   Tauolapp::Tauola::initialize();
0176   LHAPDF::initPDFSetByName(tauSpinnerPDF_);
0177   TauSpinner::initialize_spinner(ipp_, ipol_, nonSM2_, nonSMN_, cmsE_);
0178 }
0179
0180 void TauSpinnerTableProducer::produce(edm::Event &event, const edm::EventSetup &setup) {
0181   // Input gen-particles collection
0182   auto const &genParts = event.get(genPartsToken_);
0183
0184   // Output table
0185   auto wtTable = std::make_unique<nanoaod::FlatTable>(1, name_, true);
0186   wtTable->setDoc("TauSpinner weights");
0187
0188   // Search for boson
0189   edm::RefVector<edm::View<reco::GenParticle>> bosons;
0190   getBosons(bosons, genParts);
0191   if (bosons.size() !=
0192       1) {  // no boson found or more than one found, produce empty table (expected for non HTT samples)
0193     event.put(std::move(wtTable));
0194     return;
0195   }
0196
0197   // Search for taus from boson decay
0198   reco::GenParticleRefVector taus;
0199   getTaus(taus, *bosons[0]);
0200   if (taus.size() != 2) {  // boson does not decay to tau pair, produce empty table (expected for non HTT samples)
0201     event.put(std::move(wtTable));
0202     return;
0203   }
0204
0205   // Get tau daughters and convert all particles to TauSpinner format
0206   TauSpinner::SimpleParticle simple_boson = convertToSimplePart(*bosons[0]);
0207   std::array<TauSpinner::SimpleParticle, 2> simple_taus;
0208   std::array<std::vector<TauSpinner::SimpleParticle>, 2> simple_tau_daughters;
0209   bool supportedDecays = true;
0210   for (size_t tau_idx = 0; tau_idx < 2; ++tau_idx) {
0211     simple_taus[tau_idx] = convertToSimplePart(*taus[tau_idx]);
0212     reco::GenParticleRefVector tau_daughters;
0213     supportedDecays &= getTauDaughters(tau_daughters, *taus[tau_idx]);
0214     for (const auto &daughterRef : tau_daughters)
0215       simple_tau_daughters[tau_idx].push_back(convertToSimplePart(*daughterRef));
0216   }
0217
0218   // Compute TauSpinner weights and fill table
0219   std::array<double, 2> weights;
0220   for (const auto &theta : theta_vec_) {
0221     // Can make this more general by having boson pdgid as input or have option for set boson type
0222     TauSpinner::setHiggsParametersTR(-cos(2 * M_PI * theta.second),
0223                                      cos(2 * M_PI * theta.second),
0224                                      -sin(2 * M_PI * theta.second),
0225                                      -sin(2 * M_PI * theta.second));
0226     for (size_t i = 0; i < weights.size(); ++i) {
0227       Tauolapp::Tauola::setNewCurrents(i);
0228       weights[i] =
0229           supportedDecays
0230               ? TauSpinner::calculateWeightFromParticlesH(
0231                     simple_boson, simple_taus[0], simple_taus[1], simple_tau_daughters[0], simple_tau_daughters[1])
0232               : default_weight_;
0233     }
0234     // Nominal weights for setNewCurrents(0)
0235     wtTable->addColumnValue<double>(
0236         "weight_cp_" + theta.first, weights[0], "TauSpinner weight for theta_CP=" + theta.first);
0237     // Weights for alternative hadronic currents (can be used for uncertainty estimates)
0238     wtTable->addColumnValue<double>(
0239         "weight_cp_" + theta.first + "_alt",
0240         weights[1],
0241         "TauSpinner weight for theta_CP=" + theta.first + " (alternative hadronic currents)");
0242   }
0243
0244   event.put(std::move(wtTable));
0245 }
0246
0247 #include "FWCore/Framework/interface/MakerMacros.h"
0248 DEFINE_FWK_MODULE(TauSpinnerTableProducer);