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