Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 10:43:11

0001 // system include files
0002 #include <atomic>
0003 #include <cmath>
0004 #include <memory>
0005 #include <string>
0006 #include <vector>
0007 
0008 // Root objects
0009 #include "TROOT.h"
0010 #include "TSystem.h"
0011 #include "TFile.h"
0012 #include "TProfile.h"
0013 #include "TDirectory.h"
0014 #include "TTree.h"
0015 #include "TLorentzVector.h"
0016 #include "TInterpreter.h"
0017 
0018 #include "DataFormats/HcalCalibObjects/interface/HcalIsoTrkCalibVariables.h"
0019 #include "DataFormats/HcalCalibObjects/interface/HcalIsoTrkEventVariables.h"
0020 
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0023 
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029 
0030 //#define EDM_ML_DEBUG
0031 
0032 class HcalIsoTrackAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0033 public:
0034   explicit HcalIsoTrackAnalyzer(edm::ParameterSet const&);
0035   ~HcalIsoTrackAnalyzer() override {}
0036 
0037   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0038 
0039 private:
0040   void analyze(edm::Event const&, edm::EventSetup const&) override;
0041   void beginJob() override;
0042   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0043   void endRun(edm::Run const&, edm::EventSetup const&) override;
0044 
0045   const double pTrackLow_, pTrackHigh_;
0046   const int useRaw_, dataType_;
0047   const edm::InputTag labelIsoTkVar_, labelIsoTkEvt_;
0048   const std::vector<int> debEvents_;
0049   edm::EDGetTokenT<HcalIsoTrkCalibVariablesCollection> tokIsoTrkVar_;
0050   edm::EDGetTokenT<HcalIsoTrkEventVariablesCollection> tokIsoTrkEvt_;
0051   unsigned int nRun_, nRange_, nLow_, nHigh_;
0052 
0053   TTree *tree, *tree2;
0054   int t_Run, t_Event, t_DataType, t_ieta, t_iphi;
0055   int t_goodPV, t_nVtx, t_nTrk;
0056   double t_EventWeight, t_p, t_pt, t_phi;
0057   double t_l1pt, t_l1eta, t_l1phi;
0058   double t_l3pt, t_l3eta, t_l3phi;
0059   double t_mindR1, t_mindR2;
0060   double t_eMipDR, t_eMipDR2, t_eMipDR3, t_eMipDR4;
0061   double t_eMipDR5, t_hmaxNearP, t_gentrackP;
0062   double t_emaxNearP, t_eAnnular, t_hAnnular;
0063   double t_eHcal, t_eHcal10, t_eHcal30, t_rhoh;
0064   bool t_selectTk, t_qltyFlag, t_qltyMissFlag, t_qltyPVFlag;
0065   std::vector<unsigned int> t_DetIds, t_DetIds1, t_DetIds3;
0066   std::vector<double> t_HitEnergies, t_HitEnergies1, t_HitEnergies3;
0067   std::vector<bool> t_trgbits;
0068 
0069   unsigned int t_RunNo, t_EventNo;
0070   bool t_TrigPass, t_TrigPassSel, t_L1Bit;
0071   int t_Tracks, t_TracksProp, t_TracksSaved;
0072   int t_TracksLoose, t_TracksTight, t_allvertex;
0073   std::vector<int> t_ietaAll, t_ietaGood, t_trackType;
0074   std::vector<bool> t_hltbits;
0075 };
0076 
0077 HcalIsoTrackAnalyzer::HcalIsoTrackAnalyzer(const edm::ParameterSet& iConfig)
0078     : pTrackLow_(iConfig.getParameter<double>("momentumLow")),
0079       pTrackHigh_(iConfig.getParameter<double>("momentumHigh")),
0080       useRaw_(iConfig.getUntrackedParameter<int>("useRaw", 0)),
0081       dataType_(iConfig.getUntrackedParameter<int>("dataType", 0)),
0082       labelIsoTkVar_(iConfig.getParameter<edm::InputTag>("isoTrackVarLabel")),
0083       labelIsoTkEvt_(iConfig.getParameter<edm::InputTag>("isoTrackEvtLabel")),
0084       debEvents_(iConfig.getParameter<std::vector<int>>("debugEvents")),
0085       tokIsoTrkVar_(consumes<HcalIsoTrkCalibVariablesCollection>(labelIsoTkVar_)),
0086       tokIsoTrkEvt_(consumes<HcalIsoTrkEventVariablesCollection>(labelIsoTkEvt_)),
0087       nRun_(0),
0088       nRange_(0),
0089       nLow_(0),
0090       nHigh_(0) {
0091   usesResource(TFileService::kSharedResource);
0092 
0093   //now do whatever initialization is needed
0094   edm::LogVerbatim("HcalIsoTrack") << "Labels used " << labelIsoTkVar_ << " " << labelIsoTkEvt_;
0095 
0096   edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n\t momentumLow_ " << pTrackLow_
0097                                    << "\t momentumHigh_ " << pTrackHigh_ << "\t useRaw_ " << useRaw_
0098                                    << "\t dataType_      " << dataType_ << " and " << debEvents_.size()
0099                                    << " events to be debugged";
0100 }
0101 
0102 void HcalIsoTrackAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0103   t_Run = iEvent.id().run();
0104   t_Event = iEvent.id().event();
0105   t_DataType = dataType_;
0106 #ifdef EDM_ML_DEBUG
0107   bool debug = (debEvents_.empty())
0108                    ? true
0109                    : (std::find(debEvents_.begin(), debEvents_.end(), iEvent.id().event()) != debEvents_.end());
0110   if (debug)
0111     edm::LogVerbatim("HcalIsoTrack") << "Run " << t_Run << " Event " << t_Event << " type " << t_DataType
0112                                      << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0113                                      << iEvent.bunchCrossing();
0114 #endif
0115 
0116   // Fill from IsoTrkCalibVariables collection
0117   auto const& isotrkCalibColl = iEvent.getHandle(tokIsoTrkVar_);
0118   if (isotrkCalibColl.isValid()) {
0119     auto isotrkCalib = *isotrkCalibColl.product();
0120 #ifdef EDM_ML_DEBUG
0121     if (debug)
0122       edm::LogVerbatim("HcalIsoTrack") << "Finds HcalIsoTrkCalibVariablesCollection with " << isotrkCalib.size()
0123                                        << " entries";
0124     int k(0);
0125 #endif
0126     for (const auto& itr : isotrkCalib) {
0127       t_ieta = itr.ieta_;
0128       t_iphi = itr.iphi_;
0129       t_goodPV = itr.goodPV_;
0130       t_nVtx = itr.nVtx_;
0131       t_nTrk = itr.nTrk_;
0132       t_EventWeight = itr.eventWeight_;
0133       t_p = itr.p_;
0134       t_pt = itr.pt_;
0135       t_phi = itr.phi_;
0136 #ifdef EDM_ML_DEBUG
0137       ++k;
0138       if (debug)
0139         edm::LogVerbatim("HcalIsoTrack") << "Track " << k << " p:pt:phi " << t_p << ":" << t_pt << ":" << t_phi
0140                                          << " nvtx:ntrk:goodPV:wt " << t_nVtx << ":" << t_nTrk << ":" << t_goodPV << ":"
0141                                          << t_EventWeight << " ieta:iphi " << t_ieta << ":" << t_iphi;
0142 #endif
0143       t_l1pt = itr.l1pt_;
0144       t_l1eta = itr.l1eta_;
0145       t_l1phi = itr.l1phi_;
0146       t_l3pt = itr.l3pt_;
0147       t_l3eta = itr.l3eta_;
0148       t_l3phi = itr.l3phi_;
0149       t_mindR1 = itr.mindR1_;
0150       t_mindR2 = itr.mindR2_;
0151 #ifdef EDM_ML_DEBUG
0152       if (debug)
0153         edm::LogVerbatim("HcalIsoTrack") << "L1 pt:eta:phi " << t_l1pt << ":" << t_l1eta << ":" << t_l1phi
0154                                          << " L3 pt:eta:phi " << t_l3pt << ":" << t_l3eta << ":" << t_l3phi << " R1:R2 "
0155                                          << t_mindR1 << ":" << t_mindR2;
0156 #endif
0157       t_eMipDR = itr.eMipDR_[0];
0158       t_eMipDR2 = itr.eMipDR_[1];
0159       t_eMipDR3 = itr.eMipDR_[2];
0160       t_eMipDR4 = itr.eMipDR_[3];
0161       t_eMipDR5 = itr.eMipDR_[4];
0162 #ifdef EDM_ML_DEBUG
0163       if (debug)
0164         edm::LogVerbatim("HcalIsoTrack") << "eMIPDR 1:2:3:4:5 " << t_eMipDR << ":" << t_eMipDR2 << ":" << t_eMipDR3
0165                                          << ":" << t_eMipDR4 << ":" << t_eMipDR5;
0166 #endif
0167       t_hmaxNearP = itr.hmaxNearP_;
0168       t_emaxNearP = itr.emaxNearP_;
0169       t_eAnnular = itr.eAnnular_;
0170       t_hAnnular = itr.hAnnular_;
0171 #ifdef EDM_ML_DEBUG
0172       if (debug)
0173         edm::LogVerbatim("HcalIsoTrack") << "emaxNearP:hmaxNearP " << t_emaxNearP << ":" << t_hmaxNearP
0174                                          << " eAnnlar:hAnnular" << t_eAnnular << ":" << t_hAnnular;
0175 #endif
0176       t_gentrackP = itr.gentrackP_;
0177       t_rhoh = itr.rhoh_;
0178       t_selectTk = itr.selectTk_;
0179       t_qltyFlag = itr.qltyFlag_;
0180       t_qltyMissFlag = itr.qltyMissFlag_;
0181       t_qltyPVFlag = itr.qltyPVFlag_;
0182 #ifdef EDM_ML_DEBUG
0183       if (debug)
0184         edm::LogVerbatim("HcalIsoTrack") << "gentrackP " << t_gentrackP << " rhoh " << t_rhoh
0185                                          << " qltyFlag:qltyMissFlag:qltyPVFlag:selectTk " << t_qltyFlag << ":"
0186                                          << t_qltyMissFlag << ":" << t_qltyPVFlag << ":" << t_selectTk;
0187 #endif
0188       t_trgbits = itr.trgbits_;
0189       t_DetIds = itr.detIds_;
0190       t_DetIds1 = itr.detIds1_;
0191       t_DetIds3 = itr.detIds3_;
0192       if (useRaw_ == 1) {
0193         t_eHcal = itr.eHcalAux_;
0194         t_eHcal10 = itr.eHcal10Aux_;
0195         t_eHcal30 = itr.eHcal30Aux_;
0196         t_HitEnergies = itr.hitEnergiesAux_;
0197         t_HitEnergies1 = itr.hitEnergies1Aux_;
0198         t_HitEnergies3 = itr.hitEnergies3Aux_;
0199       } else if (useRaw_ == 2) {
0200         t_eHcal = itr.eHcalRaw_;
0201         t_eHcal10 = itr.eHcal10Raw_;
0202         t_eHcal30 = itr.eHcal30Raw_;
0203         t_HitEnergies = itr.hitEnergiesRaw_;
0204         t_HitEnergies1 = itr.hitEnergies1Raw_;
0205         t_HitEnergies3 = itr.hitEnergies3Raw_;
0206       } else {
0207         t_eHcal = itr.eHcal_;
0208         t_eHcal10 = itr.eHcal10_;
0209         t_eHcal30 = itr.eHcal30_;
0210         t_HitEnergies = itr.hitEnergies_;
0211         t_HitEnergies1 = itr.hitEnergies1_;
0212         t_HitEnergies3 = itr.hitEnergies3_;
0213       }
0214 #ifdef EDM_ML_DEBUG
0215       if (debug)
0216         edm::LogVerbatim("HcalIsoTrack") << "eHcal:eHcal10:eHCal30 " << t_eHcal << ":" << t_eHcal10 << t_eHcal30;
0217 #endif
0218       tree->Fill();
0219       edm::LogVerbatim("HcalIsoTrackX") << "Run " << t_Run << " Event " << t_Event << " p " << t_p;
0220 
0221       if (t_p < pTrackLow_) {
0222         ++nLow_;
0223       } else if (t_p < pTrackHigh_) {
0224         ++nHigh_;
0225       } else {
0226         ++nRange_;
0227       }
0228     }
0229   } else {
0230     edm::LogVerbatim("HcalIsoTrack") << "Cannot find HcalIsoTrkCalibVariablesCollection";
0231   }
0232 
0233   // Fill from IsoTrkEventVariables collection
0234   auto const& isotrkEventColl = iEvent.getHandle(tokIsoTrkEvt_);
0235   if (isotrkEventColl.isValid()) {
0236     auto isotrkEvent = isotrkEventColl.product();
0237 #ifdef EDM_ML_DEBUG
0238     if (debug)
0239       edm::LogVerbatim("HcalIsoTrack") << "Finds HcalIsoTrkEventVariablesCollection with " << isotrkEvent->size()
0240                                        << " entries";
0241 #endif
0242     auto itr = isotrkEvent->begin();
0243     if (itr != isotrkEvent->end()) {
0244       t_RunNo = iEvent.id().run();
0245       t_EventNo = iEvent.id().event();
0246       t_TrigPass = itr->trigPass_;
0247       t_TrigPassSel = itr->trigPassSel_;
0248       t_L1Bit = itr->l1Bit_;
0249       t_Tracks = itr->tracks_;
0250       t_TracksProp = itr->tracksProp_;
0251       t_TracksSaved = itr->tracksSaved_;
0252       t_TracksLoose = itr->tracksLoose_;
0253       t_TracksTight = itr->tracksTight_;
0254       t_allvertex = itr->allvertex_;
0255       t_ietaAll = itr->ietaAll_;
0256       t_ietaGood = itr->ietaGood_;
0257       t_trackType = itr->trackType_;
0258       t_hltbits = itr->hltbits_;
0259       tree2->Fill();
0260     }
0261   } else {
0262     edm::LogVerbatim("HcalIsoTrack") << "Cannot find HcalIsoTrkEventVariablesCollections";
0263   }
0264 }
0265 
0266 void HcalIsoTrackAnalyzer::beginJob() {
0267   edm::Service<TFileService> fs;
0268   tree = fs->make<TTree>("CalibTree", "CalibTree");
0269 
0270   tree->Branch("t_Run", &t_Run, "t_Run/I");
0271   tree->Branch("t_Event", &t_Event, "t_Event/I");
0272   tree->Branch("t_DataType", &t_DataType, "t_DataType/I");
0273   tree->Branch("t_ieta", &t_ieta, "t_ieta/I");
0274   tree->Branch("t_iphi", &t_iphi, "t_iphi/I");
0275   tree->Branch("t_EventWeight", &t_EventWeight, "t_EventWeight/D");
0276   tree->Branch("t_nVtx", &t_nVtx, "t_nVtx/I");
0277   tree->Branch("t_nTrk", &t_nTrk, "t_nTrk/I");
0278   tree->Branch("t_goodPV", &t_goodPV, "t_goodPV/I");
0279   tree->Branch("t_l1pt", &t_l1pt, "t_l1pt/D");
0280   tree->Branch("t_l1eta", &t_l1eta, "t_l1eta/D");
0281   tree->Branch("t_l1phi", &t_l1phi, "t_l1phi/D");
0282   tree->Branch("t_l3pt", &t_l3pt, "t_l3pt/D");
0283   tree->Branch("t_l3eta", &t_l3eta, "t_l3eta/D");
0284   tree->Branch("t_l3phi", &t_l3phi, "t_l3phi/D");
0285   tree->Branch("t_p", &t_p, "t_p/D");
0286   tree->Branch("t_pt", &t_pt, "t_pt/D");
0287   tree->Branch("t_phi", &t_phi, "t_phi/D");
0288   tree->Branch("t_mindR1", &t_mindR1, "t_mindR1/D");
0289   tree->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
0290   tree->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
0291   tree->Branch("t_eMipDR2", &t_eMipDR2, "t_eMipDR2/D");
0292   tree->Branch("t_eMipDR3", &t_eMipDR3, "t_eMipDR3/D");
0293   tree->Branch("t_eMipDR4", &t_eMipDR4, "t_eMipDR4/D");
0294   tree->Branch("t_eMipDR5", &t_eMipDR5, "t_eMipDR5/D");
0295   tree->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
0296   tree->Branch("t_eHcal10", &t_eHcal10, "t_eHcal10/D");
0297   tree->Branch("t_eHcal30", &t_eHcal30, "t_eHcal30/D");
0298   tree->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
0299   tree->Branch("t_emaxNearP", &t_emaxNearP, "t_emaxNearP/D");
0300   tree->Branch("t_eAnnular", &t_eAnnular, "t_eAnnular/D");
0301   tree->Branch("t_hAnnular", &t_hAnnular, "t_hAnnular/D");
0302   tree->Branch("t_rhoh", &t_rhoh, "t_rhoh/D");
0303   tree->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
0304   tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
0305   tree->Branch("t_qltyMissFlag", &t_qltyMissFlag, "t_qltyMissFlag/O");
0306   tree->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O");
0307   tree->Branch("t_gentrackP", &t_gentrackP, "t_gentrackP/D");
0308 
0309   tree->Branch("t_DetIds", &t_DetIds);
0310   tree->Branch("t_HitEnergies", &t_HitEnergies);
0311   tree->Branch("t_trgbits", &t_trgbits);
0312   tree->Branch("t_DetIds1", &t_DetIds1);
0313   tree->Branch("t_DetIds3", &t_DetIds3);
0314   tree->Branch("t_HitEnergies1", &t_HitEnergies1);
0315   tree->Branch("t_HitEnergies3", &t_HitEnergies3);
0316 
0317   tree2 = fs->make<TTree>("EventInfo", "Event Information");
0318 
0319   tree2->Branch("t_RunNo", &t_RunNo, "t_RunNo/i");
0320   tree2->Branch("t_EventNo", &t_EventNo, "t_EventNo/i");
0321   tree2->Branch("t_Tracks", &t_Tracks, "t_Tracks/I");
0322   tree2->Branch("t_TracksProp", &t_TracksProp, "t_TracksProp/I");
0323   tree2->Branch("t_TracksSaved", &t_TracksSaved, "t_TracksSaved/I");
0324   tree2->Branch("t_TracksLoose", &t_TracksLoose, "t_TracksLoose/I");
0325   tree2->Branch("t_TracksTight", &t_TracksTight, "t_TracksTight/I");
0326   tree2->Branch("t_TrigPass", &t_TrigPass, "t_TrigPass/O");
0327   tree2->Branch("t_TrigPassSel", &t_TrigPassSel, "t_TrigPassSel/O");
0328   tree2->Branch("t_L1Bit", &t_L1Bit, "t_L1Bit/O");
0329   tree2->Branch("t_allvertex", &t_allvertex, "t_allvertex/I");
0330   tree2->Branch("t_ietaAll", &t_ietaAll);
0331   tree2->Branch("t_ietaGood", &t_ietaGood);
0332   tree2->Branch("t_trackType", &t_trackType);
0333   tree2->Branch("t_hltbits", &t_hltbits);
0334 }
0335 
0336 // ------------ method called when starting to processes a run  ------------
0337 
0338 // ------------ method called when ending the processing of a run  ------------
0339 void HcalIsoTrackAnalyzer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0340   nRun_++;
0341   edm::LogVerbatim("HcalIsoTrack") << "endRun[" << nRun_ << "] " << iRun.run() << " with " << nLow_
0342                                    << " events with p < " << pTrackLow_ << ", " << nHigh_ << " events with p > "
0343                                    << pTrackHigh_ << ", and " << nRange_ << " events in the right momentum range";
0344 }
0345 
0346 void HcalIsoTrackAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0347   edm::ParameterSetDescription desc;
0348   desc.add<double>("momentumLow", 40.0);
0349   desc.add<double>("momentumHigh", 60.0);
0350   desc.addUntracked<int>("useRaw", 0);
0351   desc.addUntracked<int>("dataType", 0);
0352   desc.add<edm::InputTag>("isoTrackVarLabel", edm::InputTag("alcaHcalIsotrkProducer", "HcalIsoTrack"));
0353   desc.add<edm::InputTag>("isoTrackEvtLabel", edm::InputTag("alcaHcalIsotrkProducer", "HcalIsoTrackEvent"));
0354   std::vector<int> events;
0355   desc.add<std::vector<int>>("debugEvents", events);
0356   descriptions.add("hcalIsoTrackAnalyzer", desc);
0357 }
0358 
0359 //define this as a plug-in
0360 DEFINE_FWK_MODULE(HcalIsoTrackAnalyzer);