plugins/TauSpinner/TauSpinnerCMS.cc

0001 #include "GeneratorInterface/TauolaInterface/interface/TauSpinnerCMS.h"
0002
0003 //MC-TESTER header files
0004 #include "Tauola/Tauola.h"
0005 #include "TauSpinner/tau_reweight_lib.h"
0006 #include "TauSpinner/Tauola_wrapper.h"
0007 #include "GeneratorInterface/TauolaInterface/interface/read_particles_from_HepMC.h"
0008 #include "TLorentzVector.h"
0009
0010 #include "CLHEP/Random/RandomEngine.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0015
0016 using namespace edm;
0017 using namespace TauSpinner;
0018
0019 CLHEP::HepRandomEngine *TauSpinnerCMS::fRandomEngine = nullptr;
0020 bool TauSpinnerCMS::isTauSpinnerConfigure = false;
0021 bool TauSpinnerCMS::fInitialized = false;
0022
0023 TauSpinnerCMS::TauSpinnerCMS(const ParameterSet &pset)
0024     : EDProducer(),
0025       isReco_(pset.getParameter<bool>("isReco")),
0026       isTauolaConfigured_(pset.getParameter<bool>("isTauolaConfigured")),
0027       isLHPDFConfigured_(pset.getParameter<bool>("isLHPDFConfigured")),
0028       LHAPDFname_(pset.getUntrackedParameter("LHAPDFname", (string)("MSTW2008nnlo90cl.LHgrid"))),
0029       CMSEnergy_(pset.getParameter<double>("CMSEnergy"))  //GeV
0030       ,
0031       gensrc_(pset.getParameter<edm::InputTag>("gensrc")),
0032       MotherPDGID_(pset.getUntrackedParameter("MotherPDGID", (int)(-1))),
0033       Ipol_(pset.getUntrackedParameter("Ipol", (int)(0))),
0034       nonSM2_(pset.getUntrackedParameter("nonSM2", (int)(0))),
0035       nonSMN_(pset.getUntrackedParameter("nonSMN", (int)(0))),
0036       roundOff_(pset.getUntrackedParameter("roundOff", (double)(0.01))) {
0037   //usesResource(edm::uniqueSharedResourceName());
0038   usesResource(edm::SharedResourceNames::kTauola);
0039
0040   produces<bool>("TauSpinnerWTisValid").setBranchAlias("TauSpinnerWTisValid");
0041   produces<double>("TauSpinnerWT").setBranchAlias("TauSpinnerWT");
0042   produces<double>("TauSpinnerWTFlip").setBranchAlias("TauSpinnerWTFlip");
0043   produces<double>("TauSpinnerWThplus").setBranchAlias("TauSpinnerWThplus");
0044   produces<double>("TauSpinnerWThminus").setBranchAlias("TauSpinnerWThminus");
0045
0046   if (isReco_) {
0047     GenParticleCollectionToken_ = consumes<reco::GenParticleCollection>(gensrc_);
0048   } else {
0049     hepmcCollectionToken_ = consumes<HepMCProduct>(gensrc_);
0050   }
0051 }
0052
0053 void TauSpinnerCMS::beginJob() {}
0054 void TauSpinnerCMS::endJob() {}
0055
0056 void TauSpinnerCMS::initialize() {
0057   // Now for Tauola and TauSpinner
0058   Tauolapp::Tauola::setRandomGenerator(TauSpinnerCMS::flat);
0059   if (!isTauolaConfigured_) {
0060     Tauolapp::Tauola::initialize();
0061   }
0062   if (!isLHPDFConfigured_) {
0063     LHAPDF::initPDFSetByName(LHAPDFname_);
0064   }
0065   if (!isTauSpinnerConfigure) {
0066     isTauSpinnerConfigure = true;
0067     bool Ipp = true;  // for pp collisions
0068     // Initialize TauSpinner
0069     //Ipol - polarization of input sample
0070     //nonSM2 - nonstandard model calculations
0071     //nonSMN
0072     TauSpinner::initialize_spinner(Ipp, Ipol_, nonSM2_, nonSMN_, CMSEnergy_);
0073   }
0074   fInitialized = true;
0075 }
0076
0077 void TauSpinnerCMS::produce(edm::Event &e, const edm::EventSetup &iSetup) {
0078   RandomEngineSentry<TauSpinnerCMS> randomEngineSentry(this, e.streamID());
0079   if (!fInitialized)
0080     initialize();
0081
0082   Tauolapp::Tauola::setRandomGenerator(
0083       TauSpinnerCMS::flat);  // rest tauola++ random number incase other modules use tauola++
0084   Tauolapp::jaki_.ktom = 1;  // rest for when you run after tauola
0085   double WT = 1.0;
0086   double WTFlip = 1.0;
0087   double polSM = -999;  //range [-1,1]
0088   SimpleParticle X, tau, tau2;
0089   std::vector<SimpleParticle> tau_daughters, tau_daughters2;
0090   int stat(0);
0091   if (isReco_) {
0092     stat = readParticlesfromReco(e, X, tau, tau2, tau_daughters, tau_daughters2);
0093   } else {
0094     edm::Handle<HepMCProduct> evt;
0095     e.getByToken(hepmcCollectionToken_, evt);
0096     //Get EVENT
0097     HepMC::GenEvent *Evt = new HepMC::GenEvent(*(evt->GetEvent()));
0098     stat = readParticlesFromHepMC(Evt, X, tau, tau2, tau_daughters, tau_daughters2);
0099   }
0100   if (MotherPDGID_ < 0 || abs(X.pdgid()) == MotherPDGID_) {
0101     if (stat != 1) {
0102       // Determine the weight
0103       if (abs(X.pdgid()) == 24 || abs(X.pdgid()) == 37) {
0104         TLorentzVector tau_1r(0, 0, 0, 0);
0105         TLorentzVector tau_1(tau.px(), tau.py(), tau.pz(), tau.e());
0106         for (unsigned int i = 0; i < tau_daughters.size(); i++) {
0107           tau_1r += TLorentzVector(
0108               tau_daughters.at(i).px(), tau_daughters.at(i).py(), tau_daughters.at(i).pz(), tau_daughters.at(i).e());
0109         }
0110         if (fabs(tau_1r.M() - tau_1.M()) < roundOff_) {
0111           WT = TauSpinner::calculateWeightFromParticlesWorHpn(
0112               X, tau, tau2, tau_daughters);  // note that tau2 is tau neutrino
0113           polSM = getTauSpin();
0114           WTFlip = (2.0 - WT) / WT;
0115         }
0116       } else if (X.pdgid() == 25 || X.pdgid() == 36 || X.pdgid() == 22 || X.pdgid() == 23) {
0117         TLorentzVector tau_1r(0, 0, 0, 0), tau_2r(0, 0, 0, 0);
0118         TLorentzVector tau_1(tau.px(), tau.py(), tau.pz(), tau.e()), tau_2(tau2.px(), tau2.py(), tau2.pz(), tau2.e());
0119         for (unsigned int i = 0; i < tau_daughters.size(); i++) {
0120           tau_1r += TLorentzVector(
0121               tau_daughters.at(i).px(), tau_daughters.at(i).py(), tau_daughters.at(i).pz(), tau_daughters.at(i).e());
0122         }
0123         for (unsigned int i = 0; i < tau_daughters2.size(); i++) {
0124           tau_2r += TLorentzVector(tau_daughters2.at(i).px(),
0125                                    tau_daughters2.at(i).py(),
0126                                    tau_daughters2.at(i).pz(),
0127                                    tau_daughters2.at(i).e());
0128         }
0129
0130         if (fabs(tau_1r.M() - tau_1.M()) < roundOff_ && fabs(tau_2r.M() - tau_2.M()) < roundOff_) {
0131           WT = TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters, tau_daughters2);
0132           //std::cout << "WT " << WT << std::endl;
0133           polSM = getTauSpin();
0134           if (X.pdgid() == 25 || X.pdgid() == 22 || X.pdgid() == 23) {
0135             if (X.pdgid() == 25)
0136               X.setPdgid(23);
0137             if (X.pdgid() == 22 || X.pdgid() == 23)
0138               X.setPdgid(25);
0139
0140             double WTother = TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters, tau_daughters2);
0141             WTFlip = WTother / WT;
0142           }
0143         }
0144       } else {
0145         cout << "TauSpinner: WARNING: Unexpected PDG for tau mother: " << X.pdgid() << endl;
0146       }
0147     }
0148   }
0149   bool isValid = true;
0150   if (!(0 <= WT && WT < 10)) {
0151     isValid = false;
0152     WT = 1.0;
0153     WTFlip = 1.0;
0154   }
0155   std::unique_ptr<bool> TauSpinnerWeightisValid(new bool);
0156   *TauSpinnerWeightisValid = isValid;
0157   e.put(std::move(TauSpinnerWeightisValid), "TauSpinnerWTisValid");
0158
0159   // regular weight
0160   std::unique_ptr<double> TauSpinnerWeight(new double);
0161   *TauSpinnerWeight = WT;
0162   e.put(std::move(TauSpinnerWeight), "TauSpinnerWT");
0163
0164   // flipped weight (ie Z->H or H->Z)
0165   std::unique_ptr<double> TauSpinnerWeightFlip(new double);
0166   *TauSpinnerWeightFlip = WTFlip;
0167   e.put(std::move(TauSpinnerWeightFlip), "TauSpinnerWTFlip");
0168
0169   // h+ polarization
0170   double WThplus = WT;
0171   if (polSM < 0.0 && polSM != -999 && isValid)
0172     WThplus = 0;
0173   std::unique_ptr<double> TauSpinnerWeighthplus(new double);
0174   *TauSpinnerWeighthplus = WThplus;
0175   e.put(std::move(TauSpinnerWeighthplus), "TauSpinnerWThplus");
0176
0177   // h- polarization
0178   double WThminus = WT;
0179   if (polSM > 0.0 && polSM != -999 && isValid)
0180     WThminus = 0;
0181   std::unique_ptr<double> TauSpinnerWeighthminus(new double);
0182   *TauSpinnerWeighthminus = WThminus;
0183   e.put(std::move(TauSpinnerWeighthminus), "TauSpinnerWThminus");
0184   return;
0185 }
0186
0187 int TauSpinnerCMS::readParticlesfromReco(edm::Event &e,
0188                                          SimpleParticle &X,
0189                                          SimpleParticle &tau,
0190                                          SimpleParticle &tau2,
0191                                          std::vector<SimpleParticle> &tau_daughters,
0192                                          std::vector<SimpleParticle> &tau2_daughters) {
0193   edm::Handle<reco::GenParticleCollection> genParticles;
0194   e.getByToken(GenParticleCollectionToken_, genParticles);
0195   for (reco::GenParticleCollection::const_iterator itr = genParticles->begin(); itr != genParticles->end(); ++itr) {
0196     int pdgid = abs(itr->pdgId());
0197     if (pdgid == 24 || pdgid == 37 || pdgid == 25 || pdgid == 36 || pdgid == 22 || pdgid == 23) {
0198       const reco::GenParticle *hx = &(*itr);
0199       if (!isFirst(hx))
0200         continue;
0201       GetLastSelf(hx);
0202       const reco::GenParticle *recotau1 = nullptr;
0203       const reco::GenParticle *recotau2 = nullptr;
0204       unsigned int ntau(0), ntauornu(0);
0205       for (unsigned int i = 0; i < itr->numberOfDaughters(); i++) {
0206         const reco::Candidate *dau = itr->daughter(i);
0207         if (abs(dau->pdgId()) != pdgid) {
0208           if (abs(dau->pdgId()) == 15 || abs(dau->pdgId()) == 16) {
0209             if (ntau == 0 && abs(dau->pdgId()) == 15) {
0210               recotau1 = static_cast<const reco::GenParticle *>(dau);
0211               GetLastSelf(recotau1);
0212               ntau++;
0213             } else if ((ntau == 1 && abs(dau->pdgId()) == 15) || abs(dau->pdgId()) == 16) {
0214               recotau2 = static_cast<const reco::GenParticle *>(dau);
0215               if (abs(dau->pdgId()) == 15) {
0216                 ntau++;
0217                 GetLastSelf(recotau2);
0218               }
0219             }
0220             ntauornu++;
0221           }
0222         }
0223       }
0224       if ((ntau == 2 && ntauornu == 2) || (ntau == 1 && ntauornu == 2)) {
0225         X.setPx(itr->p4().Px());
0226         X.setPy(itr->p4().Py());
0227         X.setPz(itr->p4().Pz());
0228         X.setE(itr->p4().E());
0229         X.setPdgid(itr->pdgId());
0230         tau.setPx(recotau1->p4().Px());
0231         tau.setPy(recotau1->p4().Py());
0232         tau.setPz(recotau1->p4().Pz());
0233         tau.setE(recotau1->p4().E());
0234         tau.setPdgid(recotau1->pdgId());
0235         GetRecoDaughters(recotau1, tau_daughters, recotau1->pdgId());
0236         tau2.setPx(recotau2->p4().Px());
0237         tau2.setPy(recotau2->p4().Py());
0238         tau2.setPz(recotau2->p4().Pz());
0239         tau2.setE(recotau2->p4().E());
0240         tau2.setPdgid(recotau2->pdgId());
0241         if (ntau == 2)
0242           GetRecoDaughters(recotau2, tau2_daughters, recotau2->pdgId());
0243         return 0;
0244       }
0245     }
0246   }
0247   return 1;
0248 }
0249
0250 void TauSpinnerCMS::GetLastSelf(const reco::GenParticle *Particle) {
0251   for (unsigned int i = 0; i < Particle->numberOfDaughters(); i++) {
0252     const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(Particle->daughter(i));
0253     if (Particle->pdgId() == dau->pdgId()) {
0254       Particle = dau;
0255       GetLastSelf(Particle);
0256     }
0257   }
0258 }
0259
0260 bool TauSpinnerCMS::isFirst(const reco::GenParticle *Particle) {
0261   for (unsigned int i = 0; i < Particle->numberOfMothers(); i++) {
0262     const reco::GenParticle *moth = static_cast<const reco::GenParticle *>(Particle->mother(i));
0263     if (Particle->pdgId() == moth->pdgId()) {
0264       return false;
0265     }
0266   }
0267   return true;
0268 }
0269
0270 void TauSpinnerCMS::GetRecoDaughters(const reco::GenParticle *Particle,
0271                                      std::vector<SimpleParticle> &daughters,
0272                                      int parentpdgid) {
0273   if (Particle->numberOfDaughters() == 0 || abs(Particle->pdgId()) == 111) {
0274     SimpleParticle tp(
0275         Particle->p4().Px(), Particle->p4().Py(), Particle->p4().Pz(), Particle->p4().E(), Particle->pdgId());
0276     daughters.push_back(tp);
0277     return;
0278   }
0279   for (unsigned int i = 0; i < Particle->numberOfDaughters(); i++) {
0280     const reco::Candidate *dau = Particle->daughter(i);
0281     GetRecoDaughters(static_cast<const reco::GenParticle *>(dau), daughters, Particle->pdgId());
0282   }
0283 }
0284
0285 double TauSpinnerCMS::flat() {
0286   if (!fRandomEngine) {
0287     throw cms::Exception("LogicError")
0288         << "TauSpinnerCMS::flat: Attempt to generate random number when engine pointer is null\n"
0289         << "This might mean that the code was modified to generate a random number outside the\n"
0290         << "event and beginLuminosityBlock methods, which is not allowed.\n";
0291   }
0292   return fRandomEngine->flat();
0293 }
0294
0295 DEFINE_FWK_MODULE(TauSpinnerCMS);