Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:59:02

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