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