File indexing completed on 2024-04-06 12:29:58
0001
0002 #include <memory>
0003 #include <string>
0004 #include <vector>
0005
0006
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
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
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
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
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()) {
0219 int vertIndex = simTrkItr->vertIndex();
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) {
0231 if (!trkItr->noGenpart()) {
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
0246 DEFINE_FWK_MODULE(XtalDedxAnalysis);