TopKinFitter/plugins/TtFullHadKinFitProducer.cc

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 "TopQuarkAnalysis/TopKinFitter/interface/TtFullHadKinFitter.h"
0006 #include "AnalysisDataFormats/TopObjects/interface/TtFullHadEvtPartons.h"
0007
0008 /*
0009   \class   TtFullHadKinFitProducer TtFullHadKinFitProducer.h "TopQuarkAnalysis/TopKinFitter/plugins/TtFullHadKinFitProducer.h"
0010
0011   \brief   Retrieve kinFit result from TtFullHadKinFitter and put it into the event
0012
0013   Get jet collection and if wanted match from the event content and do the kinematic fit
0014   of the event with this objects using the kinFit class from TtFullHadKinFitter and put
0015   the result into the event content
0016
0017 **/
0018
0019 class TtFullHadKinFitProducer : public edm::stream::EDProducer<> {
0020 public:
0021   /// default constructor
0022   explicit TtFullHadKinFitProducer(const edm::ParameterSet& cfg);
0023
0024 private:
0025   /// produce fitted object collections and meta data describing fit quality
0026   void produce(edm::Event& event, const edm::EventSetup& setup) override;
0027
0028 private:
0029   /// input tag for jets
0030   edm::EDGetTokenT<std::vector<pat::Jet> > jetsToken_;
0031   /// input tag for matches (in case the fit should be performed on certain matches)
0032   edm::EDGetTokenT<std::vector<std::vector<int> > > matchToken_;
0033   /// switch to tell whether all possible combinations should be used for the fit
0034   /// or only a certain combination
0035   bool useOnlyMatch_;
0036   /// input tag for b-tagging algorithm
0037   std::string bTagAlgo_;
0038   /// min value of bTag for a b-jet
0039   double minBTagValueBJet_;
0040   /// max value of bTag for a non-b-jet
0041   double maxBTagValueNonBJet_;
0042   /// switch to tell whether to use b-tagging or not
0043   bool useBTagging_;
0044   /// minimal number of b-jets
0045   unsigned int bTags_;
0046   /// correction level for jets
0047   std::string jetCorrectionLevel_;
0048   /// maximal number of jets (-1 possible to indicate 'all')
0049   int maxNJets_;
0050   /// maximal number of combinations to be written to the event
0051   int maxNComb_;
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   /// numbering of different possible jet parametrizations
0059   unsigned int jetParam_;
0060   /// numbering of different possible kinematic constraints
0061   std::vector<unsigned> constraints_;
0062   /// W mass value used for constraints
0063   double mW_;
0064   /// top mass value used for constraints
0065   double mTop_;
0066   /// store the resolutions for the jets
0067   std::vector<edm::ParameterSet> udscResolutions_, bResolutions_;
0068   /// scale factors for jet energy resolution
0069   std::vector<double> jetEnergyResolutionScaleFactors_;
0070   std::vector<double> jetEnergyResolutionEtaBinning_;
0071
0072 public:
0073   /// kinematic fit interface
0074   std::unique_ptr<TtFullHadKinFitter::KinFit> kinFitter;
0075 };
0076
0077 static const unsigned int nPartons = 6;
0078
0079 /// default constructor
0080 TtFullHadKinFitProducer::TtFullHadKinFitProducer(const edm::ParameterSet& cfg)
0081     : jetsToken_(consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"))),
0082       matchToken_(mayConsume<std::vector<std::vector<int> > >(cfg.getParameter<edm::InputTag>("match"))),
0083       useOnlyMatch_(cfg.getParameter<bool>("useOnlyMatch")),
0084       bTagAlgo_(cfg.getParameter<std::string>("bTagAlgo")),
0085       minBTagValueBJet_(cfg.getParameter<double>("minBTagValueBJet")),
0086       maxBTagValueNonBJet_(cfg.getParameter<double>("maxBTagValueNonBJet")),
0087       useBTagging_(cfg.getParameter<bool>("useBTagging")),
0088       bTags_(cfg.getParameter<unsigned int>("bTags")),
0089       jetCorrectionLevel_(cfg.getParameter<std::string>("jetCorrectionLevel")),
0090       maxNJets_(cfg.getParameter<int>("maxNJets")),
0091       maxNComb_(cfg.getParameter<int>("maxNComb")),
0092       maxNrIter_(cfg.getParameter<unsigned int>("maxNrIter")),
0093       maxDeltaS_(cfg.getParameter<double>("maxDeltaS")),
0094       maxF_(cfg.getParameter<double>("maxF")),
0095       jetParam_(cfg.getParameter<unsigned>("jetParametrisation")),
0096       constraints_(cfg.getParameter<std::vector<unsigned> >("constraints")),
0097       mW_(cfg.getParameter<double>("mW")),
0098       mTop_(cfg.getParameter<double>("mTop")),
0099       jetEnergyResolutionScaleFactors_(cfg.getParameter<std::vector<double> >("jetEnergyResolutionScaleFactors")),
0100       jetEnergyResolutionEtaBinning_(cfg.getParameter<std::vector<double> >("jetEnergyResolutionEtaBinning")) {
0101   if (cfg.exists("udscResolutions") && cfg.exists("bResolutions")) {
0102     udscResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet> >("udscResolutions");
0103     bResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet> >("bResolutions");
0104   } else if (cfg.exists("udscResolutions") || cfg.exists("bResolutions")) {
0105     if (cfg.exists("udscResolutions"))
0106       throw cms::Exception("Configuration")
0107           << "Parameter 'bResolutions' is needed if parameter 'udscResolutions' is defined!\n";
0108     else
0109       throw cms::Exception("Configuration")
0110           << "Parameter 'udscResolutions' is needed if parameter 'bResolutions' is defined!\n";
0111   }
0112
0113   // define kinematic fit interface
0114   kinFitter = std::make_unique<TtFullHadKinFitter::KinFit>(useBTagging_,
0115                                                            bTags_,
0116                                                            bTagAlgo_,
0117                                                            minBTagValueBJet_,
0118                                                            maxBTagValueNonBJet_,
0119                                                            udscResolutions_,
0120                                                            bResolutions_,
0121                                                            jetEnergyResolutionScaleFactors_,
0122                                                            jetEnergyResolutionEtaBinning_,
0123                                                            jetCorrectionLevel_,
0124                                                            maxNJets_,
0125                                                            maxNComb_,
0126                                                            maxNrIter_,
0127                                                            maxDeltaS_,
0128                                                            maxF_,
0129                                                            jetParam_,
0130                                                            constraints_,
0131                                                            mW_,
0132                                                            mTop_);
0133
0134   // produces the following collections
0135   produces<std::vector<pat::Particle> >("PartonsB");
0136   produces<std::vector<pat::Particle> >("PartonsBBar");
0137   produces<std::vector<pat::Particle> >("PartonsLightQ");
0138   produces<std::vector<pat::Particle> >("PartonsLightQBar");
0139   produces<std::vector<pat::Particle> >("PartonsLightP");
0140   produces<std::vector<pat::Particle> >("PartonsLightPBar");
0141
0142   produces<std::vector<std::vector<int> > >();
0143   produces<std::vector<double> >("Chi2");
0144   produces<std::vector<double> >("Prob");
0145   produces<std::vector<int> >("Status");
0146 }
0147
0148 /// produce fitted object collections and meta data describing fit quality
0149 void TtFullHadKinFitProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0150   // get jet collection
0151   const edm::Handle<std::vector<pat::Jet> >& jets = event.getHandle(jetsToken_);
0152
0153   // get match in case that useOnlyMatch_ is true
0154   std::vector<int> match;
0155   bool invalidMatch = false;
0156   if (useOnlyMatch_) {
0157     kinFitter->setUseOnlyMatch(true);
0158     // in case that only a ceratin match should be used, get match here
0159     const edm::Handle<std::vector<std::vector<int> > >& matches = event.getHandle(matchToken_);
0160     match = *(matches->begin());
0161     // check if match is valid
0162     if (match.size() != nPartons) {
0163       invalidMatch = true;
0164     } else {
0165       for (unsigned int idx = 0; idx < match.size(); ++idx) {
0166         if (match[idx] < 0 || match[idx] >= (int)jets->size()) {
0167           invalidMatch = true;
0168           break;
0169         }
0170       }
0171     }
0172     /// set match to be used
0173     kinFitter->setMatch(match);
0174   }
0175
0176   /// set the validity of a match
0177   kinFitter->setMatchInvalidity(invalidMatch);
0178
0179   std::list<TtFullHadKinFitter::KinFitResult> fitResults = kinFitter->fit(*jets);
0180
0181   // pointer for output collections
0182   std::unique_ptr<std::vector<pat::Particle> > pPartonsB(new std::vector<pat::Particle>);
0183   std::unique_ptr<std::vector<pat::Particle> > pPartonsBBar(new std::vector<pat::Particle>);
0184   std::unique_ptr<std::vector<pat::Particle> > pPartonsLightQ(new std::vector<pat::Particle>);
0185   std::unique_ptr<std::vector<pat::Particle> > pPartonsLightQBar(new std::vector<pat::Particle>);
0186   std::unique_ptr<std::vector<pat::Particle> > pPartonsLightP(new std::vector<pat::Particle>);
0187   std::unique_ptr<std::vector<pat::Particle> > pPartonsLightPBar(new std::vector<pat::Particle>);
0188   // pointer for meta information
0189   std::unique_ptr<std::vector<std::vector<int> > > pCombi(new std::vector<std::vector<int> >);
0190   std::unique_ptr<std::vector<double> > pChi2(new std::vector<double>);
0191   std::unique_ptr<std::vector<double> > pProb(new std::vector<double>);
0192   std::unique_ptr<std::vector<int> > pStatus(new std::vector<int>);
0193
0194   unsigned int iComb = 0;
0195   for (std::list<TtFullHadKinFitter::KinFitResult>::const_iterator res = fitResults.begin(); res != fitResults.end();
0196        ++res) {
0197     if (maxNComb_ >= 1 && iComb == (unsigned int)maxNComb_) {
0198       break;
0199     }
0200     ++iComb;
0201
0202     pPartonsB->push_back(res->B);
0203     pPartonsBBar->push_back(res->BBar);
0204     pPartonsLightQ->push_back(res->LightQ);
0205     pPartonsLightQBar->push_back(res->LightQBar);
0206     pPartonsLightP->push_back(res->LightP);
0207     pPartonsLightPBar->push_back(res->LightPBar);
0208
0209     pCombi->push_back(res->JetCombi);
0210     pChi2->push_back(res->Chi2);
0211     pProb->push_back(res->Prob);
0212     pStatus->push_back(res->Status);
0213   }
0214
0215   event.put(std::move(pCombi));
0216   event.put(std::move(pPartonsB), "PartonsB");
0217   event.put(std::move(pPartonsBBar), "PartonsBBar");
0218   event.put(std::move(pPartonsLightQ), "PartonsLightQ");
0219   event.put(std::move(pPartonsLightQBar), "PartonsLightQBar");
0220   event.put(std::move(pPartonsLightP), "PartonsLightP");
0221   event.put(std::move(pPartonsLightPBar), "PartonsLightPBar");
0222   event.put(std::move(pChi2), "Chi2");
0223   event.put(std::move(pProb), "Prob");
0224   event.put(std::move(pStatus), "Status");
0225 }
0226
0227 #include "FWCore/Framework/interface/MakerMacros.h"
0228 DEFINE_FWK_MODULE(TtFullHadKinFitProducer);