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