File indexing completed on 2022-02-28 01:32:03
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 "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0019 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0020
0021 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0022 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0023 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0024
0025
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/TrackReco/interface/HitPattern.h"
0029 #include "DataFormats/TrackReco/interface/TrackBase.h"
0030
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032 #include "DataFormats/VertexReco/interface/Vertex.h"
0033 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0034
0035 #include "DataFormats/MuonReco/interface/Muon.h"
0036 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0037 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0038
0039
0040 #include "DataFormats/Common/interface/TriggerResults.h"
0041 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0042 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0043 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0044 #include "DataFormats/Math/interface/LorentzVector.h"
0045 #include "DataFormats/Math/interface/deltaR.h"
0046
0047 #include "L1Trigger/L1TGlobal/interface/L1TGlobalUtil.h"
0048 #include "DataFormats/L1TGlobal/interface/GlobalAlgBlk.h"
0049 #include "DataFormats/L1TGlobal/interface/GlobalExtBlk.h"
0050 #include "DataFormats/L1Trigger/interface/BXVector.h"
0051
0052 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0053 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0054 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0055
0056
0057 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0058 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0059 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0060
0061 #include "FWCore/Framework/interface/Frameworkfwd.h"
0062 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0063
0064 #include "FWCore/Framework/interface/Event.h"
0065 #include "FWCore/Framework/interface/MakerMacros.h"
0066 #include "FWCore/ServiceRegistry/interface/Service.h"
0067 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0068 #include "FWCore/Common/interface/TriggerNames.h"
0069 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0070 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0071
0072 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0073 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0074 #include "Calibration/IsolatedParticles/interface/eCone.h"
0075 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0076 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0077 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0078
0079 #include "MagneticField/Engine/interface/MagneticField.h"
0080 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0081 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0082 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0083 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0084 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0085 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0086 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0087 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0088 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0089 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0090 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0091
0092
0093
0094 class HcalIsoTrkAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0095 public:
0096 explicit HcalIsoTrkAnalyzer(edm::ParameterSet const&);
0097 ~HcalIsoTrkAnalyzer() override {}
0098
0099 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0100
0101 private:
0102 void analyze(edm::Event const&, edm::EventSetup const&) override;
0103 void beginJob() override;
0104 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0105 void endRun(edm::Run const&, edm::EventSetup const&) override;
0106
0107 std::array<int, 3> fillTree(std::vector<math::XYZTLorentzVector>& vecL1,
0108 std::vector<math::XYZTLorentzVector>& vecL3,
0109 math::XYZPoint& leadPV,
0110 std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
0111 std::vector<spr::propagatedTrackID>& trkCaloDets,
0112 const CaloGeometry* geo,
0113 const CaloTopology* topo,
0114 const HcalTopology* theHBHETopology,
0115 const EcalChannelStatus* theEcalChStatus,
0116 const EcalSeverityLevelAlgo* theEcalSevlv,
0117 edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0118 edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0119 edm::Handle<HBHERecHitCollection>& hbhe,
0120 edm::Handle<CaloTowerCollection>& towerHandle,
0121 edm::Handle<reco::GenParticleCollection>& genParticles,
0122 const HcalRespCorrs* respCorrs,
0123 const edm::Handle<reco::MuonCollection>& muonh);
0124 double dR(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
0125 double trackP(const reco::Track*, const edm::Handle<reco::GenParticleCollection>&);
0126 double rhoh(const edm::Handle<CaloTowerCollection>&);
0127 double eThreshold(const DetId& id, const CaloGeometry* geo) const;
0128 DetId newId(const DetId&);
0129 void storeEnergy(int indx,
0130 const HcalRespCorrs* respCorrs,
0131 const std::vector<DetId>& ids,
0132 std::vector<double>& edet,
0133 double& eHcal,
0134 std::vector<unsigned int>* detIds,
0135 std::vector<double>* hitEnergies);
0136 bool notaMuon(const reco::Track* pTrack0, const edm::Handle<reco::MuonCollection>& muonh);
0137
0138 l1t::L1TGlobalUtil* l1GtUtils_;
0139 HLTConfigProvider hltConfig_;
0140 const std::vector<std::string> trigNames_;
0141 spr::trackSelectionParameters selectionParameter_;
0142 const std::string theTrackQuality_;
0143 const std::string processName_, l1Filter_;
0144 const std::string l2Filter_, l3Filter_;
0145 const double a_coneR_, a_mipR_, a_mipR2_, a_mipR3_;
0146 const double a_mipR4_, a_mipR5_, pTrackMin_, eEcalMax_;
0147 const double maxRestrictionP_, slopeRestrictionP_;
0148 const double hcalScale_, eIsolate1_, eIsolate2_;
0149 const double pTrackLow_, pTrackHigh_;
0150 const int prescaleLow_, prescaleHigh_;
0151 const int useRaw_, dataType_, mode_;
0152 const bool ignoreTrigger_, useL1Trigger_;
0153 const bool unCorrect_, collapseDepth_;
0154 const double hitEthrEB_, hitEthrEE0_, hitEthrEE1_;
0155 const double hitEthrEE2_, hitEthrEE3_;
0156 const double hitEthrEELo_, hitEthrEEHi_;
0157 const edm::InputTag triggerEvent_, theTriggerResultsLabel_;
0158 const std::string labelGenTrack_, labelRecVtx_, labelEB_;
0159 const std::string labelEE_, labelHBHE_, labelTower_, l1TrigName_;
0160 const std::vector<int> oldID_, newDepth_;
0161 const bool hep17_;
0162 const std::vector<int> debEvents_;
0163 const bool usePFThresh_;
0164 const std::string labelBS_, modnam_, prdnam_, labelMuon_;
0165 const edm::InputTag algTag_, extTag_;
0166
0167 const edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0168 const edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0169 const edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
0170 const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0171 const edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0172 const edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0173 const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0174 const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0175 const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0176 const edm::EDGetTokenT<CaloTowerCollection> tok_cala_;
0177 const edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0178 const edm::EDGetTokenT<BXVector<GlobalAlgBlk>> tok_alg_;
0179 const edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0180
0181 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0182 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0183 const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0184 const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0185 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0186 const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0187 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0188 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_resp_;
0189 const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ecalPFRecHitThresholds_;
0190
0191 unsigned int nRun_, nLow_, nHigh_;
0192 double a_charIsoR_, a_coneR1_, a_coneR2_;
0193 const HcalDDDRecConstants* hdc_;
0194 const EcalPFRecHitThresholds* eThresholds_;
0195
0196 std::vector<double> etabins_, phibins_;
0197 std::vector<int> oldDet_, oldEta_, oldDepth_;
0198 double etadist_, phidist_, etahalfdist_, phihalfdist_;
0199
0200 TTree *tree, *tree2;
0201 unsigned int t_RunNo, t_EventNo;
0202 int t_Run, t_Event, t_DataType, t_ieta, t_iphi;
0203 int t_goodPV, t_nVtx, t_nTrk;
0204 double t_EventWeight, t_p, t_pt, t_phi;
0205 double t_l1pt, t_l1eta, t_l1phi;
0206 double t_l3pt, t_l3eta, t_l3phi;
0207 double t_mindR1, t_mindR2;
0208 double t_eMipDR, t_eMipDR2, t_eMipDR3, t_eMipDR4;
0209 double t_eMipDR5, t_hmaxNearP, t_gentrackP;
0210 double t_emaxNearP, t_eAnnular, t_hAnnular;
0211 double t_eHcal, t_eHcal10, t_eHcal30, t_rhoh;
0212 bool t_selectTk, t_qltyFlag, t_qltyMissFlag;
0213 bool t_qltyPVFlag, t_TrigPass, t_TrigPassSel;
0214 std::vector<unsigned int>*t_DetIds, *t_DetIds1, *t_DetIds3;
0215 std::vector<double>*t_HitEnergies, *t_HitEnergies1, *t_HitEnergies3;
0216 std::vector<bool>*t_trgbits, *t_hltbits;
0217 bool t_L1Bit;
0218 int t_Tracks, t_TracksProp, t_TracksSaved;
0219 int t_TracksLoose, t_TracksTight, t_allvertex;
0220 std::vector<int>*t_ietaAll, *t_ietaGood, *t_trackType;
0221
0222 bool debug_;
0223 };
0224
0225 HcalIsoTrkAnalyzer::HcalIsoTrkAnalyzer(const edm::ParameterSet& iConfig)
0226 : trigNames_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0227 theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0228 processName_(iConfig.getParameter<std::string>("processName")),
0229 l1Filter_(iConfig.getParameter<std::string>("l1Filter")),
0230 l2Filter_(iConfig.getParameter<std::string>("l2Filter")),
0231 l3Filter_(iConfig.getParameter<std::string>("l3Filter")),
0232 a_coneR_(iConfig.getParameter<double>("coneRadius")),
0233 a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0234 a_mipR2_(iConfig.getParameter<double>("coneRadiusMIP2")),
0235 a_mipR3_(iConfig.getParameter<double>("coneRadiusMIP3")),
0236 a_mipR4_(iConfig.getParameter<double>("coneRadiusMIP4")),
0237 a_mipR5_(iConfig.getParameter<double>("coneRadiusMIP5")),
0238 pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
0239 eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
0240 maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
0241 slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
0242 hcalScale_(iConfig.getUntrackedParameter<double>("hHcalScale", 1.0)),
0243 eIsolate1_(iConfig.getParameter<double>("isolationEnergyTight")),
0244 eIsolate2_(iConfig.getParameter<double>("isolationEnergyLoose")),
0245 pTrackLow_(iConfig.getParameter<double>("momentumLow")),
0246 pTrackHigh_(iConfig.getParameter<double>("momentumHigh")),
0247 prescaleLow_(iConfig.getParameter<int>("prescaleLow")),
0248 prescaleHigh_(iConfig.getParameter<int>("prescaleHigh")),
0249 useRaw_(iConfig.getUntrackedParameter<int>("useRaw", 0)),
0250 dataType_(iConfig.getUntrackedParameter<int>("dataType", 0)),
0251 mode_(iConfig.getUntrackedParameter<int>("outMode", 11)),
0252 ignoreTrigger_(iConfig.getUntrackedParameter<bool>("ignoreTriggers", false)),
0253 useL1Trigger_(iConfig.getUntrackedParameter<bool>("useL1Trigger", false)),
0254 unCorrect_(iConfig.getUntrackedParameter<bool>("unCorrect", false)),
0255 collapseDepth_(iConfig.getUntrackedParameter<bool>("collapseDepth", false)),
0256 hitEthrEB_(iConfig.getParameter<double>("EBHitEnergyThreshold")),
0257 hitEthrEE0_(iConfig.getParameter<double>("EEHitEnergyThreshold0")),
0258 hitEthrEE1_(iConfig.getParameter<double>("EEHitEnergyThreshold1")),
0259 hitEthrEE2_(iConfig.getParameter<double>("EEHitEnergyThreshold2")),
0260 hitEthrEE3_(iConfig.getParameter<double>("EEHitEnergyThreshold3")),
0261 hitEthrEELo_(iConfig.getParameter<double>("EEHitEnergyThresholdLow")),
0262 hitEthrEEHi_(iConfig.getParameter<double>("EEHitEnergyThresholdHigh")),
0263 triggerEvent_(iConfig.getParameter<edm::InputTag>("labelTriggerEvent")),
0264 theTriggerResultsLabel_(iConfig.getParameter<edm::InputTag>("labelTriggerResult")),
0265 labelGenTrack_(iConfig.getParameter<std::string>("labelTrack")),
0266 labelRecVtx_(iConfig.getParameter<std::string>("labelVertex")),
0267 labelEB_(iConfig.getParameter<std::string>("labelEBRecHit")),
0268 labelEE_(iConfig.getParameter<std::string>("labelEERecHit")),
0269 labelHBHE_(iConfig.getParameter<std::string>("labelHBHERecHit")),
0270 labelTower_(iConfig.getParameter<std::string>("labelCaloTower")),
0271 l1TrigName_(iConfig.getUntrackedParameter<std::string>("l1TrigName", "L1_SingleJet60")),
0272 oldID_(iConfig.getUntrackedParameter<std::vector<int>>("oldID")),
0273 newDepth_(iConfig.getUntrackedParameter<std::vector<int>>("newDepth")),
0274 hep17_(iConfig.getUntrackedParameter<bool>("hep17")),
0275 debEvents_(iConfig.getParameter<std::vector<int>>("debugEvents")),
0276 usePFThresh_(iConfig.getParameter<bool>("usePFThreshold")),
0277 labelBS_(iConfig.getParameter<std::string>("labelBeamSpot")),
0278 modnam_(iConfig.getUntrackedParameter<std::string>("moduleName", "")),
0279 prdnam_(iConfig.getUntrackedParameter<std::string>("producerName", "")),
0280 labelMuon_(iConfig.getParameter<std::string>("labelMuon")),
0281 algTag_(iConfig.getParameter<edm::InputTag>("algInputTag")),
0282 extTag_(iConfig.getParameter<edm::InputTag>("extInputTag")),
0283 tok_trigEvt_(consumes<trigger::TriggerEvent>(triggerEvent_)),
0284 tok_trigRes_(consumes<edm::TriggerResults>(theTriggerResultsLabel_)),
0285 tok_parts_(consumes<reco::GenParticleCollection>(edm::InputTag("genParticles"))),
0286 tok_genTrack_(consumes<reco::TrackCollection>(labelGenTrack_)),
0287 tok_bs_(consumes<reco::BeamSpot>(labelBS_)),
0288 tok_recVtx_((modnam_.empty()) ? consumes<reco::VertexCollection>(labelRecVtx_)
0289 : consumes<reco::VertexCollection>(edm::InputTag(modnam_, labelRecVtx_, prdnam_))),
0290 tok_EB_((modnam_.empty()) ? consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEB_))
0291 : consumes<EcalRecHitCollection>(edm::InputTag(modnam_, labelEB_, prdnam_))),
0292 tok_EE_((modnam_.empty()) ? consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEE_))
0293 : consumes<EcalRecHitCollection>(edm::InputTag(modnam_, labelEE_, prdnam_))),
0294 tok_hbhe_((modnam_.empty()) ? consumes<HBHERecHitCollection>(labelHBHE_)
0295 : consumes<HBHERecHitCollection>(edm::InputTag(modnam_, labelHBHE_, prdnam_))),
0296 tok_cala_(consumes<CaloTowerCollection>(labelTower_)),
0297 tok_ew_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0298 tok_alg_(consumes<BXVector<GlobalAlgBlk>>(algTag_)),
0299 tok_Muon_(consumes<reco::MuonCollection>(labelMuon_)),
0300 tok_ddrec_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0301 tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0302 tok_ecalChStatus_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd>()),
0303 tok_sevlv_(esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>()),
0304 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0305 tok_caloTopology_(esConsumes<CaloTopology, CaloTopologyRecord>()),
0306 tok_htopo_(esConsumes<HcalTopology, HcalRecNumberingRecord>()),
0307 tok_resp_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()),
0308 tok_ecalPFRecHitThresholds_(esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>()),
0309 nRun_(0),
0310 nLow_(0),
0311 nHigh_(0),
0312 hdc_(nullptr) {
0313 usesResource(TFileService::kSharedResource);
0314
0315
0316 const double isolationRadius(28.9), innerR(10.0), outerR(30.0);
0317 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0318 selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0319 ;
0320 selectionParameter_.minQuality = trackQuality_;
0321 selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
0322 selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
0323 selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0324 selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0325 selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0326 selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0327 selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0328 selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
0329 a_charIsoR_ = a_coneR_ + isolationRadius;
0330 a_coneR1_ = a_coneR_ + innerR;
0331 a_coneR2_ = a_coneR_ + outerR;
0332
0333
0334
0335
0336
0337
0338 for (unsigned int k = 0; k < oldID_.size(); ++k) {
0339 oldDet_.emplace_back((oldID_[k] / 10000) % 10);
0340 oldEta_.emplace_back((oldID_[k] / 100) % 100);
0341 oldDepth_.emplace_back(oldID_[k] % 100);
0342 }
0343
0344 l1GtUtils_ =
0345 new l1t::L1TGlobalUtil(iConfig, consumesCollector(), *this, algTag_, extTag_, l1t::UseEventSetupIn::Event);
0346
0347
0348 if (modnam_.empty()) {
0349 edm::LogVerbatim("HcalIsoTrack") << "Labels used " << triggerEvent_ << " " << theTriggerResultsLabel_ << " "
0350 << labelBS_ << " " << labelRecVtx_ << " " << labelGenTrack_ << " "
0351 << edm::InputTag("ecalRecHit", labelEB_) << " "
0352 << edm::InputTag("ecalRecHit", labelEE_) << " " << labelHBHE_ << " " << labelTower_
0353 << " " << labelMuon_;
0354 } else {
0355 edm::LogVerbatim("HcalIsoTrack") << "Labels used " << triggerEvent_ << " " << theTriggerResultsLabel_ << " "
0356 << labelBS_ << " " << edm::InputTag(modnam_, labelRecVtx_, prdnam_) << " "
0357 << labelGenTrack_ << " " << edm::InputTag(modnam_, labelEB_, prdnam_) << " "
0358 << edm::InputTag(modnam_, labelEE_, prdnam_) << " "
0359 << edm::InputTag(modnam_, labelHBHE_, prdnam_) << " " << labelTower_ << " "
0360 << labelMuon_;
0361 }
0362
0363 edm::LogVerbatim("HcalIsoTrack")
0364 << "Parameters read from config file \n"
0365 << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality " << theTrackQuality_ << "\t minQuality "
0366 << selectionParameter_.minQuality << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
0367 << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2 << "\t maxDpOverP "
0368 << selectionParameter_.maxDpOverP << "\t minOuterHit " << selectionParameter_.minOuterHit << "\t minLayerCrossed "
0369 << selectionParameter_.minLayerCrossed << "\t maxInMiss " << selectionParameter_.maxInMiss << "\t maxOutMiss "
0370 << selectionParameter_.maxOutMiss << "\t a_coneR " << a_coneR_ << ":" << a_coneR1_ << ":" << a_coneR2_
0371 << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR " << a_mipR_ << "\t a_mipR2 " << a_mipR2_ << "\t a_mipR3 "
0372 << a_mipR3_ << "\t a_mipR4 " << a_mipR4_ << "\t a_mipR5 " << a_mipR5_ << "\n pTrackMin_ " << pTrackMin_
0373 << "\t eEcalMax_ " << eEcalMax_ << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
0374 << slopeRestrictionP_ << "\t eIsolateStrong_ " << eIsolate1_ << "\t eIsolateSoft_ " << eIsolate2_
0375 << "\t hcalScale_ " << hcalScale_ << "\n\t momentumLow_ " << pTrackLow_ << "\t prescaleLow_ " << prescaleLow_
0376 << "\t momentumHigh_ " << pTrackHigh_ << "\t prescaleHigh_ " << prescaleHigh_ << "\n\t useRaw_ " << useRaw_
0377 << "\t ignoreTrigger_ " << ignoreTrigger_ << "\n\t useL1Trigegr_ " << useL1Trigger_ << "\t dataType_ "
0378 << dataType_ << "\t mode_ " << mode_ << "\t unCorrect_ " << unCorrect_ << "\t collapseDepth_ "
0379 << collapseDepth_ << "\t L1TrigName_ " << l1TrigName_ << "\nThreshold flag used " << usePFThresh_
0380 << " value for EB " << hitEthrEB_ << " EE " << hitEthrEE0_ << ":" << hitEthrEE1_ << ":" << hitEthrEE2_ << ":"
0381 << hitEthrEE3_ << ":" << hitEthrEELo_ << ":" << hitEthrEEHi_ << " and " << debEvents_.size()
0382 << " events to be debugged";
0383 edm::LogVerbatim("HcalIsoTrack") << "Process " << processName_ << " L1Filter:" << l1Filter_
0384 << " L2Filter:" << l2Filter_ << " L3Filter:" << l3Filter_;
0385 for (unsigned int k = 0; k < trigNames_.size(); ++k) {
0386 edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
0387 }
0388 edm::LogVerbatim("HcalIsoTrack") << oldID_.size() << " DetIDs to be corrected with HEP17 flag:" << hep17_;
0389 for (unsigned int k = 0; k < oldID_.size(); ++k)
0390 edm::LogVerbatim("HcalIsoTrack") << "[" << k << "] Det " << oldDet_[k] << " EtaAbs " << oldEta_[k] << " Depth "
0391 << oldDepth_[k] << ":" << newDepth_[k];
0392
0393 for (int i = 0; i < 10; i++)
0394 phibins_.push_back(-M_PI + 0.1 * (2 * i + 1) * M_PI);
0395 for (int i = 0; i < 8; ++i)
0396 etabins_.push_back(-2.1 + 0.6 * i);
0397 etadist_ = etabins_[1] - etabins_[0];
0398 phidist_ = phibins_[1] - phibins_[0];
0399 etahalfdist_ = 0.5 * etadist_;
0400 phihalfdist_ = 0.5 * phidist_;
0401 edm::LogVerbatim("HcalIsoTrack") << "EtaDist " << etadist_ << " " << etahalfdist_ << " PhiDist " << phidist_ << " "
0402 << phihalfdist_;
0403 unsigned int k1(0), k2(0);
0404 for (auto phi : phibins_) {
0405 edm::LogVerbatim("HcalIsoTrack") << "phibin_[" << k1 << "] " << phi;
0406 ++k1;
0407 }
0408 for (auto eta : etabins_) {
0409 edm::LogVerbatim("HcalIsoTrack") << "etabin_[" << k2 << "] " << eta;
0410 ++k2;
0411 }
0412 }
0413
0414 void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0415 t_Run = iEvent.id().run();
0416 t_Event = iEvent.id().event();
0417 t_DataType = dataType_;
0418 debug_ = (debEvents_.empty())
0419 ? true
0420 : (std::find(debEvents_.begin(), debEvents_.end(), iEvent.id().event()) != debEvents_.end());
0421 #ifdef EDM_ML_DEBUG
0422 if (debug_)
0423 edm::LogVerbatim("HcalIsoTrack") << "Run " << t_Run << " Event " << t_Event << " type " << t_DataType
0424 << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0425 << iEvent.bunchCrossing();
0426 #endif
0427
0428 const MagneticField* bField = &iSetup.getData(tok_bFieldH_);
0429 const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0430 const EcalSeverityLevelAlgo* theEcalSevlv = &iSetup.getData(tok_sevlv_);
0431 eThresholds_ = &iSetup.getData(tok_ecalPFRecHitThresholds_);
0432
0433
0434 const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0435 const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0436 const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo_);
0437
0438
0439 const HcalRespCorrs* respCorrs = &iSetup.getData(tok_resp_);
0440
0441
0442 edm::Handle<reco::GenParticleCollection> genParticles = iEvent.getHandle(tok_parts_);
0443
0444 bool okC(true);
0445
0446 edm::Handle<reco::TrackCollection> trkCollection = iEvent.getHandle(tok_genTrack_);
0447 if (!trkCollection.isValid()) {
0448 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0449 okC = false;
0450 }
0451
0452
0453 const edm::Handle<reco::MuonCollection> muonh = iEvent.getHandle(tok_Muon_);
0454
0455
0456 t_EventWeight = 1.0;
0457 const edm::Handle<GenEventInfoProduct> genEventInfo = iEvent.getHandle(tok_ew_);
0458 if (genEventInfo.isValid())
0459 t_EventWeight = genEventInfo->weight();
0460
0461
0462 const edm::Handle<reco::VertexCollection> recVtxs = iEvent.getHandle(tok_recVtx_);
0463 const edm::Handle<reco::BeamSpot> beamSpotH = iEvent.getHandle(tok_bs_);
0464 math::XYZPoint leadPV(0, 0, 0);
0465 t_goodPV = t_nVtx = 0;
0466 if (recVtxs.isValid() && !(recVtxs->empty())) {
0467 t_nVtx = recVtxs->size();
0468 for (unsigned int k = 0; k < recVtxs->size(); ++k) {
0469 if (!((*recVtxs)[k].isFake()) && ((*recVtxs)[k].ndof() > 4)) {
0470 if (t_goodPV == 0)
0471 leadPV = math::XYZPoint((*recVtxs)[k].x(), (*recVtxs)[k].y(), (*recVtxs)[k].z());
0472 t_goodPV++;
0473 }
0474 }
0475 }
0476 if (t_goodPV == 0 && beamSpotH.isValid()) {
0477 leadPV = beamSpotH->position();
0478 }
0479 t_allvertex = t_goodPV;
0480 #ifdef EDM_ML_DEBUG
0481 if (debug_) {
0482 edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV << " out of " << t_goodPV << " vertex";
0483 if (beamSpotH.isValid())
0484 edm::LogVerbatim("HcalIsoTrack") << " Beam Spot " << beamSpotH->position();
0485 }
0486 #endif
0487
0488 edm::Handle<EcalRecHitCollection> barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0489 if (!barrelRecHitsHandle.isValid()) {
0490 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0491 okC = false;
0492 }
0493 edm::Handle<EcalRecHitCollection> endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0494 if (!endcapRecHitsHandle.isValid()) {
0495 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0496 okC = false;
0497 }
0498 edm::Handle<HBHERecHitCollection> hbhe = iEvent.getHandle(tok_hbhe_);
0499 if (!hbhe.isValid()) {
0500 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0501 okC = false;
0502 }
0503 edm::Handle<CaloTowerCollection> caloTower = iEvent.getHandle(tok_cala_);
0504
0505
0506 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0507 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0508 std::vector<spr::propagatedTrackID> trkCaloDets;
0509 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, false);
0510 std::vector<math::XYZTLorentzVector> vecL1, vecL3;
0511 t_RunNo = iEvent.id().run();
0512 t_EventNo = iEvent.id().event();
0513 t_Tracks = trkCollection->size();
0514 t_TracksProp = trkCaloDirections.size();
0515 t_ietaAll->clear();
0516 t_ietaGood->clear();
0517 t_trackType->clear();
0518 t_trgbits->clear();
0519 t_hltbits->clear();
0520 #ifdef EDM_ML_DEBUG
0521 if (debug_)
0522 edm::LogVerbatim("HcalIsoTrack") << "# of propagated tracks " << t_TracksProp << " out of " << t_Tracks
0523 << " with Trigger " << ignoreTrigger_;
0524 #endif
0525
0526
0527 t_trgbits->assign(trigNames_.size(), false);
0528 t_hltbits->assign(trigNames_.size(), false);
0529 t_TracksSaved = t_TracksLoose = t_TracksTight = 0;
0530 t_L1Bit = true;
0531 t_TrigPass = false;
0532
0533 if (!ignoreTrigger_) {
0534
0535 l1GtUtils_->retrieveL1(iEvent, iSetup, tok_alg_);
0536 const std::vector<std::pair<std::string, bool>>& finalDecisions = l1GtUtils_->decisionsFinal();
0537 for (const auto& decision : finalDecisions) {
0538 if (decision.first.find(l1TrigName_) != std::string::npos) {
0539 t_L1Bit = decision.second;
0540 break;
0541 }
0542 }
0543 #ifdef EDM_ML_DEBUG
0544 if (debug_)
0545 edm::LogVerbatim("HcalIsoTrack") << "Trigger Information for " << l1TrigName_ << " is " << t_L1Bit
0546 << " from a list of " << finalDecisions.size() << " decisions";
0547 #endif
0548
0549
0550 const edm::Handle<edm::TriggerResults> triggerResults = iEvent.getHandle(tok_trigRes_);
0551 if (triggerResults.isValid()) {
0552 const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0553 const std::vector<std::string>& names = triggerNames.triggerNames();
0554 if (!trigNames_.empty()) {
0555 for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0556 int hlt = triggerResults->accept(iHLT);
0557 for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0558 if (names[iHLT].find(trigNames_[i]) != std::string::npos) {
0559 t_trgbits->at(i) = (hlt > 0);
0560 t_hltbits->at(i) = (hlt > 0);
0561 if (hlt > 0)
0562 t_TrigPass = true;
0563 #ifdef EDM_ML_DEBUG
0564 if (debug_)
0565 edm::LogVerbatim("HcalIsoTrack")
0566 << "This trigger " << names[iHLT] << " Flag " << hlt << ":" << t_trgbits->at(i);
0567 #endif
0568 }
0569 }
0570 }
0571 }
0572 }
0573 #ifdef EDM_ML_DEBUG
0574 if (debug_)
0575 edm::LogVerbatim("HcalIsoTrack") << "HLT Information shows " << t_TrigPass << ":" << trigNames_.empty() << ":"
0576 << okC;
0577 #endif
0578 }
0579
0580 std::array<int, 3> ntksave{{0, 0, 0}};
0581 if (ignoreTrigger_ || useL1Trigger_) {
0582 t_l1pt = t_l1eta = t_l1phi = 0;
0583 t_l3pt = t_l3eta = t_l3phi = 0;
0584 if (ignoreTrigger_ || t_L1Bit)
0585 ntksave = fillTree(vecL1,
0586 vecL3,
0587 leadPV,
0588 trkCaloDirections,
0589 trkCaloDets,
0590 geo,
0591 caloTopology,
0592 theHBHETopology,
0593 theEcalChStatus,
0594 theEcalSevlv,
0595 barrelRecHitsHandle,
0596 endcapRecHitsHandle,
0597 hbhe,
0598 caloTower,
0599 genParticles,
0600 respCorrs,
0601 muonh);
0602 t_TracksSaved = ntksave[0];
0603 t_TracksLoose = ntksave[1];
0604 t_TracksTight = ntksave[2];
0605 } else {
0606 trigger::TriggerEvent triggerEvent;
0607 const edm::Handle<trigger::TriggerEvent>& triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
0608 if (!triggerEventHandle.isValid()) {
0609 edm::LogWarning("HcalIsoTrack") << "Error! Can't get the product " << triggerEvent_.label();
0610 } else if (okC) {
0611 const edm::Handle<edm::TriggerResults> triggerResults = iEvent.getHandle(tok_trigRes_);
0612 triggerEvent = *(triggerEventHandle.product());
0613 const trigger::TriggerObjectCollection& TOC(triggerEvent.getObjects());
0614 bool done(false);
0615 if (triggerResults.isValid()) {
0616 std::vector<std::string> modules;
0617 const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0618 const std::vector<std::string>& names = triggerNames.triggerNames();
0619 for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0620 bool ok = (t_TrigPass) || (trigNames_.empty());
0621 if (ok) {
0622 unsigned int triggerindx = hltConfig_.triggerIndex(names[iHLT]);
0623 const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
0624 std::vector<math::XYZTLorentzVector> vecL2;
0625 vecL1.clear();
0626 vecL3.clear();
0627
0628 for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0629 std::vector<int> Keys;
0630 std::string label = triggerEvent.filterTag(ifilter).label();
0631
0632 for (unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
0633 if (label.find(moduleLabels[imodule]) != std::string::npos) {
0634 #ifdef EDM_ML_DEBUG
0635 if (debug_)
0636 edm::LogVerbatim("HcalIsoTrack") << "FilterName " << label;
0637 #endif
0638 for (unsigned int ifiltrKey = 0; ifiltrKey < triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
0639 Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
0640 const trigger::TriggerObject& TO(TOC[Keys[ifiltrKey]]);
0641 math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
0642 if (label.find(l2Filter_) != std::string::npos) {
0643 vecL2.push_back(v4);
0644 } else if (label.find(l3Filter_) != std::string::npos) {
0645 vecL3.push_back(v4);
0646 } else if ((label.find(l1Filter_) != std::string::npos) || (l1Filter_.empty())) {
0647 vecL1.push_back(v4);
0648 }
0649 #ifdef EDM_ML_DEBUG
0650 if (debug_)
0651 edm::LogVerbatim("HcalIsoTrack")
0652 << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta() << " phi " << TO.phi()
0653 << " mass " << TO.mass() << " Id " << TO.id();
0654 #endif
0655 }
0656 #ifdef EDM_ML_DEBUG
0657 if (debug_)
0658 edm::LogVerbatim("HcalIsoTrack")
0659 << "sizes " << vecL1.size() << ":" << vecL2.size() << ":" << vecL3.size();
0660 #endif
0661 }
0662 }
0663 }
0664
0665 math::XYZTLorentzVector mindRvec1;
0666 double mindR1(999);
0667 for (unsigned int i = 0; i < vecL2.size(); i++) {
0668 double dr = dR(vecL1[0], vecL2[i]);
0669 #ifdef EDM_ML_DEBUG
0670 if (debug_)
0671 edm::LogVerbatim("HcalIsoTrack") << "lvl2[" << i << "] dR " << dr;
0672 #endif
0673 if (dr < mindR1) {
0674 mindR1 = dr;
0675 mindRvec1 = vecL2[i];
0676 }
0677 }
0678 #ifdef EDM_ML_DEBUG
0679 if (debug_)
0680 edm::LogVerbatim("HcalIsoTrack") << "L2 object closest to L1 " << mindRvec1 << " at Dr " << mindR1;
0681 #endif
0682
0683 if (!vecL1.empty()) {
0684 t_l1pt = vecL1[0].pt();
0685 t_l1eta = vecL1[0].eta();
0686 t_l1phi = vecL1[0].phi();
0687 } else {
0688 t_l1pt = t_l1eta = t_l1phi = 0;
0689 }
0690 if (!vecL3.empty()) {
0691 t_l3pt = vecL3[0].pt();
0692 t_l3eta = vecL3[0].eta();
0693 t_l3phi = vecL3[0].phi();
0694 } else {
0695 t_l3pt = t_l3eta = t_l3phi = 0;
0696 }
0697
0698 if (!done) {
0699 ntksave = fillTree(vecL1,
0700 vecL3,
0701 leadPV,
0702 trkCaloDirections,
0703 trkCaloDets,
0704 geo,
0705 caloTopology,
0706 theHBHETopology,
0707 theEcalChStatus,
0708 theEcalSevlv,
0709 barrelRecHitsHandle,
0710 endcapRecHitsHandle,
0711 hbhe,
0712 caloTower,
0713 genParticles,
0714 respCorrs,
0715 muonh);
0716 t_TracksSaved += ntksave[0];
0717 t_TracksLoose += ntksave[1];
0718 t_TracksTight += ntksave[2];
0719 done = true;
0720 }
0721 }
0722 }
0723 }
0724 }
0725 }
0726 #ifdef EDM_ML_DEBUG
0727 if (debug_)
0728 edm::LogVerbatim("HcalIsoTrack") << "Final results on selected tracks " << t_TracksSaved << ":" << t_TracksLoose
0729 << ":" << t_TracksTight;
0730 #endif
0731 t_TrigPassSel = (t_TracksSaved > 0);
0732 tree2->Fill();
0733 }
0734
0735 void HcalIsoTrkAnalyzer::beginJob() {
0736 edm::Service<TFileService> fs;
0737 tree = fs->make<TTree>("CalibTree", "CalibTree");
0738
0739 tree->Branch("t_Run", &t_Run, "t_Run/I");
0740 tree->Branch("t_Event", &t_Event, "t_Event/I");
0741 tree->Branch("t_DataType", &t_DataType, "t_DataType/I");
0742 tree->Branch("t_ieta", &t_ieta, "t_ieta/I");
0743 tree->Branch("t_iphi", &t_iphi, "t_iphi/I");
0744 tree->Branch("t_EventWeight", &t_EventWeight, "t_EventWeight/D");
0745 tree->Branch("t_nVtx", &t_nVtx, "t_nVtx/I");
0746 tree->Branch("t_nTrk", &t_nTrk, "t_nTrk/I");
0747 tree->Branch("t_goodPV", &t_goodPV, "t_goodPV/I");
0748 tree->Branch("t_l1pt", &t_l1pt, "t_l1pt/D");
0749 tree->Branch("t_l1eta", &t_l1eta, "t_l1eta/D");
0750 tree->Branch("t_l1phi", &t_l1phi, "t_l1phi/D");
0751 tree->Branch("t_l3pt", &t_l3pt, "t_l3pt/D");
0752 tree->Branch("t_l3eta", &t_l3eta, "t_l3eta/D");
0753 tree->Branch("t_l3phi", &t_l3phi, "t_l3phi/D");
0754 tree->Branch("t_p", &t_p, "t_p/D");
0755 tree->Branch("t_pt", &t_pt, "t_pt/D");
0756 tree->Branch("t_phi", &t_phi, "t_phi/D");
0757 tree->Branch("t_mindR1", &t_mindR1, "t_mindR1/D");
0758 tree->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
0759 tree->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
0760 tree->Branch("t_eMipDR2", &t_eMipDR2, "t_eMipDR2/D");
0761 tree->Branch("t_eMipDR3", &t_eMipDR3, "t_eMipDR3/D");
0762 tree->Branch("t_eMipDR4", &t_eMipDR4, "t_eMipDR4/D");
0763 tree->Branch("t_eMipDR5", &t_eMipDR5, "t_eMipDR5/D");
0764 tree->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
0765 tree->Branch("t_eHcal10", &t_eHcal10, "t_eHcal10/D");
0766 tree->Branch("t_eHcal30", &t_eHcal30, "t_eHcal30/D");
0767 tree->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
0768 tree->Branch("t_emaxNearP", &t_emaxNearP, "t_emaxNearP/D");
0769 tree->Branch("t_eAnnular", &t_eAnnular, "t_eAnnular/D");
0770 tree->Branch("t_hAnnular", &t_hAnnular, "t_hAnnular/D");
0771 tree->Branch("t_rhoh", &t_rhoh, "t_rhoh/D");
0772 tree->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
0773 tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
0774 tree->Branch("t_qltyMissFlag", &t_qltyMissFlag, "t_qltyMissFlag/O");
0775 tree->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O");
0776 tree->Branch("t_gentrackP", &t_gentrackP, "t_gentrackP/D");
0777
0778 t_DetIds = new std::vector<unsigned int>();
0779 t_DetIds1 = new std::vector<unsigned int>();
0780 t_DetIds3 = new std::vector<unsigned int>();
0781 t_HitEnergies = new std::vector<double>();
0782 t_HitEnergies1 = new std::vector<double>();
0783 t_HitEnergies3 = new std::vector<double>();
0784 t_trgbits = new std::vector<bool>();
0785 tree->Branch("t_DetIds", "std::vector<unsigned int>", &t_DetIds);
0786 tree->Branch("t_HitEnergies", "std::vector<double>", &t_HitEnergies);
0787 tree->Branch("t_trgbits", "std::vector<bool>", &t_trgbits);
0788 tree->Branch("t_DetIds1", "std::vector<unsigned int>", &t_DetIds1);
0789 tree->Branch("t_DetIds3", "std::vector<unsigned int>", &t_DetIds3);
0790 tree->Branch("t_HitEnergies1", "std::vector<double>", &t_HitEnergies1);
0791 tree->Branch("t_HitEnergies3", "std::vector<double>", &t_HitEnergies3);
0792
0793 tree2 = fs->make<TTree>("EventInfo", "Event Information");
0794
0795 tree2->Branch("t_RunNo", &t_RunNo, "t_RunNo/i");
0796 tree2->Branch("t_EventNo", &t_EventNo, "t_EventNo/i");
0797 tree2->Branch("t_Tracks", &t_Tracks, "t_Tracks/I");
0798 tree2->Branch("t_TracksProp", &t_TracksProp, "t_TracksProp/I");
0799 tree2->Branch("t_TracksSaved", &t_TracksSaved, "t_TracksSaved/I");
0800 tree2->Branch("t_TracksLoose", &t_TracksLoose, "t_TracksLoose/I");
0801 tree2->Branch("t_TracksTight", &t_TracksTight, "t_TracksTight/I");
0802 tree2->Branch("t_TrigPass", &t_TrigPass, "t_TrigPass/O");
0803 tree2->Branch("t_TrigPassSel", &t_TrigPassSel, "t_TrigPassSel/O");
0804 tree2->Branch("t_L1Bit", &t_L1Bit, "t_L1Bit/O");
0805 tree2->Branch("t_allvertex", &t_allvertex, "t_allvertex/I");
0806 t_hltbits = new std::vector<bool>();
0807 t_ietaAll = new std::vector<int>();
0808 t_ietaGood = new std::vector<int>();
0809 t_trackType = new std::vector<int>();
0810 tree2->Branch("t_ietaAll", "std::vector<int>", &t_ietaAll);
0811 tree2->Branch("t_ietaGood", "std::vector<int>", &t_ietaGood);
0812 tree2->Branch("t_trackType", "std::vector<int>", &t_trackType);
0813 tree2->Branch("t_hltbits", "std::vector<bool>", &t_hltbits);
0814 }
0815
0816
0817 void HcalIsoTrkAnalyzer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0818 hdc_ = &iSetup.getData(tok_ddrec_);
0819
0820 if (!ignoreTrigger_) {
0821 bool changed_(true);
0822 bool flag = hltConfig_.init(iRun, iSetup, processName_, changed_);
0823 edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " process " << processName_
0824 << " init flag " << flag << " change flag " << changed_;
0825
0826 if (changed_) {
0827 #ifdef EDM_ML_DEBUG
0828 edm::LogVerbatim("HcalIsoTrack") << "New trigger menu found !!!";
0829 #endif
0830 const unsigned int n(hltConfig_.size());
0831 for (unsigned itrig = 0; itrig < trigNames_.size(); itrig++) {
0832 unsigned int triggerindx = hltConfig_.triggerIndex(trigNames_[itrig]);
0833 if (triggerindx >= n) {
0834 edm::LogWarning("HcalIsoTrack") << trigNames_[itrig] << " " << triggerindx << " does not exist in "
0835 << "the current menu";
0836 #ifdef EDM_ML_DEBUG
0837 } else {
0838 edm::LogVerbatim("HcalIsoTrack") << trigNames_[itrig] << " " << triggerindx << " exists";
0839 #endif
0840 }
0841 }
0842 }
0843 }
0844 }
0845
0846
0847 void HcalIsoTrkAnalyzer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0848 nRun_++;
0849 edm::LogVerbatim("HcalIsoTrack") << "endRun[" << nRun_ << "] " << iRun.run();
0850 }
0851
0852 void HcalIsoTrkAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0853 edm::ParameterSetDescription desc;
0854 std::vector<std::string> trig = {"HLT_PFJet40",
0855 "HLT_PFJet60",
0856 "HLT_PFJet80",
0857 "HLT_PFJet140",
0858 "HLT_PFJet200",
0859 "HLT_PFJet260",
0860 "HLT_PFJet320",
0861 "HLT_PFJet400",
0862 "HLT_PFJet450",
0863 "HLT_PFJet500"};
0864 desc.add<std::vector<std::string>>("triggers", trig);
0865 desc.add<std::string>("processName", "HLT");
0866 desc.add<std::string>("l1Filter", "");
0867 desc.add<std::string>("l2Filter", "L2Filter");
0868 desc.add<std::string>("l3Filter", "Filter");
0869
0870 desc.add<std::string>("trackQuality", "highPurity");
0871 desc.add<double>("minTrackPt", 1.0);
0872 desc.add<double>("maxDxyPV", 0.02);
0873 desc.add<double>("maxDzPV", 0.02);
0874 desc.add<double>("maxChi2", 5.0);
0875 desc.add<double>("maxDpOverP", 0.1);
0876 desc.add<int>("minOuterHit", 4);
0877 desc.add<int>("minLayerCrossed", 8);
0878 desc.add<int>("maxInMiss", 0);
0879 desc.add<int>("maxOutMiss", 0);
0880
0881 desc.add<double>("minimumTrackP", 10.0);
0882 desc.add<double>("coneRadius", 34.98);
0883
0884 desc.add<double>("coneRadiusMIP", 14.0);
0885 desc.add<double>("coneRadiusMIP2", 18.0);
0886 desc.add<double>("coneRadiusMIP3", 20.0);
0887 desc.add<double>("coneRadiusMIP4", 22.0);
0888 desc.add<double>("coneRadiusMIP5", 24.0);
0889 desc.add<double>("maximumEcalEnergy", 2.0);
0890
0891 desc.add<double>("maxTrackP", 8.0);
0892 desc.add<double>("slopeTrackP", 0.05090504066);
0893 desc.add<double>("isolationEnergyTight", 2.0);
0894 desc.add<double>("isolationEnergyLoose", 10.0);
0895
0896 desc.add<double>("EBHitEnergyThreshold", 0.08);
0897 desc.add<double>("EEHitEnergyThreshold0", 0.30);
0898 desc.add<double>("EEHitEnergyThreshold1", 0.00);
0899 desc.add<double>("EEHitEnergyThreshold2", 0.00);
0900 desc.add<double>("EEHitEnergyThreshold3", 0.00);
0901 desc.add<double>("EEHitEnergyThresholdLow", 0.30);
0902 desc.add<double>("EEHitEnergyThresholdHigh", 0.30);
0903
0904 desc.add<double>("momentumLow", 40.0);
0905 desc.add<double>("momentumHigh", 60.0);
0906 desc.add<int>("prescaleLow", 1);
0907 desc.add<int>("prescaleHigh", 1);
0908
0909 desc.add<edm::InputTag>("labelTriggerEvent", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
0910 desc.add<edm::InputTag>("labelTriggerResult", edm::InputTag("TriggerResults", "", "HLT"));
0911 desc.add<std::string>("labelTrack", "generalTracks");
0912 desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0913 desc.add<std::string>("labelEBRecHit", "EcalRecHitsEB");
0914 desc.add<std::string>("labelEERecHit", "EcalRecHitsEE");
0915 desc.add<std::string>("labelHBHERecHit", "hbhereco");
0916 desc.add<std::string>("labelBeamSpot", "offlineBeamSpot");
0917 desc.add<std::string>("labelCaloTower", "towerMaker");
0918 desc.add<std::string>("labelMuon", "muons");
0919 desc.add<edm::InputTag>("algInputTag", edm::InputTag("gtStage2Digis"));
0920 desc.add<edm::InputTag>("extInputTag", edm::InputTag("gtStage2Digis"));
0921 desc.addUntracked<std::string>("moduleName", "");
0922 desc.addUntracked<std::string>("producerName", "");
0923
0924
0925 desc.addUntracked<int>("useRaw", 0);
0926 desc.addUntracked<bool>("ignoreTriggers", false);
0927 desc.addUntracked<bool>("useL1Trigger", false);
0928 desc.addUntracked<double>("hcalScale", 1.0);
0929 desc.addUntracked<int>("dataType", 0);
0930 desc.addUntracked<bool>("unCorrect", false);
0931 desc.addUntracked<bool>("collapseDepth", false);
0932 desc.addUntracked<std::string>("l1TrigName", "L1_SingleJet60");
0933 desc.addUntracked<int>("outMode", 11);
0934 std::vector<int> dummy;
0935 desc.addUntracked<std::vector<int>>("oldID", dummy);
0936 desc.addUntracked<std::vector<int>>("newDepth", dummy);
0937 desc.addUntracked<bool>("hep17", false);
0938 std::vector<int> events;
0939 desc.add<std::vector<int>>("debugEvents", events);
0940 desc.add<bool>("usePFThreshold", true);
0941 descriptions.add("hcalIsoTrkAnalyzer", desc);
0942 }
0943
0944 std::array<int, 3> HcalIsoTrkAnalyzer::fillTree(std::vector<math::XYZTLorentzVector>& vecL1,
0945 std::vector<math::XYZTLorentzVector>& vecL3,
0946 math::XYZPoint& leadPV,
0947 std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
0948 std::vector<spr::propagatedTrackID>& trkCaloDets,
0949 const CaloGeometry* geo,
0950 const CaloTopology* caloTopology,
0951 const HcalTopology* theHBHETopology,
0952 const EcalChannelStatus* theEcalChStatus,
0953 const EcalSeverityLevelAlgo* theEcalSevlv,
0954 edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0955 edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0956 edm::Handle<HBHERecHitCollection>& hbhe,
0957 edm::Handle<CaloTowerCollection>& tower,
0958 edm::Handle<reco::GenParticleCollection>& genParticles,
0959 const HcalRespCorrs* respCorrs,
0960 const edm::Handle<reco::MuonCollection>& muonh) {
0961 int nSave(0), nLoose(0), nTight(0);
0962
0963 unsigned int nTracks(0), nselTracks(0);
0964 t_nTrk = trkCaloDirections.size();
0965 t_rhoh = (tower.isValid()) ? rhoh(tower) : 0;
0966 for (const auto& trkDetItr : trkCaloDirections) {
0967 const reco::Track* pTrack = &(*(trkDetItr.trkItr));
0968 math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0969 t_p = pTrack->p();
0970 t_pt = pTrack->pt();
0971 t_phi = pTrack->phi();
0972
0973 #ifdef EDM_ML_DEBUG
0974 if (debug_)
0975 edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) : " << t_pt << "|"
0976 << pTrack->eta() << "|" << t_phi << "|" << t_p;
0977 #endif
0978 t_mindR2 = 999;
0979 for (unsigned int k = 0; k < vecL3.size(); ++k) {
0980 double dr = dR(vecL3[k], v4);
0981 if (dr < t_mindR2) {
0982 t_mindR2 = dr;
0983 }
0984 }
0985 t_mindR1 = (!vecL1.empty()) ? dR(vecL1[0], v4) : 999;
0986 #ifdef EDM_ML_DEBUG
0987 if (debug_)
0988 edm::LogVerbatim("HcalIsoTrack") << "Closest L3 object at dr : " << t_mindR2 << " and from L1 " << t_mindR1;
0989 #endif
0990
0991 t_ieta = t_iphi = 0;
0992 if (trkDetItr.okHCAL) {
0993 HcalDetId detId = (HcalDetId)(trkDetItr.detIdHCAL);
0994 t_ieta = detId.ieta();
0995 t_iphi = detId.iphi();
0996 if (t_p > 40.0 && t_p <= 60.0)
0997 t_ietaAll->emplace_back(t_ieta);
0998 }
0999
1000 t_selectTk = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
1001 spr::trackSelectionParameters oneCutParameters = selectionParameter_;
1002 oneCutParameters.maxDxyPV = 10;
1003 oneCutParameters.maxDzPV = 100;
1004 oneCutParameters.maxInMiss = 2;
1005 oneCutParameters.maxOutMiss = 2;
1006 bool qltyFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
1007 oneCutParameters = selectionParameter_;
1008 oneCutParameters.maxDxyPV = 10;
1009 oneCutParameters.maxDzPV = 100;
1010 t_qltyMissFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
1011 oneCutParameters = selectionParameter_;
1012 oneCutParameters.maxInMiss = 2;
1013 oneCutParameters.maxOutMiss = 2;
1014 t_qltyPVFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
1015 double eIsolation = maxRestrictionP_ * exp(slopeRestrictionP_ * std::abs((double)t_ieta));
1016 if (eIsolation < eIsolate1_)
1017 eIsolation = eIsolate1_;
1018 if (eIsolation < eIsolate2_)
1019 eIsolation = eIsolate2_;
1020 #ifdef EDM_ML_DEBUG
1021 if (debug_)
1022 edm::LogVerbatim("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr.okECAL << "|"
1023 << trkDetItr.okHCAL << " eIsolation " << eIsolation;
1024 #endif
1025 t_qltyFlag = (qltyFlag && trkDetItr.okECAL && trkDetItr.okHCAL);
1026 bool notMuon = (muonh.isValid()) ? notaMuon(pTrack, muonh) : true;
1027 if (t_qltyFlag && notMuon) {
1028 nselTracks++;
1029 int nNearTRKs(0);
1030
1031 std::vector<DetId> eIds;
1032 std::vector<double> eHit;
1033 t_eMipDR = spr::eCone_ecal(geo,
1034 barrelRecHitsHandle,
1035 endcapRecHitsHandle,
1036 trkDetItr.pointHCAL,
1037 trkDetItr.pointECAL,
1038 a_mipR_,
1039 trkDetItr.directionECAL,
1040 eIds,
1041 eHit);
1042 double eEcal(0);
1043 for (unsigned int k = 0; k < eIds.size(); ++k) {
1044 if (eHit[k] > eThreshold(eIds[k], geo))
1045 eEcal += eHit[k];
1046 }
1047 #ifdef EDM_ML_DEBUG
1048 if (debug_)
1049 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << t_eMipDR << ":" << eEcal;
1050 #endif
1051 t_eMipDR = eEcal;
1052
1053
1054 std::vector<DetId> eIds2;
1055 std::vector<double> eHit2;
1056 t_eMipDR2 = spr::eCone_ecal(geo,
1057 barrelRecHitsHandle,
1058 endcapRecHitsHandle,
1059 trkDetItr.pointHCAL,
1060 trkDetItr.pointECAL,
1061 a_mipR2_,
1062 trkDetItr.directionECAL,
1063 eIds2,
1064 eHit2);
1065 double eEcal2(0);
1066 for (unsigned int k = 0; k < eIds2.size(); ++k) {
1067 if (eHit2[k] > eThreshold(eIds2[k], geo))
1068 eEcal2 += eHit2[k];
1069 }
1070 #ifdef EDM_ML_DEBUG
1071 if (debug_)
1072 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << t_eMipDR2 << ":" << eEcal2;
1073 #endif
1074 t_eMipDR2 = eEcal2;
1075
1076
1077 std::vector<DetId> eIds3;
1078 std::vector<double> eHit3;
1079 t_eMipDR3 = spr::eCone_ecal(geo,
1080 barrelRecHitsHandle,
1081 endcapRecHitsHandle,
1082 trkDetItr.pointHCAL,
1083 trkDetItr.pointECAL,
1084 a_mipR3_,
1085 trkDetItr.directionECAL,
1086 eIds3,
1087 eHit3);
1088 double eEcal3(0);
1089 for (unsigned int k = 0; k < eIds3.size(); ++k) {
1090 if (eHit3[k] > eThreshold(eIds3[k], geo))
1091 eEcal3 += eHit3[k];
1092 }
1093 #ifdef EDM_ML_DEBUG
1094 if (debug_)
1095 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << t_eMipDR3 << ":" << eEcal3;
1096 #endif
1097 t_eMipDR3 = eEcal3;
1098
1099
1100 std::vector<DetId> eIds4;
1101 std::vector<double> eHit4;
1102 t_eMipDR4 = spr::eCone_ecal(geo,
1103 barrelRecHitsHandle,
1104 endcapRecHitsHandle,
1105 trkDetItr.pointHCAL,
1106 trkDetItr.pointECAL,
1107 a_mipR4_,
1108 trkDetItr.directionECAL,
1109 eIds4,
1110 eHit4);
1111 double eEcal4(0);
1112 for (unsigned int k = 0; k < eIds4.size(); ++k) {
1113 if (eHit4[k] > eThreshold(eIds4[k], geo))
1114 eEcal4 += eHit4[k];
1115 }
1116 #ifdef EDM_ML_DEBUG
1117 if (debug_)
1118 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << t_eMipDR4 << ":" << eEcal4;
1119 #endif
1120 t_eMipDR4 = eEcal4;
1121
1122
1123 std::vector<DetId> eIds5;
1124 std::vector<double> eHit5;
1125 t_eMipDR5 = spr::eCone_ecal(geo,
1126 barrelRecHitsHandle,
1127 endcapRecHitsHandle,
1128 trkDetItr.pointHCAL,
1129 trkDetItr.pointECAL,
1130 a_mipR5_,
1131 trkDetItr.directionECAL,
1132 eIds5,
1133 eHit5);
1134 double eEcal5(0);
1135 for (unsigned int k = 0; k < eIds5.size(); ++k) {
1136 if (eHit5[k] > eThreshold(eIds5[k], geo))
1137 eEcal5 += eHit5[k];
1138 }
1139 #ifdef EDM_ML_DEBUG
1140 if (debug_)
1141 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << t_eMipDR5 << ":" << eEcal5;
1142 #endif
1143 t_eMipDR5 = eEcal5;
1144
1145
1146 t_emaxNearP = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15, 15);
1147 const DetId cellE(trkDetItr.detIdECAL);
1148 std::pair<double, bool> e11x11P = spr::eECALmatrix(cellE,
1149 barrelRecHitsHandle,
1150 endcapRecHitsHandle,
1151 *theEcalChStatus,
1152 geo,
1153 caloTopology,
1154 theEcalSevlv,
1155 5,
1156 5,
1157 -100.0,
1158 -100.0,
1159 -100.0,
1160 100.0);
1161 std::pair<double, bool> e15x15P = spr::eECALmatrix(cellE,
1162 barrelRecHitsHandle,
1163 endcapRecHitsHandle,
1164 *theEcalChStatus,
1165 geo,
1166 caloTopology,
1167 theEcalSevlv,
1168 7,
1169 7,
1170 -100.0,
1171 -100.0,
1172 -100.0,
1173 100.0);
1174 if (e11x11P.second && e15x15P.second) {
1175 t_eAnnular = (e15x15P.first - e11x11P.first);
1176 } else {
1177 t_eAnnular = -(e15x15P.first - e11x11P.first);
1178 }
1179 t_hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
1180 const DetId cellH(trkDetItr.detIdHCAL);
1181 double h5x5 = spr::eHCALmatrix(
1182 theHBHETopology, cellH, hbhe, 2, 2, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1183 double h7x7 = spr::eHCALmatrix(
1184 theHBHETopology, cellH, hbhe, 3, 3, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1185 t_hAnnular = h7x7 - h5x5;
1186 #ifdef EDM_ML_DEBUG
1187 if (debug_)
1188 edm::LogVerbatim("HcalIsoTrack") << "max p Near (Ecal) " << t_emaxNearP << " (Hcal) " << t_hmaxNearP
1189 << " Annular E (Ecal) " << e11x11P.first << ":" << e15x15P.first << ":"
1190 << t_eAnnular << " (Hcal) " << h5x5 << ":" << h7x7 << ":" << t_hAnnular;
1191 #endif
1192 t_gentrackP = trackP(pTrack, genParticles);
1193 if (t_eMipDR < eEcalMax_ && t_hmaxNearP < eIsolation) {
1194 t_DetIds->clear();
1195 t_HitEnergies->clear();
1196 t_DetIds1->clear();
1197 t_HitEnergies1->clear();
1198 t_DetIds3->clear();
1199 t_HitEnergies3->clear();
1200 int nRecHits(-999), nRecHits1(-999), nRecHits3(-999);
1201 std::vector<DetId> ids, ids1, ids3;
1202 std::vector<double> edet0, edet1, edet3;
1203 t_eHcal = spr::eCone_hcal(geo,
1204 hbhe,
1205 trkDetItr.pointHCAL,
1206 trkDetItr.pointECAL,
1207 a_coneR_,
1208 trkDetItr.directionHCAL,
1209 nRecHits,
1210 ids,
1211 edet0,
1212 useRaw_);
1213 if (!oldID_.empty()) {
1214 for (unsigned k = 0; k < ids.size(); ++k)
1215 ids[k] = newId(ids[k]);
1216 }
1217 storeEnergy(0, respCorrs, ids, edet0, t_eHcal, t_DetIds, t_HitEnergies);
1218
1219
1220 t_eHcal10 = spr::eCone_hcal(geo,
1221 hbhe,
1222 trkDetItr.pointHCAL,
1223 trkDetItr.pointECAL,
1224 a_coneR1_,
1225 trkDetItr.directionHCAL,
1226 nRecHits1,
1227 ids1,
1228 edet1,
1229 useRaw_);
1230 if (!oldID_.empty()) {
1231 for (unsigned k = 0; k < ids1.size(); ++k)
1232 ids1[k] = newId(ids1[k]);
1233 }
1234 storeEnergy(1, respCorrs, ids1, edet1, t_eHcal10, t_DetIds1, t_HitEnergies1);
1235
1236
1237 t_eHcal30 = spr::eCone_hcal(geo,
1238 hbhe,
1239 trkDetItr.pointHCAL,
1240 trkDetItr.pointECAL,
1241 a_coneR2_,
1242 trkDetItr.directionHCAL,
1243 nRecHits3,
1244 ids3,
1245 edet3,
1246 useRaw_);
1247 if (!oldID_.empty()) {
1248 for (unsigned k = 0; k < ids3.size(); ++k)
1249 ids3[k] = newId(ids3[k]);
1250 }
1251 storeEnergy(3, respCorrs, ids3, edet3, t_eHcal30, t_DetIds3, t_HitEnergies3);
1252
1253 #ifdef EDM_ML_DEBUG
1254 if (debug_) {
1255 edm::LogVerbatim("HcalIsoTrack")
1256 << "This track : " << nTracks << " (pt|eta|phi|p) : " << t_pt << "|" << pTrack->eta() << "|" << t_phi
1257 << "|" << t_p << " Generator Level p " << t_gentrackP;
1258 edm::LogVerbatim("HcalIsoTrack")
1259 << "e_MIP " << t_eMipDR << " Chg Isolation " << t_hmaxNearP << " eHcal" << t_eHcal << " ieta " << t_ieta
1260 << " Quality " << t_qltyMissFlag << ":" << t_qltyPVFlag << ":" << t_selectTk;
1261 for (unsigned int ll = 0; ll < t_DetIds->size(); ll++) {
1262 edm::LogVerbatim("HcalIsoTrack")
1263 << "det id is = " << HcalDetId(t_DetIds->at(ll)) << " hit enery is = " << t_HitEnergies->at(ll);
1264 }
1265 for (unsigned int ll = 0; ll < t_DetIds1->size(); ll++) {
1266 edm::LogVerbatim("HcalIsoTrack")
1267 << "det id is = " << HcalDetId(t_DetIds1->at(ll)) << " hit enery is = " << t_HitEnergies1->at(ll);
1268 }
1269 for (unsigned int ll = 0; ll < t_DetIds3->size(); ll++) {
1270 edm::LogVerbatim("HcalIsoTrack")
1271 << "det id is = " << HcalDetId(t_DetIds3->at(ll)) << " hit enery is = " << t_HitEnergies3->at(ll);
1272 }
1273 }
1274 #endif
1275 bool accept(false);
1276 if (t_p > pTrackMin_) {
1277 if (t_p < pTrackLow_) {
1278 ++nLow_;
1279 if (prescaleLow_ <= 1)
1280 accept = true;
1281 else if (nLow_ % prescaleLow_ == 1)
1282 accept = true;
1283 } else if (t_p > pTrackHigh_) {
1284 ++nHigh_;
1285 if (prescaleHigh_ <= 1)
1286 accept = true;
1287 else if (nHigh_ % prescaleHigh_ == 1)
1288 accept = true;
1289 } else {
1290 accept = true;
1291 }
1292 }
1293 if (accept) {
1294 tree->Fill();
1295 edm::LogVerbatim("HcalIsoTrackX")
1296 << "Run " << t_RunNo << " Event " << t_EventNo << " Track " << nTracks << " p " << t_p;
1297 nSave++;
1298 int type(0);
1299 if (t_eMipDR < 1.0) {
1300 if (t_hmaxNearP < eIsolate2_) {
1301 ++nLoose;
1302 type = 1;
1303 }
1304 if (t_hmaxNearP < eIsolate1_) {
1305 ++nTight;
1306 type = 2;
1307 }
1308 }
1309 if (t_p > 40.0 && t_p <= 60.0 && t_selectTk) {
1310 t_ietaGood->emplace_back(t_ieta);
1311 t_trackType->emplace_back(type);
1312 }
1313 #ifdef EDM_ML_DEBUG
1314 if (debug_) {
1315 for (unsigned int k = 0; k < t_trgbits->size(); k++) {
1316 edm::LogVerbatim("HcalIsoTrack") << "trigger bit is = " << t_trgbits->at(k);
1317 }
1318 }
1319 #endif
1320 }
1321 }
1322 }
1323 ++nTracks;
1324 }
1325 std::array<int, 3> i3{{nSave, nLoose, nTight}};
1326 return i3;
1327 }
1328
1329 double HcalIsoTrkAnalyzer::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
1330 return reco::deltaR(vec1.eta(), vec1.phi(), vec2.eta(), vec2.phi());
1331 }
1332
1333 double HcalIsoTrkAnalyzer::trackP(const reco::Track* pTrack,
1334 const edm::Handle<reco::GenParticleCollection>& genParticles) {
1335 double pmom = -1.0;
1336 if (genParticles.isValid()) {
1337 double mindR(999.9);
1338 for (const auto& p : (*genParticles)) {
1339 double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), p.momentum().Eta(), p.momentum().Phi());
1340 if (dR < mindR) {
1341 mindR = dR;
1342 pmom = p.momentum().R();
1343 }
1344 }
1345 }
1346 return pmom;
1347 }
1348
1349 double HcalIsoTrkAnalyzer::rhoh(const edm::Handle<CaloTowerCollection>& tower) {
1350 std::vector<double> sumPFNallSMDQH2;
1351 sumPFNallSMDQH2.reserve(phibins_.size() * etabins_.size());
1352
1353 for (auto eta : etabins_) {
1354 for (auto phi : phibins_) {
1355 double hadder = 0;
1356 for (const auto& pf_it : (*tower)) {
1357 if (fabs(eta - pf_it.eta()) > etahalfdist_)
1358 continue;
1359 if (fabs(reco::deltaPhi(phi, pf_it.phi())) > phihalfdist_)
1360 continue;
1361 hadder += pf_it.hadEt();
1362 }
1363 sumPFNallSMDQH2.emplace_back(hadder);
1364 }
1365 }
1366
1367 double evt_smdq(0);
1368 std::sort(sumPFNallSMDQH2.begin(), sumPFNallSMDQH2.end());
1369 if (sumPFNallSMDQH2.size() % 2)
1370 evt_smdq = sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 1) / 2];
1371 else
1372 evt_smdq = (sumPFNallSMDQH2[sumPFNallSMDQH2.size() / 2] + sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 2) / 2]) / 2.;
1373 double rhoh = evt_smdq / (etadist_ * phidist_);
1374 #ifdef EDM_ML_DEBUG
1375 if (debug_)
1376 edm::LogVerbatim("HcalIsoTrack") << "Rho " << evt_smdq << ":" << rhoh;
1377 #endif
1378 return rhoh;
1379 }
1380
1381 double HcalIsoTrkAnalyzer::eThreshold(const DetId& id, const CaloGeometry* geo) const {
1382 double eThr(hitEthrEB_);
1383 if (usePFThresh_) {
1384 eThr = static_cast<double>((*eThresholds_)[id]);
1385 } else {
1386 const GlobalPoint& pos = geo->getPosition(id);
1387 double eta = std::abs(pos.eta());
1388 if (id.subdetId() != EcalBarrel) {
1389 eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
1390 if (eThr < hitEthrEELo_)
1391 eThr = hitEthrEELo_;
1392 else if (eThr > hitEthrEEHi_)
1393 eThr = hitEthrEEHi_;
1394 }
1395 }
1396 return eThr;
1397 }
1398
1399 DetId HcalIsoTrkAnalyzer::newId(const DetId& id) {
1400 HcalDetId hid(id);
1401 if (hep17_ && ((hid.iphi() < 63) || (hid.iphi() > 66) || (hid.zside() < 0)))
1402 return id;
1403 for (unsigned int k = 0; k < oldID_.size(); ++k) {
1404 if ((hid.subdetId() == oldDet_[k]) && (hid.ietaAbs() == oldEta_[k]) && (hid.depth() == oldDepth_[k])) {
1405 return static_cast<DetId>(HcalDetId(hid.subdet(), hid.ieta(), hid.iphi(), newDepth_[k]));
1406 }
1407 }
1408 return id;
1409 }
1410
1411 void HcalIsoTrkAnalyzer::storeEnergy(int indx,
1412 const HcalRespCorrs* respCorrs,
1413 const std::vector<DetId>& ids,
1414 std::vector<double>& edet,
1415 double& eHcal,
1416 std::vector<unsigned int>* detIds,
1417 std::vector<double>* hitEnergies) {
1418 double ehcal(0);
1419 if (unCorrect_) {
1420 for (unsigned int k = 0; k < ids.size(); ++k) {
1421 double corr = (respCorrs->getValues(ids[k]))->getValue();
1422 if (corr != 0)
1423 edet[k] /= corr;
1424 ehcal += edet[k];
1425 }
1426 } else {
1427 for (const auto& en : edet)
1428 ehcal += en;
1429 }
1430 if ((std::abs(ehcal - eHcal) > 0.001) && (!unCorrect_))
1431 edm::LogWarning("HcalIsoTrack") << "Check inconsistent energies: " << indx << " " << eHcal << ":" << ehcal
1432 << " from " << ids.size() << " cells";
1433 eHcal = hcalScale_ * ehcal;
1434
1435 if (collapseDepth_) {
1436 std::map<HcalDetId, double> hitMap;
1437 for (unsigned int k = 0; k < ids.size(); ++k) {
1438 HcalDetId id = hdc_->mergedDepthDetId(HcalDetId(ids[k]));
1439 auto itr = hitMap.find(id);
1440 if (itr == hitMap.end()) {
1441 hitMap[id] = edet[k];
1442 } else {
1443 (itr->second) += edet[k];
1444 }
1445 }
1446 detIds->reserve(hitMap.size());
1447 hitEnergies->reserve(hitMap.size());
1448 for (const auto& hit : hitMap) {
1449 detIds->emplace_back(hit.first.rawId());
1450 hitEnergies->emplace_back(hit.second);
1451 }
1452 } else {
1453 detIds->reserve(ids.size());
1454 hitEnergies->reserve(ids.size());
1455 for (unsigned int k = 0; k < ids.size(); ++k) {
1456 detIds->emplace_back(ids[k].rawId());
1457 hitEnergies->emplace_back(edet[k]);
1458 }
1459 }
1460 #ifdef EDM_ML_DEBUG
1461 if (debug_) {
1462 edm::LogVerbatim("HcalIsoTrack") << "Input to storeEnergy with " << ids.size() << " cells";
1463 for (unsigned int k = 0; k < ids.size(); ++k)
1464 edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId(ids[k]) << " E " << edet[k];
1465 edm::LogVerbatim("HcalIsoTrack") << "Output of storeEnergy with " << detIds->size() << " cells and Etot " << eHcal;
1466 for (unsigned int k = 0; k < detIds->size(); ++k)
1467 edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId((*detIds)[k]) << " E " << (*hitEnergies)[k];
1468 }
1469 #endif
1470 }
1471
1472 bool HcalIsoTrkAnalyzer::notaMuon(const reco::Track* pTrack0, const edm::Handle<reco::MuonCollection>& muonh) {
1473 bool flag(true);
1474 for (reco::MuonCollection::const_iterator recMuon = muonh->begin(); recMuon != muonh->end(); ++recMuon) {
1475 if (recMuon->innerTrack().isNonnull()) {
1476 const reco::Track* pTrack = (recMuon->innerTrack()).get();
1477 bool mediumMuon = (((recMuon->isPFMuon()) && (recMuon->isGlobalMuon() || recMuon->isTrackerMuon())) &&
1478 (recMuon->innerTrack()->validFraction() > 0.49));
1479 if (mediumMuon) {
1480 double chiGlobal = ((recMuon->globalTrack().isNonnull()) ? recMuon->globalTrack()->normalizedChi2() : 999);
1481 bool goodGlob = (recMuon->isGlobalMuon() && chiGlobal < 3 &&
1482 recMuon->combinedQuality().chi2LocalPosition < 12 && recMuon->combinedQuality().trkKink < 20);
1483 mediumMuon = muon::segmentCompatibility(*recMuon) > (goodGlob ? 0.303 : 0.451);
1484 }
1485 if (mediumMuon) {
1486 double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), pTrack0->eta(), pTrack0->phi());
1487 if (dR < 0.1) {
1488 flag = false;
1489 break;
1490 }
1491 }
1492 }
1493 }
1494 return flag;
1495 }
1496
1497
1498 DEFINE_FWK_MODULE(HcalIsoTrkAnalyzer);
