Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:22 

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 
0005 #include "PhysicsTools/JetMCUtils/interface/combination.h"
0006 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
0007 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
0008 
0009 template <typename LeptonCollection>
0010 class TtSemiLepKinFitProducer : public edm::stream::EDProducer<> {
0011 public:
0012   explicit TtSemiLepKinFitProducer(const edm::ParameterSet&);
0013 
0014 private:
0015   // produce
0016   void produce(edm::Event&, const edm::EventSetup&) override;
0017 
0018   // convert unsigned to Param
0019   TtSemiLepKinFitter::Param param(unsigned);
0020   // convert unsigned to Param
0021   TtSemiLepKinFitter::Constraint constraint(unsigned);
0022   // convert unsigned to Param
0023   std::vector<TtSemiLepKinFitter::Constraint> constraints(std::vector<unsigned>&);
0024   // helper function for b-tagging
0025   bool doBTagging(bool& useBTag_,
0026                   const edm::Handle<std::vector<pat::Jet>>& jets,
0027                   const std::vector<int>& combi,
0028                   std::string& bTagAlgo_,
0029                   double& minBTagValueBJets_,
0030                   double& maxBTagValueNonBJets_);
0031 
0032   edm::EDGetTokenT<std::vector<pat::Jet>> jetsToken_;
0033   edm::EDGetTokenT<LeptonCollection> lepsToken_;
0034   edm::EDGetTokenT<std::vector<pat::MET>> metsToken_;
0035 
0036   edm::EDGetTokenT<std::vector<std::vector<int>>> matchToken_;
0037   /// switch to use only a combination given by another hypothesis
0038   bool useOnlyMatch_;
0039   /// input tag for b-tagging algorithm
0040   std::string bTagAlgo_;
0041   /// min value of bTag for a b-jet
0042   double minBTagValueBJet_;
0043   /// max value of bTag for a non-b-jet
0044   double maxBTagValueNonBJet_;
0045   /// switch to tell whether to use b-tagging or not
0046   bool useBTag_;
0047   /// maximal number of jets (-1 possible to indicate 'all')
0048   int maxNJets_;
0049   /// maximal number of combinations to be written to the event
0050   int maxNComb_;
0051 
0052   /// maximal number of iterations to be performed for the fit
0053   unsigned int maxNrIter_;
0054   /// maximal chi2 equivalent
0055   double maxDeltaS_;
0056   /// maximal deviation for contstraints
0057   double maxF_;
0058   unsigned int jetParam_;
0059   unsigned int lepParam_;
0060   unsigned int metParam_;
0061   /// constrains
0062   std::vector<unsigned> constraints_;
0063   double mW_;
0064   double mTop_;
0065   /// scale factors for jet energy resolution
0066   std::vector<double> jetEnergyResolutionScaleFactors_;
0067   std::vector<double> jetEnergyResolutionEtaBinning_;
0068   /// config-file-based object resolutions
0069   std::vector<edm::ParameterSet> udscResolutions_;
0070   std::vector<edm::ParameterSet> bResolutions_;
0071   std::vector<edm::ParameterSet> lepResolutions_;
0072   std::vector<edm::ParameterSet> metResolutions_;
0073 
0074   std::unique_ptr<TtSemiLepKinFitter> fitter;
0075 
0076   struct KinFitResult {
0077     int Status;
0078     double Chi2;
0079     double Prob;
0080     pat::Particle HadB;
0081     pat::Particle HadP;
0082     pat::Particle HadQ;
0083     pat::Particle LepB;
0084     pat::Particle LepL;
0085     pat::Particle LepN;
0086     std::vector<int> JetCombi;
0087     bool operator<(const KinFitResult& rhs) { return Chi2 < rhs.Chi2; };
0088   };
0089 };
0090 
0091 template <typename LeptonCollection>
0092 TtSemiLepKinFitProducer<LeptonCollection>::TtSemiLepKinFitProducer(const edm::ParameterSet& cfg)
0093     : jetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("jets"))),
0094       lepsToken_(consumes<LeptonCollection>(cfg.getParameter<edm::InputTag>("leps"))),
0095       metsToken_(consumes<std::vector<pat::MET>>(cfg.getParameter<edm::InputTag>("mets"))),
0096       matchToken_(mayConsume<std::vector<std::vector<int>>>(cfg.getParameter<edm::InputTag>("match"))),
0097       useOnlyMatch_(cfg.getParameter<bool>("useOnlyMatch")),
0098       bTagAlgo_(cfg.getParameter<std::string>("bTagAlgo")),
0099       minBTagValueBJet_(cfg.getParameter<double>("minBDiscBJets")),
0100       maxBTagValueNonBJet_(cfg.getParameter<double>("maxBDiscLightJets")),
0101       useBTag_(cfg.getParameter<bool>("useBTagging")),
0102       maxNJets_(cfg.getParameter<int>("maxNJets")),
0103       maxNComb_(cfg.getParameter<int>("maxNComb")),
0104       maxNrIter_(cfg.getParameter<unsigned>("maxNrIter")),
0105       maxDeltaS_(cfg.getParameter<double>("maxDeltaS")),
0106       maxF_(cfg.getParameter<double>("maxF")),
0107       jetParam_(cfg.getParameter<unsigned>("jetParametrisation")),
0108       lepParam_(cfg.getParameter<unsigned>("lepParametrisation")),
0109       metParam_(cfg.getParameter<unsigned>("metParametrisation")),
0110       constraints_(cfg.getParameter<std::vector<unsigned>>("constraints")),
0111       mW_(cfg.getParameter<double>("mW")),
0112       mTop_(cfg.getParameter<double>("mTop")),
0113       jetEnergyResolutionScaleFactors_(cfg.getParameter<std::vector<double>>("jetEnergyResolutionScaleFactors")),
0114       jetEnergyResolutionEtaBinning_(cfg.getParameter<std::vector<double>>("jetEnergyResolutionEtaBinning")),
0115       udscResolutions_(0),
0116       bResolutions_(0),
0117       lepResolutions_(0),
0118       metResolutions_(0) {
0119   if (cfg.exists("udscResolutions") && cfg.exists("bResolutions") && cfg.exists("lepResolutions") &&
0120       cfg.exists("metResolutions")) {
0121     udscResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet>>("udscResolutions");
0122     bResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet>>("bResolutions");
0123     lepResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet>>("lepResolutions");
0124     metResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet>>("metResolutions");
0125   } else if (cfg.exists("udscResolutions") || cfg.exists("bResolutions") || cfg.exists("lepResolutions") ||
0126              cfg.exists("metResolutions")) {
0127     throw cms::Exception("Configuration") << "Parameters 'udscResolutions', 'bResolutions', 'lepResolutions', "
0128                                              "'metResolutions' should be used together.\n";
0129   }
0130 
0131   fitter = std::make_unique<TtSemiLepKinFitter>(param(jetParam_),
0132                                                 param(lepParam_),
0133                                                 param(metParam_),
0134                                                 maxNrIter_,
0135                                                 maxDeltaS_,
0136                                                 maxF_,
0137                                                 constraints(constraints_),
0138                                                 mW_,
0139                                                 mTop_,
0140                                                 &udscResolutions_,
0141                                                 &bResolutions_,
0142                                                 &lepResolutions_,
0143                                                 &metResolutions_,
0144                                                 &jetEnergyResolutionScaleFactors_,
0145                                                 &jetEnergyResolutionEtaBinning_);
0146 
0147   produces<std::vector<pat::Particle>>("PartonsHadP");
0148   produces<std::vector<pat::Particle>>("PartonsHadQ");
0149   produces<std::vector<pat::Particle>>("PartonsHadB");
0150   produces<std::vector<pat::Particle>>("PartonsLepB");
0151   produces<std::vector<pat::Particle>>("Leptons");
0152   produces<std::vector<pat::Particle>>("Neutrinos");
0153 
0154   produces<std::vector<std::vector<int>>>();
0155   produces<std::vector<double>>("Chi2");
0156   produces<std::vector<double>>("Prob");
0157   produces<std::vector<int>>("Status");
0158 
0159   produces<int>("NumberOfConsideredJets");
0160 }
0161 
0162 template <typename LeptonCollection>
0163 bool TtSemiLepKinFitProducer<LeptonCollection>::doBTagging(bool& useBTag_,
0164                                                            const edm::Handle<std::vector<pat::Jet>>& jets,
0165                                                            const std::vector<int>& combi,
0166                                                            std::string& bTagAlgo_,
0167                                                            double& minBTagValueBJet_,
0168                                                            double& maxBTagValueNonBJet_) {
0169   if (!useBTag_) {
0170     return true;
0171   }
0172   if (useBTag_ && (*jets)[combi[TtSemiLepEvtPartons::HadB]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0173       (*jets)[combi[TtSemiLepEvtPartons::LepB]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0174       (*jets)[combi[TtSemiLepEvtPartons::LightQ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0175       (*jets)[combi[TtSemiLepEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_) {
0176     return true;
0177   } else {
0178     return false;
0179   }
0180 }
0181 
0182 template <typename LeptonCollection>
0183 void TtSemiLepKinFitProducer<LeptonCollection>::produce(edm::Event& evt, const edm::EventSetup& setup) {
0184   std::unique_ptr<std::vector<pat::Particle>> pPartonsHadP(new std::vector<pat::Particle>);
0185   std::unique_ptr<std::vector<pat::Particle>> pPartonsHadQ(new std::vector<pat::Particle>);
0186   std::unique_ptr<std::vector<pat::Particle>> pPartonsHadB(new std::vector<pat::Particle>);
0187   std::unique_ptr<std::vector<pat::Particle>> pPartonsLepB(new std::vector<pat::Particle>);
0188   std::unique_ptr<std::vector<pat::Particle>> pLeptons(new std::vector<pat::Particle>);
0189   std::unique_ptr<std::vector<pat::Particle>> pNeutrinos(new std::vector<pat::Particle>);
0190 
0191   std::unique_ptr<std::vector<std::vector<int>>> pCombi(new std::vector<std::vector<int>>);
0192   std::unique_ptr<std::vector<double>> pChi2(new std::vector<double>);
0193   std::unique_ptr<std::vector<double>> pProb(new std::vector<double>);
0194   std::unique_ptr<std::vector<int>> pStatus(new std::vector<int>);
0195 
0196   std::unique_ptr<int> pJetsConsidered(new int);
0197 
0198   const edm::Handle<std::vector<pat::Jet>>& jets = evt.getHandle(jetsToken_);
0199 
0200   const edm::Handle<std::vector<pat::MET>>& mets = evt.getHandle(metsToken_);
0201 
0202   const edm::Handle<LeptonCollection>& leps = evt.getHandle(lepsToken_);
0203 
0204   const unsigned int nPartons = 4;
0205 
0206   std::vector<int> match;
0207   bool invalidMatch = false;
0208   if (useOnlyMatch_) {
0209     *pJetsConsidered = nPartons;
0210     const edm::Handle<std::vector<std::vector<int>>>& matchHandle = evt.getHandle(matchToken_);
0211     match = *(matchHandle->begin());
0212     // check if match is valid
0213     if (match.size() != nPartons)
0214       invalidMatch = true;
0215     else {
0216       for (unsigned int idx = 0; idx < match.size(); ++idx) {
0217         if (match[idx] < 0 || match[idx] >= (int)jets->size()) {
0218           invalidMatch = true;
0219           break;
0220         }
0221       }
0222     }
0223   }
0224 
0225   // -----------------------------------------------------
0226   // skip events with no appropriate lepton candidate in
0227   // or empty MET or less jets than partons or invalid match
0228   // -----------------------------------------------------
0229 
0230   if (leps->empty() || mets->empty() || jets->size() < nPartons || invalidMatch) {
0231     // the kinFit getters return empty objects here
0232     pPartonsHadP->push_back(fitter->fittedHadP());
0233     pPartonsHadQ->push_back(fitter->fittedHadQ());
0234     pPartonsHadB->push_back(fitter->fittedHadB());
0235     pPartonsLepB->push_back(fitter->fittedLepB());
0236     pLeptons->push_back(fitter->fittedLepton());
0237     pNeutrinos->push_back(fitter->fittedNeutrino());
0238     // indices referring to the jet combination
0239     std::vector<int> invalidCombi;
0240     for (unsigned int i = 0; i < nPartons; ++i)
0241       invalidCombi.push_back(-1);
0242     pCombi->push_back(invalidCombi);
0243     // chi2
0244     pChi2->push_back(-1.);
0245     // chi2 probability
0246     pProb->push_back(-1.);
0247     // status of the fitter
0248     pStatus->push_back(-1);
0249     // number of jets
0250     *pJetsConsidered = jets->size();
0251     // feed out all products
0252     evt.put(std::move(pCombi));
0253     evt.put(std::move(pPartonsHadP), "PartonsHadP");
0254     evt.put(std::move(pPartonsHadQ), "PartonsHadQ");
0255     evt.put(std::move(pPartonsHadB), "PartonsHadB");
0256     evt.put(std::move(pPartonsLepB), "PartonsLepB");
0257     evt.put(std::move(pLeptons), "Leptons");
0258     evt.put(std::move(pNeutrinos), "Neutrinos");
0259     evt.put(std::move(pChi2), "Chi2");
0260     evt.put(std::move(pProb), "Prob");
0261     evt.put(std::move(pStatus), "Status");
0262     evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0263     return;
0264   }
0265 
0266   // -----------------------------------------------------
0267   // analyze different jet combinations using the KinFitter
0268   // (or only a given jet combination if useOnlyMatch=true)
0269   // -----------------------------------------------------
0270 
0271   std::vector<int> jetIndices;
0272   if (!useOnlyMatch_) {
0273     for (unsigned int i = 0; i < jets->size(); ++i) {
0274       if (maxNJets_ >= (int)nPartons && maxNJets_ == (int)i) {
0275         *pJetsConsidered = i;
0276         break;
0277       }
0278       jetIndices.push_back(i);
0279     }
0280   }
0281 
0282   std::vector<int> combi;
0283   for (unsigned int i = 0; i < nPartons; ++i) {
0284     if (useOnlyMatch_)
0285       combi.push_back(match[i]);
0286     else
0287       combi.push_back(i);
0288   }
0289 
0290   std::list<KinFitResult> FitResultList;
0291 
0292   do {
0293     for (int cnt = 0; cnt < TMath::Factorial(combi.size()); ++cnt) {
0294       // take into account indistinguishability of the two jets from the hadr. W decay,
0295       // reduces combinatorics by a factor of 2
0296       if ((combi[TtSemiLepEvtPartons::LightQ] < combi[TtSemiLepEvtPartons::LightQBar] || useOnlyMatch_) &&
0297           doBTagging(useBTag_, jets, combi, bTagAlgo_, minBTagValueBJet_, maxBTagValueNonBJet_)) {
0298         std::vector<pat::Jet> jetCombi;
0299         jetCombi.resize(nPartons);
0300         jetCombi[TtSemiLepEvtPartons::LightQ] = (*jets)[combi[TtSemiLepEvtPartons::LightQ]];
0301         jetCombi[TtSemiLepEvtPartons::LightQBar] = (*jets)[combi[TtSemiLepEvtPartons::LightQBar]];
0302         jetCombi[TtSemiLepEvtPartons::HadB] = (*jets)[combi[TtSemiLepEvtPartons::HadB]];
0303         jetCombi[TtSemiLepEvtPartons::LepB] = (*jets)[combi[TtSemiLepEvtPartons::LepB]];
0304 
0305         // do the kinematic fit
0306         const int status = fitter->fit(jetCombi, (*leps)[0], (*mets)[0]);
0307 
0308         if (status == 0) {  // only take into account converged fits
0309           KinFitResult result;
0310           result.Status = status;
0311           result.Chi2 = fitter->fitS();
0312           result.Prob = fitter->fitProb();
0313           result.HadB = fitter->fittedHadB();
0314           result.HadP = fitter->fittedHadP();
0315           result.HadQ = fitter->fittedHadQ();
0316           result.LepB = fitter->fittedLepB();
0317           result.LepL = fitter->fittedLepton();
0318           result.LepN = fitter->fittedNeutrino();
0319           result.JetCombi = combi;
0320 
0321           FitResultList.push_back(result);
0322         }
0323       }
0324       if (useOnlyMatch_)
0325         break;  // don't go through combinatorics if useOnlyMatch was chosen
0326       next_permutation(combi.begin(), combi.end());
0327     }
0328     if (useOnlyMatch_)
0329       break;  // don't go through combinatorics if useOnlyMatch was chosen
0330   } while (stdcomb::next_combination(jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end()));
0331 
0332   // sort results w.r.t. chi2 values
0333   FitResultList.sort();
0334 
0335   // -----------------------------------------------------
0336   // feed out result
0337   // starting with the JetComb having the smallest chi2
0338   // -----------------------------------------------------
0339 
0340   if ((unsigned)FitResultList.size() < 1) {  // in case no fit results were stored in the list (all fits aborted)
0341     pPartonsHadP->push_back(fitter->fittedHadP());
0342     pPartonsHadQ->push_back(fitter->fittedHadQ());
0343     pPartonsHadB->push_back(fitter->fittedHadB());
0344     pPartonsLepB->push_back(fitter->fittedLepB());
0345     pLeptons->push_back(fitter->fittedLepton());
0346     pNeutrinos->push_back(fitter->fittedNeutrino());
0347     // indices referring to the jet combination
0348     std::vector<int> invalidCombi;
0349     for (unsigned int i = 0; i < nPartons; ++i)
0350       invalidCombi.push_back(-1);
0351     pCombi->push_back(invalidCombi);
0352     // chi2
0353     pChi2->push_back(-1.);
0354     // chi2 probability
0355     pProb->push_back(-1.);
0356     // status of the fitter
0357     pStatus->push_back(-1);
0358   } else {
0359     unsigned int iComb = 0;
0360     for (typename std::list<KinFitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end();
0361          ++result) {
0362       if (maxNComb_ >= 1 && iComb == (unsigned int)maxNComb_)
0363         break;
0364       iComb++;
0365       // partons
0366       pPartonsHadP->push_back(result->HadP);
0367       pPartonsHadQ->push_back(result->HadQ);
0368       pPartonsHadB->push_back(result->HadB);
0369       pPartonsLepB->push_back(result->LepB);
0370       // lepton
0371       pLeptons->push_back(result->LepL);
0372       // neutrino
0373       pNeutrinos->push_back(result->LepN);
0374       // indices referring to the jet combination
0375       pCombi->push_back(result->JetCombi);
0376       // chi2
0377       pChi2->push_back(result->Chi2);
0378       // chi2 probability
0379       pProb->push_back(result->Prob);
0380       // status of the fitter
0381       pStatus->push_back(result->Status);
0382     }
0383   }
0384   evt.put(std::move(pCombi));
0385   evt.put(std::move(pPartonsHadP), "PartonsHadP");
0386   evt.put(std::move(pPartonsHadQ), "PartonsHadQ");
0387   evt.put(std::move(pPartonsHadB), "PartonsHadB");
0388   evt.put(std::move(pPartonsLepB), "PartonsLepB");
0389   evt.put(std::move(pLeptons), "Leptons");
0390   evt.put(std::move(pNeutrinos), "Neutrinos");
0391   evt.put(std::move(pChi2), "Chi2");
0392   evt.put(std::move(pProb), "Prob");
0393   evt.put(std::move(pStatus), "Status");
0394   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0395 }
0396 
0397 template <typename LeptonCollection>
0398 TtSemiLepKinFitter::Param TtSemiLepKinFitProducer<LeptonCollection>::param(unsigned val) {
0399   TtSemiLepKinFitter::Param result;
0400   switch (val) {
0401     case TtSemiLepKinFitter::kEMom:
0402       result = TtSemiLepKinFitter::kEMom;
0403       break;
0404     case TtSemiLepKinFitter::kEtEtaPhi:
0405       result = TtSemiLepKinFitter::kEtEtaPhi;
0406       break;
0407     case TtSemiLepKinFitter::kEtThetaPhi:
0408       result = TtSemiLepKinFitter::kEtThetaPhi;
0409       break;
0410     default:
0411       throw cms::Exception("Configuration") << "Chosen jet parametrization is not supported: " << val << "\n";
0412       break;
0413   }
0414   return result;
0415 }
0416 
0417 template <typename LeptonCollection>
0418 TtSemiLepKinFitter::Constraint TtSemiLepKinFitProducer<LeptonCollection>::constraint(unsigned val) {
0419   TtSemiLepKinFitter::Constraint result;
0420   switch (val) {
0421     case TtSemiLepKinFitter::kWHadMass:
0422       result = TtSemiLepKinFitter::kWHadMass;
0423       break;
0424     case TtSemiLepKinFitter::kWLepMass:
0425       result = TtSemiLepKinFitter::kWLepMass;
0426       break;
0427     case TtSemiLepKinFitter::kTopHadMass:
0428       result = TtSemiLepKinFitter::kTopHadMass;
0429       break;
0430     case TtSemiLepKinFitter::kTopLepMass:
0431       result = TtSemiLepKinFitter::kTopLepMass;
0432       break;
0433     case TtSemiLepKinFitter::kNeutrinoMass:
0434       result = TtSemiLepKinFitter::kNeutrinoMass;
0435       break;
0436     case TtSemiLepKinFitter::kEqualTopMasses:
0437       result = TtSemiLepKinFitter::kEqualTopMasses;
0438       break;
0439     case TtSemiLepKinFitter::kSumPt:
0440       result = TtSemiLepKinFitter::kSumPt;
0441       break;
0442     default:
0443       throw cms::Exception("Configuration") << "Chosen fit constraint is not supported: " << val << "\n";
0444       break;
0445   }
0446   return result;
0447 }
0448 
0449 template <typename LeptonCollection>
0450 std::vector<TtSemiLepKinFitter::Constraint> TtSemiLepKinFitProducer<LeptonCollection>::constraints(
0451     std::vector<unsigned>& val) {
0452   std::vector<TtSemiLepKinFitter::Constraint> result;
0453   for (unsigned i = 0; i < val.size(); ++i) {
0454     result.push_back(constraint(val[i]));
0455   }
0456   return result;
0457 }
0458 
0459 #include "DataFormats/PatCandidates/interface/Muon.h"
0460 #include "DataFormats/PatCandidates/interface/Electron.h"
0461 using TtSemiLepKinFitProducerMuon = TtSemiLepKinFitProducer<std::vector<pat::Muon>>;
0462 using TtSemiLepKinFitProducerElectron = TtSemiLepKinFitProducer<std::vector<pat::Electron>>;
0463 
0464 #include "FWCore/Framework/interface/MakerMacros.h"
0465 DEFINE_FWK_MODULE(TtSemiLepKinFitProducerMuon);
0466 DEFINE_FWK_MODULE(TtSemiLepKinFitProducerElectron);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!