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