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