File indexing completed on 2024-04-06 11:59:13
0001 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0002 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0003
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015
0016 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0017 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0018 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0019 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0020 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0021 #include "MagneticField/Engine/interface/MagneticField.h"
0022 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0023
0024 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0025 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0026 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0027 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0028
0029 #include <TH1F.h>
0030
0031 #include <algorithm>
0032 #include <fstream>
0033 #include <iostream>
0034 #include <memory>
0035 #include <string>
0036 #include <vector>
0037
0038 class ElectronStudy : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0039 public:
0040 explicit ElectronStudy(const edm::ParameterSet& ps);
0041 ~ElectronStudy() override {}
0042
0043 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0044
0045 private:
0046 void analyze(edm::Event const&, edm::EventSetup const&) override;
0047 void beginJob() override {}
0048
0049 static const int NEtaBins_ = 3;
0050 static const int NPBins_ = 8;
0051 double pBins_[NPBins_ + 1], etaBins_[NEtaBins_ + 1];
0052
0053 edm::EDGetTokenT<edm::PCaloHitContainer> tok_EBhit_;
0054 edm::EDGetTokenT<edm::PCaloHitContainer> tok_EEhit_;
0055 edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0056 edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0057
0058 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0059 edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0060 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0061
0062 int hotZone_, verbose_;
0063 bool histos_;
0064 TH1F* histoR1[NPBins_ + 1][NEtaBins_ + 1];
0065 TH1F* histoR2[NPBins_ + 1][NEtaBins_ + 1];
0066 TH1F* histoR3[NPBins_ + 1][NEtaBins_ + 1];
0067 TH1F* histoE1x1[NPBins_ + 1][NEtaBins_ + 1];
0068 TH1F* histoE3x3[NPBins_ + 1][NEtaBins_ + 1];
0069 TH1F* histoE5x5[NPBins_ + 1][NEtaBins_ + 1];
0070 };
0071
0072 ElectronStudy::ElectronStudy(const edm::ParameterSet& ps) {
0073 usesResource("TFileService");
0074
0075 std::string g4Label = ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits");
0076 std::string hitLabEB = ps.getUntrackedParameter<std::string>("EBCollection", "EcalHitsEB");
0077 std::string hitLabEE = ps.getUntrackedParameter<std::string>("EECollection", "EcalHitsEE");
0078
0079 tok_EBhit_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hitLabEB));
0080 tok_EEhit_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hitLabEE));
0081 tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag(g4Label));
0082 tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag(g4Label));
0083
0084 hotZone_ = ps.getUntrackedParameter<int>("HotZone", 0);
0085 verbose_ = ps.getUntrackedParameter<int>("Verbosity", 0);
0086 edm::LogInfo("ElectronStudy") << "Module Label: " << g4Label << " Hits: " << hitLabEB << ", " << hitLabEE;
0087
0088 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0089 tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0090 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0091
0092 double tempP[NPBins_ + 1] = {0.0, 10.0, 20.0, 40.0, 60.0, 100.0, 500.0, 1000.0, 10000.0};
0093 double tempEta[NEtaBins_ + 1] = {0.0, 1.2, 1.6, 3.0};
0094
0095 for (int i = 0; i < NPBins_ + 1; i++)
0096 pBins_[i] = tempP[i];
0097 for (int i = 0; i < NEtaBins_ + 1; i++)
0098 etaBins_[i] = tempEta[i];
0099
0100 edm::Service<TFileService> tfile;
0101 if (!tfile.isAvailable()) {
0102 edm::LogInfo("ElectronStudy") << "TFileService unavailable: no histograms";
0103 histos_ = false;
0104 } else {
0105 char name[20], title[200], cpbin[30], cebin[30];
0106 histos_ = true;
0107 for (unsigned int i = 0; i < NPBins_ + 1; ++i) {
0108 if (i == 0)
0109 sprintf(cpbin, " All p");
0110 else
0111 sprintf(cpbin, " p (%6.0f:%6.0f)", pBins_[i - 1], pBins_[i]);
0112 for (unsigned int j = 0; j < NEtaBins_ + 1; ++j) {
0113 if (j == 0)
0114 sprintf(cebin, " All #eta");
0115 else
0116 sprintf(cebin, " #eta (%4.1f:%4.1f)", etaBins_[j - 1], etaBins_[j]);
0117 sprintf(name, "R1%d%d", i, j);
0118 sprintf(title, "E1/E9 for %s%s", cpbin, cebin);
0119 histoR1[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
0120 histoR1[i][j]->GetXaxis()->SetTitle(title);
0121 histoR1[i][j]->GetYaxis()->SetTitle("Tracks");
0122 sprintf(name, "R2%d%d", i, j);
0123 sprintf(title, "E1/E25 for %s%s", cpbin, cebin);
0124 histoR2[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
0125 histoR2[i][j]->GetXaxis()->SetTitle(title);
0126 histoR2[i][j]->GetYaxis()->SetTitle("Tracks");
0127 sprintf(name, "R3%d%d", i, j);
0128 sprintf(title, "E9/E25 for %s%s", cpbin, cebin);
0129 histoR3[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
0130 histoR3[i][j]->GetXaxis()->SetTitle(title);
0131 histoR3[i][j]->GetYaxis()->SetTitle("Tracks");
0132 sprintf(name, "E1x1%d%d", i, j);
0133 sprintf(title, "E1/P for %s%s", cpbin, cebin);
0134 histoE1x1[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
0135 histoE1x1[i][j]->GetXaxis()->SetTitle(title);
0136 histoE1x1[i][j]->GetYaxis()->SetTitle("Tracks");
0137 sprintf(name, "E3x3%d%d", i, j);
0138 sprintf(title, "E9/P for %s%s", cpbin, cebin);
0139 histoE3x3[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
0140 histoE3x3[i][j]->GetXaxis()->SetTitle(title);
0141 histoE3x3[i][j]->GetYaxis()->SetTitle("Tracks");
0142 sprintf(name, "E5x5%d%d", i, j);
0143 sprintf(title, "E25/P for %s%s", cpbin, cebin);
0144 histoE5x5[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
0145 histoE5x5[i][j]->GetXaxis()->SetTitle(title);
0146 histoE5x5[i][j]->GetYaxis()->SetTitle("Tracks");
0147 }
0148 }
0149 }
0150 }
0151
0152 void ElectronStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0153 edm::ParameterSetDescription desc;
0154 desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0155 desc.addUntracked<std::string>("EBCollection", "EcalHitsEB");
0156 desc.addUntracked<std::string>("EECollection", "EcalHitsEE");
0157 desc.addUntracked<int>("Verbosity", 0);
0158 descriptions.add("electronStudy", desc);
0159 }
0160
0161 void ElectronStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0162 if (verbose_ > 1)
0163 edm::LogVerbatim("IsoTrack") << "Run = " << iEvent.id().run() << " Event = " << iEvent.id().event();
0164
0165
0166 const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0167 const MagneticField* bField = &iSetup.getData(tok_magField_);
0168 const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0169
0170
0171 edm::Handle<edm::PCaloHitContainer> caloHitEB;
0172 iEvent.getByToken(tok_EBhit_, caloHitEB);
0173
0174
0175 edm::Handle<edm::PCaloHitContainer> caloHitEE;
0176 iEvent.getByToken(tok_EEhit_, caloHitEE);
0177
0178
0179 edm::Handle<edm::SimTrackContainer> SimTk;
0180 iEvent.getByToken(tok_simTk_, SimTk);
0181
0182
0183 edm::Handle<edm::SimVertexContainer> SimVtx;
0184 iEvent.getByToken(tok_simVtx_, SimVtx);
0185
0186 if (verbose_ > 0)
0187 edm::LogVerbatim("IsoTrack") << "ElectronStudy: hits valid[EB]: " << caloHitEB.isValid()
0188 << " valid[EE]: " << caloHitEE.isValid();
0189
0190 if (caloHitEB.isValid() && caloHitEE.isValid()) {
0191 unsigned int indx;
0192 if (verbose_ > 2) {
0193 edm::PCaloHitContainer::const_iterator ihit;
0194 for (ihit = caloHitEB->begin(), indx = 0; ihit != caloHitEB->end(); ihit++, indx++) {
0195 EBDetId id = ihit->id();
0196 edm::LogVerbatim("IsoTrack") << "Hit[" << indx << "] " << id << " E " << ihit->energy() << " T "
0197 << ihit->time();
0198 }
0199 for (ihit = caloHitEE->begin(), indx = 0; ihit != caloHitEE->end(); ihit++, indx++) {
0200 EEDetId id = ihit->id();
0201 edm::LogVerbatim("IsoTrack") << "Hit[" << indx << "] " << id << " E " << ihit->energy() << " T "
0202 << ihit->time();
0203 }
0204 }
0205 edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin();
0206 for (indx = 0; simTrkItr != SimTk->end(); simTrkItr++, indx++) {
0207 if (verbose_ > 0)
0208 edm::LogVerbatim("IsoTrack") << "ElectronStudy: Track[" << indx << "] ID " << simTrkItr->trackId() << " type "
0209 << simTrkItr->type() << " charge " << simTrkItr->charge() << " p "
0210 << simTrkItr->momentum() << " Generator Index " << simTrkItr->genpartIndex()
0211 << " vertex " << simTrkItr->vertIndex();
0212 if (std::abs(simTrkItr->type()) == 11 && simTrkItr->vertIndex() != -1) {
0213 int thisTrk = simTrkItr->trackId();
0214 spr::propagatedTrackDirection trkD = spr::propagateCALO(thisTrk, SimTk, SimVtx, geo, bField, (verbose_ > 1));
0215 if (trkD.okECAL) {
0216 const DetId isoCell = trkD.detIdECAL;
0217 DetId hotCell = isoCell;
0218 if (hotZone_ > 0)
0219 hotCell = spr::hotCrystal(
0220 isoCell, caloHitEB, caloHitEE, geo, caloTopology, hotZone_, hotZone_, -500.0, 500.0, (verbose_ > 1));
0221 double e1x1 = spr::eECALmatrix(
0222 hotCell, caloHitEB, caloHitEE, geo, caloTopology, 0, 0, -100.0, -100.0, -500.0, 500.0, (verbose_ > 2));
0223 double e3x3 = spr::eECALmatrix(
0224 hotCell, caloHitEB, caloHitEE, geo, caloTopology, 1, 1, -100.0, -100.0, -500.0, 500.0, (verbose_ > 2));
0225 double e5x5 = spr::eECALmatrix(
0226 hotCell, caloHitEB, caloHitEE, geo, caloTopology, 2, 2, -100.0, -100.0, -500.0, 500.0, (verbose_ > 2));
0227 double p = simTrkItr->momentum().P();
0228 double eta = std::abs(simTrkItr->momentum().eta());
0229 int etaBin = -1, momBin = -1;
0230 for (int ieta = 0; ieta < NEtaBins_; ieta++) {
0231 if (eta > etaBins_[ieta] && eta < etaBins_[ieta + 1])
0232 etaBin = ieta + 1;
0233 }
0234 for (int ipt = 0; ipt < NPBins_; ipt++) {
0235 if (p > pBins_[ipt] && p < pBins_[ipt + 1])
0236 momBin = ipt + 1;
0237 }
0238 double r1 = -1, r2 = -1, r3 = -1;
0239 if (e3x3 > 0)
0240 r1 = e1x1 / e3x3;
0241 if (e5x5 > 0) {
0242 r2 = e1x1 / e5x5;
0243 r3 = e3x3 / e5x5;
0244 }
0245 if (verbose_ > 0) {
0246 edm::LogVerbatim("IsoTrack") << "ElectronStudy: p " << p << " [" << momBin << "] eta " << eta << " ["
0247 << etaBin << "] Cell 0x" << std::hex << isoCell() << std::dec;
0248 if (isoCell.subdetId() == EcalBarrel) {
0249 edm::LogVerbatim("IsoTrack") << EBDetId(isoCell);
0250 } else if (isoCell.subdetId() == EcalEndcap) {
0251 edm::LogVerbatim("IsoTrack") << EEDetId(isoCell);
0252 }
0253 edm::LogVerbatim("IsoTrack") << " e1x1 " << e1x1 << "|" << r1 << "|" << r2 << " e3x3 " << e3x3 << "|" << r3
0254 << " e5x5 " << e5x5;
0255 }
0256 if (histos_) {
0257 histoR1[0][0]->Fill(r1);
0258 histoR2[0][0]->Fill(r2);
0259 histoR3[0][0]->Fill(r3);
0260 histoE1x1[0][0]->Fill(e1x1 / p);
0261 histoE3x3[0][0]->Fill(e3x3 / p);
0262 histoE5x5[0][0]->Fill(e5x5 / p);
0263 if (momBin > 0) {
0264 histoR1[momBin][0]->Fill(r1);
0265 histoR2[momBin][0]->Fill(r2);
0266 histoR3[momBin][0]->Fill(r3);
0267 histoE1x1[momBin][0]->Fill(e1x1 / p);
0268 histoE3x3[momBin][0]->Fill(e3x3 / p);
0269 histoE5x5[momBin][0]->Fill(e5x5 / p);
0270 }
0271 if (etaBin > 0) {
0272 histoR1[0][etaBin]->Fill(r1);
0273 histoR2[0][etaBin]->Fill(r2);
0274 histoR3[0][etaBin]->Fill(r3);
0275 histoE1x1[0][etaBin]->Fill(e1x1 / p);
0276 histoE3x3[0][etaBin]->Fill(e3x3 / p);
0277 histoE5x5[0][etaBin]->Fill(e5x5 / p);
0278 if (momBin > 0) {
0279 histoR1[momBin][etaBin]->Fill(r1);
0280 histoR2[momBin][etaBin]->Fill(r2);
0281 histoR3[momBin][etaBin]->Fill(r3);
0282 histoE1x1[momBin][etaBin]->Fill(e1x1 / p);
0283 histoE3x3[momBin][etaBin]->Fill(e3x3 / p);
0284 histoE5x5[momBin][etaBin]->Fill(e5x5 / p);
0285 }
0286 }
0287 }
0288 }
0289 }
0290 }
0291 }
0292 }
0293
0294
0295 DEFINE_FWK_MODULE(ElectronStudy);
