File indexing completed on 2023-10-25 09:34:52
0001 #include <memory>
0002 #include <iostream>
0003 #include <vector>
0004
0005 #include <TTree.h>
0006
0007
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "FWCore/Common/interface/TriggerNames.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0018 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0019 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0020 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0021 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0022 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0023 #include "DataFormats/MuonReco/interface/Muon.h"
0024 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0025 #include "DataFormats/VertexReco/interface/Vertex.h"
0026 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0027 #include "DataFormats/TrackReco/interface/Track.h"
0028 #include "DataFormats/Common/interface/TriggerResults.h"
0029 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0030 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0031 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0032 #include "HLTrigger/HLTcore/interface/HLTConfigData.h"
0033
0034 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0035 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0036 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0037
0038 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0039 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0040 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0041
0042 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0043 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0044 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0045 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0046 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0047 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0048 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0049 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0050 #include "MagneticField/Engine/interface/MagneticField.h"
0051 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0052 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0053
0054 class HcalRaddamMuon : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0055 public:
0056 explicit HcalRaddamMuon(const edm::ParameterSet&);
0057 ~HcalRaddamMuon() override;
0058
0059 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062 void beginJob() override;
0063 void analyze(const edm::Event&, const edm::EventSetup&) override;
0064 void endJob() override;
0065 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0066 void endRun(edm::Run const&, edm::EventSetup const&) override;
0067 void clearVectors();
0068 int matchId(const HcalDetId&, const HcalDetId&);
0069 double activeLength(const DetId&);
0070
0071
0072 std::vector<double> PtGlob, track_cosmic_positionIX, track_cosmic_positionIY, track_cosmic_positionIZ,
0073 track_cosmic_positionOX, track_cosmic_positionOY, track_cosmic_positionOZ;
0074
0075 std::vector<double> track_cosmic_momentumOX, track_cosmic_momentumOY, track_cosmic_momentumOZ,
0076 track_cosmic_momentumIX, track_cosmic_momentumIY, track_cosmic_momentumIZ, track_cosmic_detIDinner,
0077 track_cosmic_detIDouter;
0078 std::vector<double> EtaGlob;
0079 std::vector<double> PhiGlob, chiGlobal, GlobalMuonHits, MatchedStat, GlobalTrckPt, GlobalTrckEta, GlobalTrckPhi;
0080 std::vector<double> TrackerLayer, innerTrackpt, innerTracketa, innerTrackphi;
0081 std::vector<double> NumPixelLayers, chiTracker, DxyTracker, DzTracker;
0082 std::vector<double> OuterTrackPt, OuterTrackEta, OuterTrackPhi, OuterTrackChi, OuterTrackHits, OuterTrackRHits;
0083 std::vector<double> trackerlayer_hits, No_pixelLayers, NormChi2, ImpactParameter;
0084 std::vector<double> Tight_GlobalMuonTrkFit, Tight_MuonHits, Tight_MatchedStations, Tight_LongPara, Tight_PixelHits,
0085 Tight_TrkerLayers, Tight_TransImpara, High_TrackLayers;
0086 std::vector<bool> innerTrack, OuterTrack, GlobalTrack;
0087 std::vector<double> IsolationR04, IsolationR03;
0088 std::vector<double> Energy, MuonHcalEnergy, MuonEcalEnergy, MuonHOEnergy, MuonEcal3x3Energy, MuonHcal1x1Energy, Pmuon;
0089 std::vector<unsigned int> MuonEcalDetId, MuonHcalDetId, MuonEHcalDetId, MuonHcalHot, MuonHcalHotCalo;
0090 std::vector<double> MuonHcalDepth1Energy, MuonHcalDepth2Energy, MuonHcalDepth3Energy, MuonHcalDepth4Energy,
0091 MuonHcalDepth5Energy, MuonHcalDepth6Energy, MuonHcalDepth7Energy;
0092 std::vector<double> MuonHcalDepth1HotEnergy, MuonHcalDepth2HotEnergy, MuonHcalDepth3HotEnergy,
0093 MuonHcalDepth4HotEnergy, MuonHcalDepth5HotEnergy, MuonHcalDepth6HotEnergy, MuonHcalDepth7HotEnergy;
0094
0095
0096 std::vector<double> MuonHcalDepth1EnergyCalo, MuonHcalDepth2EnergyCalo, MuonHcalDepth3EnergyCalo,
0097 MuonHcalDepth4EnergyCalo, MuonHcalDepth5EnergyCalo, MuonHcalDepth6EnergyCalo, MuonHcalDepth7EnergyCalo;
0098 std::vector<double> MuonHcalDepth1HotEnergyCalo, MuonHcalDepth2HotEnergyCalo, MuonHcalDepth3HotEnergyCalo,
0099 MuonHcalDepth4HotEnergyCalo, MuonHcalDepth5HotEnergyCalo, MuonHcalDepth6HotEnergyCalo,
0100 MuonHcalDepth7HotEnergyCalo;
0101
0102
0103 std::vector<double> MuonHcalActiveLength;
0104 std::vector<HcalDDDRecConstants::HcalActiveLength> actHB, actHE;
0105 int maxDepth_, type;
0106 edm::Service<TFileService> fs;
0107
0108 HLTConfigProvider hltConfig_;
0109 const edm::InputTag hlTriggerResults_, muonsrc_;
0110 const int verbosity_, useRaw_;
0111 const bool isAOD_;
0112 std::string hltlabel_;
0113 std::vector<std::string> all_triggers, all_triggers1, all_triggers2, all_triggers3, all_triggers4, all_triggers5;
0114
0115
0116 std::vector<bool> isHB, isHE;
0117 TTree* TREE;
0118 std::vector<bool> all_ifTriggerpassed;
0119 std::vector<bool> muon_is_good, muon_global, muon_tracker, Trk_match_MuStat;
0120 std::vector<int> hltresults;
0121 std::vector<float> energy_hb, time_hb;
0122 std::vector<std::string> hltpaths, TrigName_;
0123 std::vector<int> v_RH_h3x3_ieta;
0124 std::vector<int> v_RH_h3x3_iphi;
0125 std::vector<double> v_RH_h3x3_ene, PxGlob, PyGlob, PzGlob, Pthetha;
0126 std::vector<double> PCharge, PChi2, PD0, PD0Error, dxyWithBS, dzWithBS, PdxyTrack, PdzTrack, PNormalizedChi2, PNDoF,
0127 PValidHits, PLostHits, NPvx, NPvy, NPvz, NQOverP, NQOverPError, NTrkMomentum, NRefPointX, NRefPointY, NRefPointZ;
0128 double h3x3, h3x3Calo;
0129 unsigned int RunNumber, EventNumber, LumiNumber, BXNumber;
0130 double _RecoMuon1TrackIsoSumPtMaxCutValue_03, _RecoMuon1TrackIsoSumPtMaxCutValue_04;
0131 int ntriggers;
0132 std::vector<double> track_cosmic_xposition, track_cosmic_yposition, track_cosmic_zposition, track_cosmic_xmomentum,
0133 track_cosmic_ymomentum, track_cosmic_zmomentum, track_cosmic_rad, track_cosmic_detid;
0134
0135 edm::EDGetTokenT<edm::PCaloHitContainer> tok_hcal_;
0136 edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0137 edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0138 edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0139 edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0140 edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0141 edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0142 edm::EDGetTokenT<reco::MuonCollection> tok_muon_;
0143
0144 edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0145 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0146 edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0147 edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_topo_;
0148 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0149 edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0150 edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0151 };
0152
0153 HcalRaddamMuon::HcalRaddamMuon(const edm::ParameterSet& iConfig)
0154 : hlTriggerResults_(iConfig.getUntrackedParameter<edm::InputTag>("hlTriggerResults_")),
0155 muonsrc_(iConfig.getUntrackedParameter<edm::InputTag>("muonSource")),
0156 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0157 useRaw_(iConfig.getUntrackedParameter<int>("UseRaw", 0)),
0158 isAOD_(iConfig.getUntrackedParameter<bool>("IsAOD", false)) {
0159 usesResource(TFileService::kSharedResource);
0160
0161
0162 maxDepth_ = iConfig.getUntrackedParameter<int>("MaxDepth", 4);
0163 if (maxDepth_ > 7)
0164 maxDepth_ = 7;
0165 else if (maxDepth_ < 1)
0166 maxDepth_ = 4;
0167
0168 tok_hcal_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
0169 tok_trigRes_ = consumes<edm::TriggerResults>(hlTriggerResults_);
0170 tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0171 tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0172 if (isAOD_) {
0173 tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEB"));
0174 tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEE"));
0175 tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("reducedHcalRecHits", "hbhereco"));
0176 } else {
0177 tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0178 tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0179 tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0180 }
0181 tok_muon_ = consumes<reco::MuonCollection>(muonsrc_);
0182
0183 tok_ddrec_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>();
0184 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0185 tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0186 tok_topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0187 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0188 tok_ecalChStatus_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0189 tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0190 }
0191
0192 HcalRaddamMuon::~HcalRaddamMuon() {
0193
0194
0195 }
0196
0197
0198
0199
0200
0201
0202 void HcalRaddamMuon::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0203 clearVectors();
0204 RunNumber = iEvent.id().run();
0205 EventNumber = iEvent.id().event();
0206 LumiNumber = iEvent.id().luminosityBlock();
0207 BXNumber = iEvent.bunchCrossing();
0208
0209 edm::Handle<edm::PCaloHitContainer> calosimhits;
0210 iEvent.getByToken(tok_hcal_, calosimhits);
0211
0212 edm::Handle<edm::TriggerResults> _Triggers;
0213 iEvent.getByToken(tok_trigRes_, _Triggers);
0214
0215 if ((verbosity_ % 10) > 1)
0216 edm::LogVerbatim("HBHEMuon") << "size of all triggers " << all_triggers.size();
0217 int Ntriggers = all_triggers.size();
0218
0219 if ((verbosity_ % 10) > 1)
0220 edm::LogVerbatim("HBHEMuon") << "size of HLT MENU: " << _Triggers->size();
0221
0222 if (_Triggers.isValid()) {
0223 const edm::TriggerNames& triggerNames_ = iEvent.triggerNames(*_Triggers);
0224
0225 std::vector<int> index;
0226 for (int i = 0; i < Ntriggers; i++) {
0227 index.push_back(triggerNames_.triggerIndex(all_triggers[i]));
0228 int triggerSize = int(_Triggers->size());
0229 if ((verbosity_ % 10) > 2)
0230 edm::LogVerbatim("HBHEMuon") << "outside loop " << index[i] << "\ntriggerSize " << triggerSize;
0231 if (index[i] < triggerSize) {
0232 hltresults.push_back(_Triggers->accept(index[i]));
0233 if ((verbosity_ % 10) > 2)
0234 edm::LogVerbatim("HBHEMuon") << "trigger_info " << triggerSize << " triggerSize " << index[i]
0235 << " trigger_index " << hltresults.at(i) << " hltresult ";
0236 } else {
0237 edm::LogVerbatim("HBHEMuon") << "Requested HLT path \""
0238 << "\" does not exist";
0239 }
0240 }
0241 }
0242
0243
0244 const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0245 const MagneticField* bField = &iSetup.getData(tok_magField_);
0246 const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0247 const EcalSeverityLevelAlgo* sevlv = &iSetup.getData(tok_sevlv_);
0248 const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0249 const HcalTopology* theHBHETopology = &iSetup.getData(tok_topo_);
0250
0251 edm::Handle<reco::BeamSpot> bmspot;
0252 iEvent.getByToken(tok_bs_, bmspot);
0253
0254 edm::Handle<reco::VertexCollection> vtx;
0255 iEvent.getByToken(tok_recVtx_, vtx);
0256
0257 edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0258 edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0259 iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
0260 iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
0261
0262 edm::Handle<HBHERecHitCollection> hbhe;
0263 iEvent.getByToken(tok_hbhe_, hbhe);
0264
0265 edm::Handle<reco::MuonCollection> _Muon;
0266 iEvent.getByToken(tok_muon_, _Muon);
0267 const reco::Vertex& vertex = (*(vtx)->begin());
0268
0269 math::XYZPoint bspot;
0270 bspot = (bmspot.isValid()) ? bmspot->position() : math::XYZPoint(0, 0, 0);
0271
0272 if (_Muon.isValid()) {
0273 for (reco::MuonCollection::const_iterator RecMuon = _Muon->begin(); RecMuon != _Muon->end(); ++RecMuon) {
0274 muon_is_good.push_back(RecMuon->isPFMuon());
0275 muon_global.push_back(RecMuon->isGlobalMuon());
0276 muon_tracker.push_back(RecMuon->isTrackerMuon());
0277 PtGlob.push_back((RecMuon)->pt());
0278 EtaGlob.push_back(RecMuon->eta());
0279 PhiGlob.push_back(RecMuon->phi());
0280 Energy.push_back(RecMuon->energy());
0281 Pmuon.push_back(RecMuon->p());
0282
0283
0284 if (RecMuon->track().isNonnull()) {
0285 TrackerLayer.push_back(RecMuon->track()->hitPattern().trackerLayersWithMeasurement());
0286 } else {
0287 TrackerLayer.push_back(-1);
0288 }
0289 if (RecMuon->innerTrack().isNonnull()) {
0290 innerTrack.push_back(true);
0291 NumPixelLayers.push_back(RecMuon->innerTrack()->hitPattern().pixelLayersWithMeasurement());
0292 chiTracker.push_back(RecMuon->innerTrack()->normalizedChi2());
0293 DxyTracker.push_back(fabs(RecMuon->innerTrack()->dxy((vertex).position())));
0294 DzTracker.push_back(fabs(RecMuon->innerTrack()->dz((vertex).position())));
0295 innerTrackpt.push_back(RecMuon->innerTrack()->pt());
0296 innerTracketa.push_back(RecMuon->innerTrack()->eta());
0297 innerTrackphi.push_back(RecMuon->innerTrack()->phi());
0298 Tight_PixelHits.push_back(RecMuon->innerTrack()->hitPattern().numberOfValidPixelHits());
0299 } else {
0300 innerTrack.push_back(false);
0301 NumPixelLayers.push_back(0);
0302 chiTracker.push_back(0);
0303 DxyTracker.push_back(0);
0304 DzTracker.push_back(0);
0305 innerTrackpt.push_back(0);
0306 innerTracketa.push_back(0);
0307 innerTrackphi.push_back(0);
0308 Tight_PixelHits.push_back(0);
0309 }
0310
0311
0312 if (RecMuon->outerTrack().isNonnull()) {
0313 OuterTrack.push_back(true);
0314 OuterTrackPt.push_back(RecMuon->outerTrack()->pt());
0315 OuterTrackEta.push_back(RecMuon->outerTrack()->eta());
0316 OuterTrackPhi.push_back(RecMuon->outerTrack()->phi());
0317 OuterTrackChi.push_back(RecMuon->outerTrack()->normalizedChi2());
0318 OuterTrackHits.push_back(RecMuon->outerTrack()->numberOfValidHits());
0319 OuterTrackRHits.push_back(RecMuon->outerTrack()->recHitsSize());
0320 } else {
0321 OuterTrack.push_back(false);
0322 OuterTrackPt.push_back(0);
0323 OuterTrackEta.push_back(0);
0324 OuterTrackPhi.push_back(0);
0325 OuterTrackChi.push_back(0);
0326 OuterTrackHits.push_back(0);
0327 OuterTrackRHits.push_back(0);
0328 }
0329
0330 if (RecMuon->globalTrack().isNonnull()) {
0331 GlobalTrack.push_back(true);
0332 chiGlobal.push_back(RecMuon->globalTrack()->normalizedChi2());
0333 GlobalMuonHits.push_back(RecMuon->globalTrack()->hitPattern().numberOfValidMuonHits());
0334 MatchedStat.push_back(RecMuon->numberOfMatchedStations());
0335 GlobalTrckPt.push_back(RecMuon->globalTrack()->pt());
0336 GlobalTrckEta.push_back(RecMuon->globalTrack()->eta());
0337 GlobalTrckPhi.push_back(RecMuon->globalTrack()->phi());
0338 Tight_TransImpara.push_back(fabs(RecMuon->muonBestTrack()->dxy(vertex.position())));
0339 Tight_LongPara.push_back(fabs(RecMuon->muonBestTrack()->dz(vertex.position())));
0340 } else {
0341 GlobalTrack.push_back(false);
0342 chiGlobal.push_back(0);
0343 GlobalMuonHits.push_back(0);
0344 MatchedStat.push_back(0);
0345 GlobalTrckPt.push_back(0);
0346 GlobalTrckEta.push_back(0);
0347 GlobalTrckPhi.push_back(0);
0348 Tight_TransImpara.push_back(0);
0349 Tight_LongPara.push_back(0);
0350 }
0351
0352 IsolationR04.push_back(
0353 ((RecMuon->pfIsolationR04().sumChargedHadronPt +
0354 std::max(0.,
0355 RecMuon->pfIsolationR04().sumNeutralHadronEt + RecMuon->pfIsolationR04().sumPhotonEt -
0356 (0.5 * RecMuon->pfIsolationR04().sumPUPt))) /
0357 RecMuon->pt()));
0358
0359 IsolationR03.push_back(
0360 ((RecMuon->pfIsolationR03().sumChargedHadronPt +
0361 std::max(0.,
0362 RecMuon->pfIsolationR03().sumNeutralHadronEt + RecMuon->pfIsolationR03().sumPhotonEt -
0363 (0.5 * RecMuon->pfIsolationR03().sumPUPt))) /
0364 RecMuon->pt()));
0365
0366 MuonEcalEnergy.push_back(RecMuon->calEnergy().emS9);
0367 MuonHcalEnergy.push_back(RecMuon->calEnergy().hadS9);
0368 MuonHOEnergy.push_back(RecMuon->calEnergy().hoS9);
0369
0370 double eEcal(0), eHcal(0), activeL(0), eHcalDepth[7], eHcalDepthHot[7], eHcalDepthCalo[7], eHcalDepthHotCalo[7];
0371 unsigned int isHot = 0;
0372 unsigned int isHotCalo = 0;
0373
0374 for (int i = 0; i < 7; ++i)
0375 eHcalDepth[i] = eHcalDepthHot[i] = eHcalDepthCalo[i] = eHcalDepthHotCalo[i] = -10000;
0376
0377 if (RecMuon->innerTrack().isNonnull()) {
0378 const reco::Track* pTrack = (RecMuon->innerTrack()).get();
0379 spr::propagatedTrackID trackID = spr::propagateCALO(pTrack, geo, bField, (((verbosity_ / 100) % 10 > 0)));
0380
0381 MuonEcalDetId.push_back((trackID.detIdECAL)());
0382 MuonHcalDetId.push_back((trackID.detIdHCAL)());
0383 MuonEHcalDetId.push_back((trackID.detIdEHCAL)());
0384
0385 if (trackID.okECAL) {
0386 const DetId isoCell(trackID.detIdECAL);
0387 std::pair<double, bool> e3x3 = spr::eECALmatrix(isoCell,
0388 barrelRecHitsHandle,
0389 endcapRecHitsHandle,
0390 *theEcalChStatus,
0391 geo,
0392 caloTopology,
0393 sevlv,
0394 1,
0395 1,
0396 -100.0,
0397 -100.0,
0398 -500.0,
0399 500.0,
0400 false);
0401
0402 eEcal = e3x3.first;
0403 if (((verbosity_ / 10) % 10) > 1)
0404 edm::LogVerbatim("HBHEMuon") << "eEcal" << eEcal;
0405 }
0406
0407 if (trackID.okHCAL) {
0408 const DetId closestCell(trackID.detIdHCAL);
0409 eHcal = spr::eHCALmatrix(theHBHETopology,
0410 closestCell,
0411 hbhe,
0412 0,
0413 0,
0414 false,
0415 true,
0416 -100.0,
0417 -100.0,
0418 -100.0,
0419 -100.0,
0420 -500.,
0421 500.,
0422 useRaw_);
0423
0424 int iphi = ((HcalDetId)(closestCell)).iphi();
0425 int zside = ((HcalDetId)(closestCell)).iphi();
0426 int depthHE = theHBHETopology->dddConstants()->getMinDepth(1, 16, iphi, zside);
0427 if (((verbosity_ / 10) % 10) > 1)
0428 edm::LogVerbatim("HBHEMuon") << "eHcal " << eHcal;
0429 std::vector<std::pair<double, int> > ehdepth;
0430 spr::energyHCALCell((HcalDetId)closestCell,
0431 hbhe,
0432 ehdepth,
0433 maxDepth_,
0434 -100.0,
0435 -100.0,
0436 -100.0,
0437 -100.0,
0438 -500.0,
0439 500.0,
0440 useRaw_,
0441 depthHE,
0442 (((verbosity_ / 1000) % 10) > 0));
0443 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0444 eHcalDepth[ehdepth[i].second - 1] = ehdepth[i].first;
0445 if (((verbosity_ / 10) % 10) > 1)
0446 edm::LogVerbatim("HBHEMuon")
0447 << "eHcalDepth " << i << ":" << (ehdepth[i].second - 1) << ":" << eHcalDepth[ehdepth[i].second - 1];
0448 }
0449
0450 eHcal = spr::eHCALmatrix(theHBHETopology,
0451 closestCell,
0452 calosimhits,
0453 0,
0454 0,
0455 false,
0456 true,
0457 -100.0,
0458 -100.0,
0459 -100.0,
0460 -100.0,
0461 -500.,
0462 500.,
0463 useRaw_);
0464
0465 if (((verbosity_ / 10) % 10) > 1)
0466 edm::LogVerbatim("HBHEMuon") << "eHcal " << eHcal;
0467 const DetId closestCellCalo(trackID.detIdHCAL);
0468 iphi = ((HcalDetId)(closestCellCalo)).iphi();
0469 zside = ((HcalDetId)(closestCellCalo)).iphi();
0470 depthHE = theHBHETopology->dddConstants()->getMinDepth(1, 16, iphi, zside);
0471 std::vector<std::pair<double, int> > ehdepthCalo;
0472 spr::energyHCALCell((HcalDetId)closestCellCalo,
0473 calosimhits,
0474 ehdepthCalo,
0475 maxDepth_,
0476 -100.0,
0477 -100.0,
0478 -100.0,
0479 -100.0,
0480 -500.0,
0481 500.0,
0482 useRaw_,
0483 depthHE,
0484 (((verbosity_ / 1000) % 10) > 0));
0485 for (unsigned int i = 0; i < ehdepthCalo.size(); ++i) {
0486 eHcalDepthCalo[ehdepthCalo[i].second - 1] = ehdepthCalo[i].first;
0487 if (((verbosity_ / 10) % 10) > 1)
0488 edm::LogVerbatim("HBHEMuon") << "eHcalDepthCalo " << i << ":" << (ehdepth[i].second - 1) << ":"
0489 << eHcalDepth[ehdepth[i].second - 1];
0490 }
0491
0492 HcalDetId hcid0(closestCell.rawId());
0493 activeL = activeLength(trackID.detIdHCAL);
0494
0495 if (((verbosity_ / 10) % 10) > 0)
0496 edm::LogVerbatim("HBHEMuon") << "activeL " << activeL;
0497 HcalDetId hotCell, hotCellCalo;
0498 h3x3 = spr::eHCALmatrix(geo, theHBHETopology, closestCell, hbhe, 1, 1, hotCell, false, useRaw_, false);
0499 h3x3Calo = spr::eHCALmatrix(
0500 geo, theHBHETopology, closestCellCalo, calosimhits, 1, 1, hotCellCalo, false, useRaw_, false);
0501
0502 isHot = matchId(closestCell, hotCell);
0503 isHotCalo = matchId(closestCellCalo, hotCellCalo);
0504
0505 if (((verbosity_ / 10) % 10) > 1)
0506 edm::LogVerbatim("HBHEMuon") << "hcal 3X3 < " << h3x3 << ">"
0507 << " ClosestCell <" << (HcalDetId)(closestCell) << "> hotCell id < " << hotCell
0508 << "> isHot" << isHot;
0509 if (hotCell != HcalDetId()) {
0510 iphi = ((HcalDetId)(hotCell)).iphi();
0511 zside = ((HcalDetId)(hotCell)).iphi();
0512 depthHE = theHBHETopology->dddConstants()->getMinDepth(1, 16, iphi, zside);
0513 std::vector<std::pair<double, int> > ehdepth;
0514 spr::energyHCALCell(hotCell,
0515 hbhe,
0516 ehdepth,
0517 maxDepth_,
0518 -100.0,
0519 -100.0,
0520 -100.0,
0521 -100.0,
0522 -500.0,
0523 500.0,
0524 depthHE,
0525 false);
0526 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0527 eHcalDepthHot[ehdepth[i].second - 1] = ehdepth[i].first;
0528 if (((verbosity_ / 10) % 10) > 1)
0529 edm::LogVerbatim("HBHEMuon") << "eHcalDepthHot " << i << ":" << (ehdepth[i].second - 1) << ":"
0530 << eHcalDepthHot[ehdepth[i].second - 1];
0531 }
0532 }
0533
0534 if (hotCellCalo != HcalDetId()) {
0535 iphi = ((HcalDetId)(hotCellCalo)).iphi();
0536 zside = ((HcalDetId)(hotCellCalo)).iphi();
0537 depthHE = theHBHETopology->dddConstants()->getMinDepth(1, 16, iphi, zside);
0538 std::vector<std::pair<double, int> > ehdepthCalo;
0539
0540 spr::energyHCALCell(hotCellCalo,
0541 calosimhits,
0542 ehdepthCalo,
0543 maxDepth_,
0544 -100.0,
0545 -100.0,
0546 -100.0,
0547 -100.0,
0548 -500.0,
0549 500.0,
0550 useRaw_,
0551 depthHE,
0552 false);
0553 for (unsigned int i = 0; i < ehdepthCalo.size(); ++i) {
0554 eHcalDepthHotCalo[ehdepthCalo[i].second - 1] = ehdepthCalo[i].first;
0555 if (((verbosity_ / 10) % 10) > 1)
0556 edm::LogVerbatim("HBHEMuon") << "eHcalDepthHotCalo " << i << ":" << (ehdepth[i].second - 1) << ":"
0557 << eHcalDepthHot[ehdepth[i].second - 1];
0558 }
0559 }
0560 }
0561 } else {
0562 MuonEcalDetId.push_back(0);
0563 MuonHcalDetId.push_back(0);
0564 MuonEHcalDetId.push_back(0);
0565 }
0566
0567 MuonEcal3x3Energy.push_back(eEcal);
0568 MuonHcal1x1Energy.push_back(eHcal);
0569 MuonHcalDepth1Energy.push_back(eHcalDepth[0]);
0570 MuonHcalDepth2Energy.push_back(eHcalDepth[1]);
0571 MuonHcalDepth3Energy.push_back(eHcalDepth[2]);
0572 MuonHcalDepth4Energy.push_back(eHcalDepth[3]);
0573 MuonHcalDepth5Energy.push_back(eHcalDepth[4]);
0574 MuonHcalDepth6Energy.push_back(eHcalDepth[5]);
0575 MuonHcalDepth7Energy.push_back(eHcalDepth[6]);
0576 MuonHcalDepth1HotEnergy.push_back(eHcalDepthHot[0]);
0577 MuonHcalDepth2HotEnergy.push_back(eHcalDepthHot[1]);
0578 MuonHcalDepth3HotEnergy.push_back(eHcalDepthHot[2]);
0579 MuonHcalDepth4HotEnergy.push_back(eHcalDepthHot[3]);
0580 MuonHcalDepth5HotEnergy.push_back(eHcalDepthHot[4]);
0581 MuonHcalDepth6HotEnergy.push_back(eHcalDepthHot[5]);
0582 MuonHcalDepth7HotEnergy.push_back(eHcalDepthHot[6]);
0583 MuonHcalHot.push_back(isHot);
0584
0585
0586 MuonHcalDepth1EnergyCalo.push_back(eHcalDepthCalo[0]);
0587 MuonHcalDepth2EnergyCalo.push_back(eHcalDepthCalo[1]);
0588 MuonHcalDepth3EnergyCalo.push_back(eHcalDepthCalo[2]);
0589 MuonHcalDepth4EnergyCalo.push_back(eHcalDepthCalo[3]);
0590 MuonHcalDepth5EnergyCalo.push_back(eHcalDepthCalo[4]);
0591 MuonHcalDepth6EnergyCalo.push_back(eHcalDepthCalo[5]);
0592 MuonHcalDepth7EnergyCalo.push_back(eHcalDepthCalo[6]);
0593 MuonHcalDepth1HotEnergyCalo.push_back(eHcalDepthHotCalo[0]);
0594 MuonHcalDepth2HotEnergyCalo.push_back(eHcalDepthHotCalo[1]);
0595 MuonHcalDepth3HotEnergyCalo.push_back(eHcalDepthHotCalo[2]);
0596 MuonHcalDepth4HotEnergyCalo.push_back(eHcalDepthHotCalo[3]);
0597 MuonHcalDepth5HotEnergyCalo.push_back(eHcalDepthHotCalo[4]);
0598 MuonHcalDepth6HotEnergyCalo.push_back(eHcalDepthHotCalo[5]);
0599 MuonHcalDepth7HotEnergyCalo.push_back(eHcalDepthHotCalo[6]);
0600 MuonHcalHotCalo.push_back(isHotCalo);
0601
0602
0603 MuonHcalActiveLength.push_back(activeL);
0604 }
0605 }
0606 TREE->Fill();
0607 }
0608
0609
0610 void HcalRaddamMuon::beginJob() {
0611 TREE = fs->make<TTree>("TREE", "TREE");
0612 TREE->Branch("Event_No", &EventNumber);
0613 TREE->Branch("Run_No", &RunNumber);
0614 TREE->Branch("LumiNumber", &LumiNumber);
0615 TREE->Branch("BXNumber", &BXNumber);
0616 TREE->Branch("pt_of_muon", &PtGlob);
0617 TREE->Branch("eta_of_muon", &EtaGlob);
0618 TREE->Branch("phi_of_muon", &PhiGlob);
0619 TREE->Branch("energy_of_muon", &Energy);
0620 TREE->Branch("p_of_muon", &Pmuon);
0621 TREE->Branch("PF_Muon", &muon_is_good);
0622 TREE->Branch("Global_Muon", &muon_global);
0623 TREE->Branch("Tracker_muon", &muon_tracker);
0624
0625 TREE->Branch("hcal_3into3", &MuonHcalEnergy);
0626 TREE->Branch("hcal_1x1", &MuonHcal1x1Energy);
0627 TREE->Branch("hcal_detID", &MuonHcalDetId);
0628 TREE->Branch("hcal_edepth1", &MuonHcalDepth1Energy);
0629 TREE->Branch("hcal_edepth2", &MuonHcalDepth2Energy);
0630 TREE->Branch("hcal_edepth3", &MuonHcalDepth3Energy);
0631 TREE->Branch("hcal_edepth4", &MuonHcalDepth4Energy);
0632 TREE->Branch("hcal_edepthHot1", &MuonHcalDepth1HotEnergy);
0633 TREE->Branch("hcal_edepthHot2", &MuonHcalDepth2HotEnergy);
0634 TREE->Branch("hcal_edepthHot3", &MuonHcalDepth3HotEnergy);
0635 TREE->Branch("hcal_edepthHot4", &MuonHcalDepth4HotEnergy);
0636
0637 TREE->Branch("hcal_edepth1PSim", &MuonHcalDepth1EnergyCalo);
0638 TREE->Branch("hcal_edepth2PSim", &MuonHcalDepth2EnergyCalo);
0639 TREE->Branch("hcal_edepth3PSim", &MuonHcalDepth3EnergyCalo);
0640 TREE->Branch("hcal_edepth4PSim", &MuonHcalDepth4EnergyCalo);
0641 TREE->Branch("hcal_edepthHot1PSim", &MuonHcalDepth1HotEnergyCalo);
0642 TREE->Branch("hcal_edepthHot2PSim", &MuonHcalDepth2HotEnergyCalo);
0643 TREE->Branch("hcal_edepthHot3PSim", &MuonHcalDepth3HotEnergyCalo);
0644 TREE->Branch("hcal_edepthHot4PSim", &MuonHcalDepth4HotEnergyCalo);
0645
0646 if (maxDepth_ > 4) {
0647 TREE->Branch("hcal_edepth5PSim", &MuonHcalDepth5EnergyCalo);
0648 TREE->Branch("hcal_edepthHot5PSim", &MuonHcalDepth5HotEnergyCalo);
0649 if (maxDepth_ > 5) {
0650 TREE->Branch("hcal_edepth6PSim", &MuonHcalDepth6EnergyCalo);
0651 TREE->Branch("hcal_edepthHot6PSim", &MuonHcalDepth6HotEnergyCalo);
0652 if (maxDepth_ > 6) {
0653 TREE->Branch("hcal_edepth7PSim", &MuonHcalDepth7EnergyCalo);
0654 TREE->Branch("hcal_edepthHot7PSim", &MuonHcalDepth7HotEnergyCalo);
0655 }
0656 }
0657 }
0658
0659 TREE->Branch("TrackerLayer", &TrackerLayer);
0660 TREE->Branch("innerTrack", &innerTrack);
0661 TREE->Branch("innerTrackpt", &innerTrackpt);
0662 TREE->Branch("innerTracketa", &innerTracketa);
0663 TREE->Branch("innerTrackphi", &innerTrackphi);
0664 TREE->Branch("MatchedStat", &MatchedStat);
0665 TREE->Branch("GlobalTrckPt", &GlobalTrckPt);
0666 TREE->Branch("GlobalTrckEta", &GlobalTrckEta);
0667 TREE->Branch("GlobalTrckPhi", &GlobalTrckPhi);
0668 TREE->Branch("NumPixelLayers", &NumPixelLayers);
0669 TREE->Branch("chiTracker", &chiTracker);
0670 TREE->Branch("DxyTracker", &DxyTracker);
0671 TREE->Branch("DzTracker", &DzTracker);
0672 TREE->Branch("OuterTrack", &OuterTrack);
0673 TREE->Branch("OuterTrackPt", &OuterTrackPt);
0674 TREE->Branch("OuterTrackEta", &OuterTrackEta);
0675 TREE->Branch("OuterTrackPhi", &OuterTrackPhi);
0676 TREE->Branch("OuterTrackHits", &OuterTrackHits);
0677 TREE->Branch("OuterTrackRHits", &OuterTrackRHits);
0678 TREE->Branch("OuterTrackChi", &OuterTrackChi);
0679 TREE->Branch("GlobalTrack", &GlobalTrack);
0680 TREE->Branch("GlobTrack_Chi", &chiGlobal);
0681 TREE->Branch("Global_Muon_Hits", &GlobalMuonHits);
0682 TREE->Branch("MatchedStations", &MatchedStat);
0683 TREE->Branch("Global_Track_Pt", &GlobalTrckPt);
0684 TREE->Branch("Global_Track_Eta", &GlobalTrckEta);
0685 TREE->Branch("Global_Track_Phi", &GlobalTrckPhi);
0686
0687 TREE->Branch("Tight_LongitudinalImpactparameter", &Tight_LongPara);
0688 TREE->Branch("Tight_TransImpactparameter", &Tight_TransImpara);
0689 TREE->Branch("InnerTrackPixelHits", &Tight_PixelHits);
0690 TREE->Branch("IsolationR04", &IsolationR04);
0691 TREE->Branch("IsolationR03", &IsolationR03);
0692
0693 TREE->Branch("hcal_cellHot", &MuonHcalHot);
0694 TREE->Branch("hcal_cellHotPSim", &MuonHcalHotCalo);
0695
0696 TREE->Branch("ecal_3into3", &MuonEcalEnergy);
0697 TREE->Branch("ecal_3x3", &MuonEcal3x3Energy);
0698 TREE->Branch("ecal_detID", &MuonEcalDetId);
0699 TREE->Branch("ehcal_detID", &MuonEHcalDetId);
0700 TREE->Branch("tracker_3into3", &MuonHOEnergy);
0701 TREE->Branch("activeLength", &MuonHcalActiveLength);
0702
0703
0704 TREE->Branch("hltresults", &hltresults);
0705 TREE->Branch("all_triggers", &all_triggers);
0706 TREE->Branch("rechit_energy", &energy_hb);
0707 TREE->Branch("rechit_time", &time_hb);
0708 }
0709
0710
0711 void HcalRaddamMuon::endJob() {}
0712
0713
0714 void HcalRaddamMuon::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0715 const HcalDDDRecConstants* hdc = &iSetup.getData(tok_ddrec_);
0716 actHB.clear();
0717 actHE.clear();
0718 actHB = hdc->getThickActive(0);
0719 actHE = hdc->getThickActive(1);
0720
0721 bool changed = true;
0722 all_triggers.clear();
0723 if (hltConfig_.init(iRun, iSetup, "HLT", changed)) {
0724
0725 edm::LogVerbatim("HBHEMuon") << "HLT config with process name 'HLT' successfully extracted";
0726 std::string string_search[5] = {"HLT_IsoMu_", "HLT_L1SingleMu_", "HLT_L2Mu", "HLT_Mu", "HLT_RelIso1p0Mu"};
0727 unsigned int ntriggers = hltConfig_.size();
0728 for (unsigned int t = 0; t < ntriggers; ++t) {
0729 std::string hltname(hltConfig_.triggerName(t));
0730 for (unsigned int ik = 0; ik < 5; ++ik) {
0731 if (hltname.find(string_search[ik]) != std::string::npos) {
0732 all_triggers.push_back(hltname);
0733 break;
0734 }
0735 }
0736 }
0737 edm::LogVerbatim("HBHEMuon") << "All triggers size in begin run" << all_triggers.size();
0738 } else {
0739 edm::LogError("HBHEMuon") << "Error! HLT config extraction with process name 'HLT' failed";
0740 }
0741
0742 }
0743
0744
0745 void HcalRaddamMuon::endRun(edm::Run const&, edm::EventSetup const&) {}
0746
0747 void HcalRaddamMuon::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0748 edm::ParameterSetDescription desc;
0749 desc.addUntracked<edm::InputTag>("hlTriggerResults", edm::InputTag("TriggerResults", "", "HLT"));
0750 desc.addUntracked<edm::InputTag>("muonSource", edm::InputTag("muons"));
0751 desc.addUntracked<int>("verbosity", 0);
0752 desc.addUntracked<int>("useRaw", 0);
0753 desc.add<bool>("isAOD", false);
0754 desc.addUntracked<int>("maxDepth", 4);
0755 descriptions.add("hcalRaddamMuon", desc);
0756 }
0757
0758 void HcalRaddamMuon::clearVectors() {
0759
0760 EventNumber = -99999;
0761 RunNumber = -99999;
0762 LumiNumber = -99999;
0763 BXNumber = -99999;
0764 energy_hb.clear();
0765 time_hb.clear();
0766 muon_is_good.clear();
0767 muon_global.clear();
0768 muon_tracker.clear();
0769 PtGlob.clear();
0770 EtaGlob.clear();
0771 PhiGlob.clear();
0772 Energy.clear();
0773 Pmuon.clear();
0774 TrackerLayer.clear();
0775 innerTrack.clear();
0776 NumPixelLayers.clear();
0777 chiTracker.clear();
0778 DxyTracker.clear();
0779 DzTracker.clear();
0780 innerTrackpt.clear();
0781 innerTracketa.clear();
0782 innerTrackphi.clear();
0783 Tight_PixelHits.clear();
0784 OuterTrack.clear();
0785 OuterTrackPt.clear();
0786 OuterTrackEta.clear();
0787 OuterTrackPhi.clear();
0788 OuterTrackHits.clear();
0789 OuterTrackRHits.clear();
0790 OuterTrackChi.clear();
0791 GlobalTrack.clear();
0792 chiGlobal.clear();
0793 GlobalMuonHits.clear();
0794 MatchedStat.clear();
0795 GlobalTrckPt.clear();
0796 GlobalTrckEta.clear();
0797 GlobalTrckPhi.clear();
0798 Tight_TransImpara.clear();
0799 Tight_LongPara.clear();
0800
0801 IsolationR04.clear();
0802 IsolationR03.clear();
0803 MuonEcalEnergy.clear();
0804 MuonHcalEnergy.clear();
0805 MuonHOEnergy.clear();
0806 MuonEcalDetId.clear();
0807 MuonHcalDetId.clear();
0808 MuonEHcalDetId.clear();
0809 MuonEcal3x3Energy.clear();
0810 MuonHcal1x1Energy.clear();
0811 MuonHcalDepth1Energy.clear();
0812 MuonHcalDepth2Energy.clear();
0813 MuonHcalDepth3Energy.clear();
0814 MuonHcalDepth4Energy.clear();
0815 MuonHcalDepth5Energy.clear();
0816 MuonHcalDepth6Energy.clear();
0817 MuonHcalDepth7Energy.clear();
0818
0819 MuonHcalDepth1HotEnergy.clear();
0820 MuonHcalDepth2HotEnergy.clear();
0821 MuonHcalDepth3HotEnergy.clear();
0822 MuonHcalDepth4HotEnergy.clear();
0823 MuonHcalDepth5HotEnergy.clear();
0824 MuonHcalDepth6HotEnergy.clear();
0825 MuonHcalDepth7HotEnergy.clear();
0826 MuonHcalHot.clear();
0827 MuonHcalActiveLength.clear();
0828 hltresults.clear();
0829 }
0830
0831 int HcalRaddamMuon::matchId(const HcalDetId& id1, const HcalDetId& id2) {
0832 HcalDetId kd1(id1.subdet(), id1.ieta(), id1.iphi(), 1);
0833 HcalDetId kd2(id2.subdet(), id2.ieta(), id2.iphi(), 1);
0834 int match = ((kd1 == kd2) ? 1 : 0);
0835 return match;
0836 }
0837
0838 double HcalRaddamMuon::activeLength(const DetId& id_) {
0839 HcalDetId id(id_);
0840 int ieta = id.ietaAbs();
0841 int depth = id.depth();
0842 double lx(0);
0843 if (id.subdet() == HcalBarrel) {
0844 if (((verbosity_ / 10) % 10) > 2)
0845 edm::LogVerbatim("HBHEMuon") << "actHB.size() " << actHB.size();
0846 for (unsigned int i = 0; i < actHB.size(); ++i) {
0847 if (ieta == actHB[i].ieta && depth == actHB[i].depth) {
0848 lx = actHB[i].thick;
0849 break;
0850 }
0851 }
0852 } else {
0853 if (((verbosity_ / 10) % 10) > 2)
0854 edm::LogVerbatim("HBHEMuon") << "actHE.size() " << actHE.size();
0855 for (unsigned int i = 0; i < actHE.size(); ++i) {
0856 if (ieta == actHE[i].ieta && depth == actHE[i].depth) {
0857 lx = actHE[i].thick;
0858 break;
0859 }
0860 }
0861 }
0862 if (((verbosity_ / 10) % 10) > 2)
0863 edm::LogVerbatim("HBHEMuon") << "active thick " << lx;
0864 return lx;
0865 }
0866
0867
0868 DEFINE_FWK_MODULE(HcalRaddamMuon);