File indexing completed on 2023-03-17 10:43:11
0001
0002 #include <atomic>
0003 #include <cmath>
0004 #include <memory>
0005 #include <string>
0006 #include <vector>
0007
0008
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
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
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
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
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
0337
0338
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
0360 DEFINE_FWK_MODULE(HcalIsoTrackAnalyzer);
