Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-26 00:02:02

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   // get Geometry, B-field, Topology
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   // get PCaloHits for ecal barrel
0171   edm::Handle<edm::PCaloHitContainer> caloHitEB;
0172   iEvent.getByToken(tok_EBhit_, caloHitEB);
0173 
0174   // get PCaloHits for ecal endcap
0175   edm::Handle<edm::PCaloHitContainer> caloHitEE;
0176   iEvent.getByToken(tok_EEhit_, caloHitEE);
0177 
0178   // get sim tracks
0179   edm::Handle<edm::SimTrackContainer> SimTk;
0180   iEvent.getByToken(tok_simTk_, SimTk);
0181 
0182   // get sim vertices
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 //define this as a plug-in
0295 DEFINE_FWK_MODULE(ElectronStudy);