Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:58

0001 // system include files
0002 #include <memory>
0003 #include <string>
0004 #include <vector>
0005 
0006 // user include files
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 
0021 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 
0024 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0025 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0026 #include "SimDataFormats/Track/interface/SimTrack.h"
0027 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0028 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0029 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0030 
0031 #include <TH1F.h>
0032 #include <TH1I.h>
0033 
0034 //#define EDM_ML_DEBUG
0035 
0036 class XtalDedxAnalysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0037 public:
0038   explicit XtalDedxAnalysis(const edm::ParameterSet &);
0039   ~XtalDedxAnalysis() override {}
0040 
0041   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0042 
0043 protected:
0044   void beginJob() override {}
0045   void analyze(const edm::Event &, const edm::EventSetup &) override;
0046   void endJob() override {}
0047 
0048   void analyzeHits(std::vector<PCaloHit> &,
0049                    const edm::Handle<edm::SimTrackContainer> &,
0050                    const edm::Handle<edm::SimVertexContainer> &);
0051 
0052 private:
0053   const edm::InputTag caloHitSource_;
0054   const std::string simTkLabel_;
0055   const double energyMax_;
0056   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0057   const edm::EDGetTokenT<edm::SimTrackContainer> tok_tk_;
0058   const edm::EDGetTokenT<edm::SimVertexContainer> tok_vtx_;
0059 
0060   TH1F *meNHit_[4], *meE1T0_[4], *meE9T0_[4], *meE1T1_[4], *meE9T1_[4];
0061   TH1I *mType_;
0062 };
0063 
0064 XtalDedxAnalysis::XtalDedxAnalysis(const edm::ParameterSet &ps)
0065     : caloHitSource_(ps.getParameter<edm::InputTag>("caloHitSource")),
0066       simTkLabel_(ps.getParameter<std::string>("moduleLabelTk")),
0067       energyMax_(ps.getParameter<double>("energyMax")),
0068       tok_calo_(consumes<edm::PCaloHitContainer>(caloHitSource_)),
0069       tok_tk_(consumes<edm::SimTrackContainer>(edm::InputTag(simTkLabel_))),
0070       tok_vtx_(consumes<edm::SimVertexContainer>(edm::InputTag(simTkLabel_))) {
0071   usesResource(TFileService::kSharedResource);
0072   edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis::Source " << caloHitSource_ << " Track Label "
0073                                         << simTkLabel_ << " Energy Max " << energyMax_;
0074 
0075   // Book histograms
0076   edm::Service<TFileService> tfile;
0077 
0078   if (!tfile.isAvailable())
0079     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0080                                       << "please add it to config file";
0081   // Histograms for Hits
0082   std::string types[4] = {"total", "by dE/dx", "by delta-ray", "by bremms"};
0083   char name[20], title[80];
0084   for (int i = 0; i < 4; i++) {
0085     sprintf(name, "Hits%d", i);
0086     sprintf(title, "Number of hits (%s)", types[i].c_str());
0087     meNHit_[i] = tfile->make<TH1F>(name, title, 5000, 0., 5000.);
0088     meNHit_[i]->GetXaxis()->SetTitle(title);
0089     meNHit_[i]->GetYaxis()->SetTitle("Events");
0090     sprintf(name, "E1T0%d", i);
0091     sprintf(title, "E1 (Loss %s) in GeV", types[i].c_str());
0092     meE1T0_[i] = tfile->make<TH1F>(name, title, 5000, 0, energyMax_);
0093     meE1T0_[i]->GetXaxis()->SetTitle(title);
0094     meE1T0_[i]->GetYaxis()->SetTitle("Events");
0095     sprintf(name, "E9T0%d", i);
0096     sprintf(title, "E9 (Loss %s) in GeV", types[i].c_str());
0097     meE9T0_[i] = tfile->make<TH1F>(name, title, 5000, 0, energyMax_);
0098     meE9T0_[i]->GetXaxis()->SetTitle(title);
0099     meE9T0_[i]->GetYaxis()->SetTitle("Events");
0100     sprintf(name, "E1T1%d", i);
0101     sprintf(title, "E1 (Loss %s with t < 400 ns) in GeV", types[i].c_str());
0102     meE1T1_[i] = tfile->make<TH1F>(name, title, 5000, 0, energyMax_);
0103     meE1T1_[i]->GetXaxis()->SetTitle(title);
0104     meE1T1_[i]->GetYaxis()->SetTitle("Events");
0105     sprintf(name, "E9T1%d", i);
0106     sprintf(title, "E9 (Loss %s with t < 400 ns) in GeV", types[i].c_str());
0107     meE9T1_[i] = tfile->make<TH1F>(name, title, 5000, 0, energyMax_);
0108     meE9T1_[i]->GetXaxis()->SetTitle(title);
0109     meE9T1_[i]->GetYaxis()->SetTitle("Events");
0110   }
0111   sprintf(name, "PDGType");
0112   sprintf(title, "PDG ID of first level secondary");
0113   mType_ = tfile->make<TH1I>(name, title, 5000, -2500, 2500);
0114   mType_->GetXaxis()->SetTitle(title);
0115   mType_->GetYaxis()->SetTitle("Tracks");
0116 }
0117 
0118 void XtalDedxAnalysis::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0119   edm::ParameterSetDescription desc;
0120   desc.add<edm::InputTag>("caloHitSource", edm::InputTag("g4SimHits", "EcalHitsEB"));
0121   desc.add<std::string>("moduleLabelTk", "g4SimHits");
0122   desc.add<double>("energyMax", 2.0);
0123   descriptions.add("xtalDedxAnalysis", desc);
0124 }
0125 
0126 void XtalDedxAnalysis::analyze(const edm::Event &e, const edm::EventSetup &) {
0127 #ifdef EDM_ML_DEBUG
0128   edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis::Run = " << e.id().run() << " Event = " << e.id().event();
0129 #endif
0130   std::vector<PCaloHit> caloHits;
0131   const edm::Handle<edm::PCaloHitContainer> &pCaloHits = e.getHandle(tok_calo_);
0132 
0133   const edm::Handle<edm::SimTrackContainer> &simTk = e.getHandle(tok_tk_);
0134   const edm::Handle<edm::SimVertexContainer> &simVtx = e.getHandle(tok_vtx_);
0135 
0136   if (pCaloHits.isValid()) {
0137     caloHits.insert(caloHits.end(), pCaloHits->begin(), pCaloHits->end());
0138 #ifdef EDM_ML_DEBUG
0139     edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis: Hit buffer " << caloHits.size();
0140 #endif
0141     analyzeHits(caloHits, simTk, simVtx);
0142   }
0143 }
0144 
0145 void XtalDedxAnalysis::analyzeHits(std::vector<PCaloHit> &hits,
0146                                    const edm::Handle<edm::SimTrackContainer> &SimTk,
0147                                    const edm::Handle<edm::SimVertexContainer> &SimVtx) {
0148   edm::SimTrackContainer::const_iterator simTrkItr;
0149   int nHit = hits.size();
0150   double e10[4], e90[4], e11[4], e91[4], hit[4];
0151   for (int i = 0; i < 4; i++)
0152     e10[i] = e90[i] = e11[i] = e91[i] = hit[i] = 0;
0153   for (int i = 0; i < nHit; i++) {
0154     double energy = hits[i].energy();
0155     double time = hits[i].time();
0156     unsigned int id_ = hits[i].id();
0157     int trackID = hits[i].geantTrackId();
0158     int type = 1;
0159     for (simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0160       if (trackID == (int)(simTrkItr->trackId())) {
0161         int thePID = simTrkItr->type();
0162         if (thePID == 11)
0163           type = 2;
0164         else if (thePID != -13 && thePID != 13)
0165           type = 3;
0166         break;
0167       }
0168     }
0169     hit[0]++;
0170     hit[type]++;
0171     e90[0] += energy;
0172     e90[type] += energy;
0173     if (time < 400) {
0174       e91[0] += energy;
0175       e91[type] += energy;
0176     }
0177     if (id_ == 22) {
0178       e10[0] += energy;
0179       e10[type] += energy;
0180       if (time < 400) {
0181         e11[0] += energy;
0182         e11[type] += energy;
0183       }
0184     }
0185 #ifdef EDM_ML_DEBUG
0186     edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis:Hit[" << i << "] ID " << id_ << " E " << energy
0187                                           << " time " << time << " track " << trackID << " type " << type;
0188 #endif
0189   }
0190   for (int i = 0; i < 4; i++) {
0191 #ifdef EDM_ML_DEBUG
0192     edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis:Type(" << i << ") Hit " << hit[i] << " E10 " << e10[i]
0193                                           << " E11 " << e11[i] << " E90 " << e90[i] << " E91 " << e91[i];
0194 #endif
0195     meNHit_[i]->Fill(hit[i]);
0196     meE1T0_[i]->Fill(e10[i]);
0197     meE9T0_[i]->Fill(e90[i]);
0198     meE1T1_[i]->Fill(e11[i]);
0199     meE9T1_[i]->Fill(e91[i]);
0200   }
0201 
0202   // Type of the secondary (coming directly from a generator level track)
0203   int nvtx = 0, k1 = 0;
0204   edm::SimVertexContainer::const_iterator simVtxItr;
0205   for (simVtxItr = SimVtx->begin(); simVtxItr != SimVtx->end(); simVtxItr++)
0206     nvtx++;
0207 #ifdef EDM_ML_DEBUG
0208   int ntrk = 0;
0209   for (simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++)
0210     ntrk++;
0211   edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis: " << ntrk << " tracks and " << nvtx << " vertices";
0212 #endif
0213   for (simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++, ++k1) {
0214 #ifdef EDM_ML_DEBUG
0215     edm::LogVerbatim("CherenkovAnalysis") << "Track " << k1 << " PDGId " << simTrkItr->type() << " Vertex ID "
0216                                           << simTrkItr->vertIndex() << " Generator " << simTrkItr->noGenpart();
0217 #endif
0218     if (simTrkItr->noGenpart()) {              // This is a secondary
0219       int vertIndex = simTrkItr->vertIndex();  // Vertex index of origin
0220       if (vertIndex >= 0 && vertIndex < nvtx) {
0221         simVtxItr = SimVtx->begin();
0222         for (int iv = 0; iv < vertIndex; iv++)
0223           simVtxItr++;
0224         int parent = simVtxItr->parentIndex(), k2 = 0;
0225         for (edm::SimTrackContainer::const_iterator trkItr = SimTk->begin(); trkItr != SimTk->end(); trkItr++, ++k2) {
0226 #ifdef EDM_ML_DEBUG
0227           edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis::Track " << k2 << " ID " << trkItr->trackId()
0228                                                 << " (" << parent << ")  Generator " << trkItr->noGenpart();
0229 #endif
0230           if ((int)trkItr->trackId() == parent) {  // Parent track
0231             if (!trkItr->noGenpart()) {            // Generator level
0232 #ifdef EDM_ML_DEBUG
0233               edm::LogVerbatim("CherenkovAnalysis") << "XtalDedxAnalysis::Track found with ID " << simTrkItr->type();
0234 #endif
0235               mType_->Fill(simTrkItr->type());
0236             }
0237             break;
0238           }
0239         }
0240       }
0241     }
0242   }
0243 }
0244 
0245 // define this as a plug-in
0246 DEFINE_FWK_MODULE(XtalDedxAnalysis);