Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:12 

0001 #include "PhysicsTools/TagAndProbe/interface/BaseTreeFiller.h"
0002 #include "FWCore/ServiceRegistry/interface/Service.h"
0003 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0004 
0005 #include <TList.h>
0006 #include <TObjString.h>
0007 
0008 #include <iostream>
0009 using namespace std;
0010 
0011 tnp::ProbeVariable::~ProbeVariable() {}
0012 
0013 tnp::ProbeFlag::~ProbeFlag() {}
0014 
0015 void tnp::ProbeFlag::init(const edm::Event &iEvent) const {
0016   if (external_) {
0017     edm::Handle<edm::View<reco::Candidate> > view;
0018     iEvent.getByToken(srcToken_, view);
0019     passingProbes_.clear();
0020     for (size_t i = 0, n = view->size(); i < n; ++i)
0021       passingProbes_.push_back(view->refAt(i));
0022   }
0023 }
0024 
0025 void tnp::ProbeFlag::fill(const reco::CandidateBaseRef &probe) const {
0026   if (external_) {
0027     value_ = (std::find(passingProbes_.begin(), passingProbes_.end(), probe) != passingProbes_.end());
0028   } else {
0029     value_ = bool(cut_(*probe));
0030   }
0031 }
0032 
0033 tnp::BaseTreeFiller::BaseTreeFiller(const char *name, const edm::ParameterSet &iConfig, edm::ConsumesCollector &iC) {
0034   // make trees as requested
0035   edm::Service<TFileService> fs;
0036   tree_ = fs->make<TTree>(name, name);
0037 
0038   // add the branches
0039   addBranches_(tree_, iConfig, iC, "");
0040 
0041   // set up weights, if needed
0042   if (iConfig.existsAs<double>("eventWeight")) {
0043     weightMode_ = Fixed;
0044     weight_ = iConfig.getParameter<double>("eventWeight");
0045   } else if (iConfig.existsAs<edm::InputTag>("eventWeight")) {
0046     weightMode_ = External;
0047     weightSrcToken_ = iC.consumes<GenEventInfoProduct>(iConfig.getParameter<edm::InputTag>("eventWeight"));
0048     tree_->Branch("psWeight", &psWeight_, "psWeight[5]/F");
0049   } else {
0050     weightMode_ = None;
0051   }
0052   if (weightMode_ != None) {
0053     tree_->Branch("weight", &weight_, "weight/F");
0054     tree_->Branch("totWeight", &totWeight_, "totWeight/F");
0055   }
0056 
0057   LHEinfo_ = iConfig.existsAs<edm::InputTag>("LHEWeightSrc");
0058   if (LHEinfo_) {
0059     _LHECollection = iC.consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("LHEWeightSrc"));
0060     tree_->Branch("lheWeight", &lheWeight_, "lheWeight[9]/F");
0061     tree_->Branch("lhe_ht", &lhe_ht_, "lhe_ht/F");
0062   }
0063 
0064   storePUweight_ = iConfig.existsAs<edm::InputTag>("PUWeightSrc") ? true : false;
0065   if (storePUweight_) {
0066     PUweightSrcToken_ = iC.consumes<double>(iConfig.getParameter<edm::InputTag>("PUWeightSrc"));
0067     tree_->Branch("PUweight", &PUweight_, "PUweight/F");
0068   }
0069 
0070   if (iConfig.existsAs<edm::InputTag>("pileupInfoTag"))
0071     pileupInfoToken_ =
0072         iC.consumes<std::vector<PileupSummaryInfo> >(iConfig.getParameter<edm::InputTag>("pileupInfoTag"));
0073 
0074   addRunLumiInfo_ = iConfig.existsAs<bool>("addRunLumiInfo") ? iConfig.getParameter<bool>("addRunLumiInfo") : false;
0075   if (addRunLumiInfo_) {
0076     tree_->Branch("run", &run_, "run/i");
0077     tree_->Branch("lumi", &lumi_, "lumi/i");
0078     tree_->Branch("event", &event_, "event/l");
0079     tree_->Branch("truePU", &truePU_, "truePU/I");
0080   }
0081   addEventVariablesInfo_ =
0082       iConfig.existsAs<bool>("addEventVariablesInfo") ? iConfig.getParameter<bool>("addEventVariablesInfo") : false;
0083   if (addEventVariablesInfo_) {
0084     /// FC (EGM) - add possibility to customize collections (can run other miniAOD)
0085     edm::InputTag bsIT = iConfig.existsAs<edm::InputTag>("beamSpot") ? iConfig.getParameter<edm::InputTag>("beamSpot")
0086                                                                      : edm::InputTag("offlineBeamSpot");
0087     edm::InputTag vtxIT = iConfig.existsAs<edm::InputTag>("vertexCollection")
0088                               ? iConfig.getParameter<edm::InputTag>("vertexCollection")
0089                               : edm::InputTag("offlinePrimaryVertices");
0090     edm::InputTag pfMetIT = iConfig.existsAs<edm::InputTag>("pfMet") ? iConfig.getParameter<edm::InputTag>("pfMet")
0091                                                                      : edm::InputTag("pfMet");
0092     edm::InputTag tcMetIT = iConfig.existsAs<edm::InputTag>("tcMet") ? iConfig.getParameter<edm::InputTag>("tcMet")
0093                                                                      : edm::InputTag("tcMet");
0094     edm::InputTag clMetIT =
0095         iConfig.existsAs<edm::InputTag>("clMet") ? iConfig.getParameter<edm::InputTag>("clMet") : edm::InputTag("met");
0096 
0097     recVtxsToken_ = iC.consumes<reco::VertexCollection>(vtxIT);
0098     beamSpotToken_ = iC.consumes<reco::BeamSpot>(bsIT);
0099     pfmetToken_ = iC.mayConsume<reco::PFMETCollection>(pfMetIT);
0100     pfmetTokenMiniAOD_ = iC.mayConsume<pat::METCollection>(pfMetIT);
0101     addCaloMet_ = iConfig.existsAs<bool>("addCaloMet") ? iConfig.getParameter<bool>("addCaloMet") : true;
0102     tree_->Branch("event_nPV", &mNPV_, "mNPV/I");
0103     if (addCaloMet_) {
0104       metToken_ = iC.mayConsume<reco::CaloMETCollection>(clMetIT);
0105       tcmetToken_ = iC.mayConsume<reco::METCollection>(tcMetIT);
0106       tree_->Branch("event_met_calomet", &mMET_, "mMET/F");
0107       tree_->Branch("event_met_calosumet", &mSumET_, "mSumET/F");
0108       tree_->Branch("event_met_calometsignificance", &mMETSign_, "mMETSign/F");
0109       tree_->Branch("event_met_tcmet", &mtcMET_, "mtcMET/F");
0110       tree_->Branch("event_met_tcsumet", &mtcSumET_, "mtcSumET/F");
0111       tree_->Branch("event_met_tcmetsignificance", &mtcMETSign_, "mtcMETSign/F");
0112     }
0113     tree_->Branch("event_met_pfmet", &mpfMET_, "mpfMET/F");
0114     tree_->Branch("event_met_pfphi", &mpfPhi_, "mpfPhi/F");
0115     tree_->Branch("event_met_pfsumet", &mpfSumET_, "mpfSumET/F");
0116 
0117     tree_->Branch("event_met_pfmetsignificance", &mpfMETSign_, "mpfMETSign/F");
0118     tree_->Branch("event_PrimaryVertex_x", &mPVx_, "mPVx/F");
0119     tree_->Branch("event_PrimaryVertex_y", &mPVy_, "mPVy/F");
0120     tree_->Branch("event_PrimaryVertex_z", &mPVz_, "mPVz/F");
0121     tree_->Branch("event_BeamSpot_x", &mBSx_, "mBSx/F");
0122     tree_->Branch("event_BeamSpot_y", &mBSy_, "mBSy/F");
0123     tree_->Branch("event_BeamSpot_z", &mBSz_, "mBSz/F");
0124   }
0125 
0126   addRho_ = iConfig.existsAs<edm::InputTag>("rho") ? true : false;
0127   if (addRho_) {
0128     rhoToken_ = iC.consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
0129     tree_->Branch("event_rho", &rho_, "rho/F");
0130   }
0131 }
0132 
0133 tnp::BaseTreeFiller::BaseTreeFiller(BaseTreeFiller &main,
0134                                     const edm::ParameterSet &iConfig,
0135                                     edm::ConsumesCollector &&iC,
0136                                     const std::string &branchNamePrefix)
0137     : addRunLumiInfo_(false), addEventVariablesInfo_(false), tree_(nullptr) {
0138   addRunLumiInfo_ = main.addRunLumiInfo_;
0139   storePUweight_ = main.storePUweight_;
0140   addBranches_(main.tree_, iConfig, iC, branchNamePrefix);
0141 }
0142 
0143 void tnp::BaseTreeFiller::addBranches_(TTree *tree,
0144                                        const edm::ParameterSet &iConfig,
0145                                        edm::ConsumesCollector &iC,
0146                                        const std::string &branchNamePrefix) {
0147   // set up variables
0148   edm::ParameterSet variables = iConfig.getParameter<edm::ParameterSet>("variables");
0149   //.. the ones that are strings
0150   std::vector<std::string> stringVars = variables.getParameterNamesForType<std::string>();
0151   for (std::vector<std::string>::const_iterator it = stringVars.begin(), ed = stringVars.end(); it != ed; ++it) {
0152     vars_.push_back(tnp::ProbeVariable(branchNamePrefix + *it, variables.getParameter<std::string>(*it)));
0153   }
0154   //.. the ones that are InputTags
0155   std::vector<std::string> inputTagVars = variables.getParameterNamesForType<edm::InputTag>();
0156   for (std::vector<std::string>::const_iterator it = inputTagVars.begin(), ed = inputTagVars.end(); it != ed; ++it) {
0157     vars_.push_back(tnp::ProbeVariable(branchNamePrefix + *it,
0158                                        iC.consumes<edm::ValueMap<float> >(variables.getParameter<edm::InputTag>(*it))));
0159   }
0160   // set up flags
0161   edm::ParameterSet flags = iConfig.getParameter<edm::ParameterSet>("flags");
0162   //.. the ones that are strings
0163   std::vector<std::string> stringFlags = flags.getParameterNamesForType<std::string>();
0164   for (std::vector<std::string>::const_iterator it = stringFlags.begin(), ed = stringFlags.end(); it != ed; ++it) {
0165     flags_.push_back(tnp::ProbeFlag(branchNamePrefix + *it, flags.getParameter<std::string>(*it)));
0166   }
0167   //.. the ones that are InputTags
0168   std::vector<std::string> inputTagFlags = flags.getParameterNamesForType<edm::InputTag>();
0169   for (std::vector<std::string>::const_iterator it = inputTagFlags.begin(), ed = inputTagFlags.end(); it != ed; ++it) {
0170     flags_.push_back(tnp::ProbeFlag(branchNamePrefix + *it,
0171                                     iC.consumes<edm::View<reco::Candidate> >(flags.getParameter<edm::InputTag>(*it))));
0172   }
0173 
0174   // then make all the variables in the trees
0175   for (std::vector<tnp::ProbeVariable>::iterator it = vars_.begin(), ed = vars_.end(); it != ed; ++it) {
0176     tree->Branch(it->name().c_str(), it->address(), (it->name() + "/F").c_str());
0177   }
0178 
0179   for (std::vector<tnp::ProbeFlag>::iterator it = flags_.begin(), ed = flags_.end(); it != ed; ++it) {
0180     tree->Branch(it->name().c_str(), it->address(), (it->name() + "/I").c_str());
0181   }
0182 }
0183 
0184 tnp::BaseTreeFiller::~BaseTreeFiller() {}
0185 
0186 void tnp::BaseTreeFiller::init(const edm::Event &iEvent) const {
0187   run_ = iEvent.id().run();
0188   lumi_ = iEvent.id().luminosityBlock();
0189   event_ = iEvent.id().event();
0190 
0191   truePU_ = 0;
0192   if (!iEvent.isRealData() and !pileupInfoToken_.isUninitialized()) {
0193     edm::Handle<std::vector<PileupSummaryInfo> > PupInfo;
0194     iEvent.getByToken(pileupInfoToken_, PupInfo);
0195     truePU_ = PupInfo->begin()->getTrueNumInteractions();
0196   }
0197 
0198   totWeight_ = 1.;
0199   for (std::vector<tnp::ProbeVariable>::const_iterator it = vars_.begin(), ed = vars_.end(); it != ed; ++it) {
0200     it->init(iEvent);
0201   }
0202   for (std::vector<tnp::ProbeFlag>::const_iterator it = flags_.begin(), ed = flags_.end(); it != ed; ++it) {
0203     it->init(iEvent);
0204   }
0205   for (int i = 0; i < 5; i++) {
0206     psWeight_[i] = 1.;  // init
0207   }
0208   if (weightMode_ == External) {
0209     // edm::Handle<double> weight;
0210     //        iEvent.getByToken(weightSrcToken_, weight);
0211     //        weight_ = *weight;
0212     edm::Handle<GenEventInfoProduct> weight;
0213     iEvent.getByToken(weightSrcToken_, weight);
0214     weight_ = weight->weight();
0215     totWeight_ *= weight_;
0216     if (weight->weights().size() >= 10) {
0217       int k = 1;
0218       for (int i = 6; i < 10; i++) {
0219         // hardcoded Pythia 8 isrDefHi,fsrDefHi,isrDefLo,fsrDefLo
0220         psWeight_[k] = weight->weights().at(i) / weight->weight();
0221         k++;
0222       }
0223     }
0224   }
0225 
0226   for (unsigned int i = 0; i < 9; i++) {
0227     lheWeight_[i] = 1.;  // init
0228   }
0229   lhe_ht_ = 0.;
0230   if (LHEinfo_ and !_LHECollection.isUninitialized()) {
0231     edm::Handle<LHEEventProduct> lheEventHandle;
0232     iEvent.getByToken(_LHECollection, lheEventHandle);
0233     for (unsigned int i = 0; i < 9; i++) {
0234       lheWeight_[i] = lheEventHandle->weights().at(i).wgt / lheEventHandle->originalXWGTUP();
0235     }
0236     for (int i = 0; i < lheEventHandle->hepeup().NUP; i++) {
0237       int id = lheEventHandle->hepeup().IDUP[i];
0238       int st = lheEventHandle->hepeup().ISTUP[i];
0239 
0240       // calculate HT at LHE level
0241       if ((abs(id) < 6 || id == 21) && st > 0) {
0242         lhe_ht_ += sqrt(pow(lheEventHandle->hepeup().PUP[i][0], 2) + pow(lheEventHandle->hepeup().PUP[i][1], 2));
0243       }
0244     }
0245   }
0246 
0247   ///// ********** Pileup weight: needed for MC re-weighting for PU *************
0248   PUweight_ = 1;
0249   if (storePUweight_ and !PUweightSrcToken_.isUninitialized()) {
0250     edm::Handle<double> weightPU;
0251     bool isPresent = iEvent.getByToken(PUweightSrcToken_, weightPU);
0252     if (isPresent)
0253       PUweight_ = float(*weightPU);
0254     totWeight_ *= PUweight_;
0255   }
0256 
0257   if (addEventVariablesInfo_) {
0258     /// *********** store some event variables: MET, SumET ******
0259     //////////// Primary vertex //////////////
0260     edm::Handle<reco::VertexCollection> recVtxs;
0261     iEvent.getByToken(recVtxsToken_, recVtxs);
0262     mNPV_ = 0;
0263     mPVx_ = 100.0;
0264     mPVy_ = 100.0;
0265     mPVz_ = 100.0;
0266 
0267     for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
0268       if (!((*recVtxs)[ind].isFake()) && ((*recVtxs)[ind].ndof() > 4) && (fabs((*recVtxs)[ind].z()) <= 24.0) &&
0269           ((*recVtxs)[ind].position().Rho() <= 2.0)) {
0270         mNPV_++;
0271         if (mNPV_ == 1) {  // store the first good primary vertex
0272           mPVx_ = (*recVtxs)[ind].x();
0273           mPVy_ = (*recVtxs)[ind].y();
0274           mPVz_ = (*recVtxs)[ind].z();
0275         }
0276       }
0277     }
0278 
0279     //////////// Beam spot //////////////
0280     edm::Handle<reco::BeamSpot> beamSpot;
0281     iEvent.getByToken(beamSpotToken_, beamSpot);
0282     mBSx_ = beamSpot->position().X();
0283     mBSy_ = beamSpot->position().Y();
0284     mBSz_ = beamSpot->position().Z();
0285 
0286     if (addCaloMet_) {
0287       ////////////// CaloMET //////
0288       edm::Handle<reco::CaloMETCollection> met;
0289       iEvent.getByToken(metToken_, met);
0290       if (met->empty()) {
0291         mMET_ = -1;
0292         mSumET_ = -1;
0293         mMETSign_ = -1;
0294       } else {
0295         mMET_ = (*met)[0].et();
0296         mSumET_ = (*met)[0].sumEt();
0297         mMETSign_ = (*met)[0].significance();
0298       }
0299 
0300       /////// TcMET information /////
0301       edm::Handle<reco::METCollection> tcmet;
0302       iEvent.getByToken(tcmetToken_, tcmet);
0303       if (tcmet->empty()) {
0304         mtcMET_ = -1;
0305         mtcSumET_ = -1;
0306         mtcMETSign_ = -1;
0307       } else {
0308         mtcMET_ = (*tcmet)[0].et();
0309         mtcSumET_ = (*tcmet)[0].sumEt();
0310         mtcMETSign_ = (*tcmet)[0].significance();
0311       }
0312     }
0313 
0314     /////// PfMET information /////
0315     edm::Handle<reco::PFMETCollection> pfmet;
0316     iEvent.getByToken(pfmetToken_, pfmet);
0317     if (pfmet.isValid()) {
0318       if (pfmet->empty()) {
0319         mpfMET_ = -1;
0320         mpfSumET_ = -1;
0321         mpfMETSign_ = -1;
0322       } else {
0323         mpfMET_ = (*pfmet)[0].et();
0324         mpfPhi_ = (*pfmet)[0].phi();
0325         mpfSumET_ = (*pfmet)[0].sumEt();
0326         mpfMETSign_ = (*pfmet)[0].significance();
0327       }
0328     } else {
0329       edm::Handle<pat::METCollection> pfmet2;
0330       iEvent.getByToken(pfmetTokenMiniAOD_, pfmet2);
0331       const pat::MET &met = pfmet2->front();
0332       mpfMET_ = met.pt();
0333       mpfPhi_ = met.phi();
0334       mpfSumET_ = met.sumEt();
0335       mpfMETSign_ = met.significance();
0336     }
0337 
0338     if (addRho_) {
0339       edm::Handle<double> rhos;
0340       iEvent.getByToken(rhoToken_, rhos);
0341       rho_ = (float)*rhos;
0342     }
0343   }
0344 }
0345 
0346 void tnp::BaseTreeFiller::fill(const reco::CandidateBaseRef &probe) const {
0347   for (auto const &var : vars_)
0348     var.fill(probe);
0349   for (auto const &flag : flags_)
0350     flag.fill(probe);
0351 
0352   if (tree_)
0353     tree_->Fill();
0354 }
0355 void tnp::BaseTreeFiller::writeProvenance(const edm::ParameterSet &pset) const {
0356   TList *list = tree_->GetUserInfo();
0357   list->Add(new TObjString(pset.dump().c_str()));
0358 }
[ 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!