File indexing completed on 2024-04-06 12:29:50
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009
0010 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0011 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0012 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0013 #include "SimDataFormats/Track/interface/SimTrack.h"
0014 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0015
0016
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023
0024 #include <iostream>
0025 #include <fstream>
0026 #include <vector>
0027 #include <map>
0028 #include <string>
0029 #include <cmath>
0030
0031 #include "TH1F.h"
0032 #include "TProfile.h"
0033
0034 class HFPMTHitAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0035 public:
0036 explicit HFPMTHitAnalyzer(const edm::ParameterSet &);
0037 ~HFPMTHitAnalyzer() override {}
0038 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0039
0040 private:
0041 void beginJob() override;
0042 void beginRun(edm::Run const &, edm::EventSetup const &) override {}
0043 void analyze(const edm::Event &, const edm::EventSetup &) override;
0044 void endJob() override;
0045 void endRun(edm::Run const &, edm::EventSetup const &) override {}
0046 void analyzeHits(std::vector<PCaloHit> &, const std::vector<SimTrack> &);
0047
0048
0049 const std::string g4Label_, hcalHits_;
0050 const edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
0051 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0052 const edm::EDGetTokenT<edm::SimTrackContainer> tok_track_;
0053
0054 int event_no;
0055
0056
0057 TH1F *h_HFDepHit, *hHF_MC_e, *h_HFEta[3], *h_HFPhi[3];
0058
0059 TH1F *hHF_e_1[3], *hHF_em_1[3], *hHF_had_1[3];
0060 TH1F *hHF_e_2[3], *hHF_em_2[3], *hHF_had_2[3];
0061 TH1F *hHF_e_12[3], *hHF_em_12[3], *hHF_had_12[3];
0062
0063 TH1F *hHF1_time[3], *hHF1_time_Ewt[3];
0064 TH1F *hHF2_time[3], *hHF2_time_Ewt[3];
0065 TH1F *hHF12_time[3], *hHF12_time_Ewt[3];
0066 };
0067
0068 HFPMTHitAnalyzer::HFPMTHitAnalyzer(const edm::ParameterSet &iConfig)
0069 : g4Label_(iConfig.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
0070 hcalHits_(iConfig.getUntrackedParameter<std::string>("HitCollection", "HcalHits")),
0071 tok_evt_(consumes<edm::HepMCProduct>(
0072 edm::InputTag(iConfig.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")))),
0073 tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_))),
0074 tok_track_(consumes<edm::SimTrackContainer>(edm::InputTag(g4Label_))) {
0075 usesResource(TFileService::kSharedResource);
0076 }
0077
0078 void HFPMTHitAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0079 edm::ParameterSetDescription desc;
0080 desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
0081 desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0082 desc.addUntracked<std::string>("HitCollection", "HcalHits");
0083 descriptions.add("HFPMTHitAnalyzer", desc);
0084 }
0085
0086 void HFPMTHitAnalyzer::beginJob() {
0087 event_no = 0;
0088 char name[20], title[120], sub[11];
0089
0090 edm::Service<TFileService> fs;
0091 if (!fs.isAvailable())
0092 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0093 << "please add it to config file";
0094
0095 TFileDirectory HFHitsDir = fs->mkdir("HFPMTHits");
0096 h_HFDepHit = HFHitsDir.make<TH1F>("Hit20", "Depths in HF", 50, 0., 50.);
0097 h_HFDepHit->GetXaxis()->SetTitle("Depths in HF");
0098 for (int i = 0; i < 3; i++) {
0099 if (i == 0)
0100 sprintf(sub, "(PMT)");
0101 else if (i == 1)
0102 sprintf(sub, "(Bundle)");
0103 else
0104 sprintf(sub, "(Jungle)");
0105 sprintf(name, "Eta%d", i);
0106 sprintf(title, "Eta Index of hits in %s", sub);
0107 h_HFEta[i] = HFHitsDir.make<TH1F>(name, title, 100, 0, 100.);
0108 h_HFEta[i]->GetXaxis()->SetTitle(title);
0109 sprintf(name, "Phi%d", i);
0110 sprintf(title, "Phi Index of hits in %s", sub);
0111 h_HFPhi[i] = HFHitsDir.make<TH1F>(name, title, 100, 0, 100.);
0112 h_HFPhi[i]->GetXaxis()->SetTitle(title);
0113 }
0114 TFileDirectory HFSourcePart = fs->mkdir("HFMCinfo");
0115 hHF_MC_e = HFSourcePart.make<TH1F>("MCEnergy", "Energy of Generated Particle", 1000, 0., 500.);
0116 hHF_MC_e->GetXaxis()->SetTitle("Energy of Generated Particle");
0117
0118
0119 TFileDirectory HFPCaloHitEnergyDir = fs->mkdir("HFPCaloHitEnergy");
0120 for (int i = 0; i < 3; i++) {
0121 if (i == 0)
0122 sprintf(sub, "(Absorber)");
0123 else if (i == 1)
0124 sprintf(sub, "(PMT)");
0125 else
0126 sprintf(sub, "(All)");
0127 sprintf(name, "Energy1%d", i);
0128 sprintf(title, "Energy in depth 1 %s", sub);
0129 hHF_e_1[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
0130 hHF_e_1[i]->GetXaxis()->SetTitle(title);
0131 sprintf(name, "Energy2%d", i);
0132 sprintf(title, "Energy in depth 2 %s", sub);
0133 hHF_e_2[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
0134 hHF_e_2[i]->GetXaxis()->SetTitle(title);
0135 sprintf(name, "Energy12%d", i);
0136 sprintf(title, "Energy in depths 1,2 %s", sub);
0137 hHF_e_12[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
0138 hHF_e_12[i]->GetXaxis()->SetTitle(title);
0139 sprintf(name, "Em1%d", i);
0140 sprintf(title, "EM energy in depth 1 %s", sub);
0141 hHF_em_1[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
0142 hHF_em_1[i]->GetXaxis()->SetTitle(title);
0143 sprintf(name, "Em2%d", i);
0144 sprintf(title, "EM energy in depth 2 %s", sub);
0145 hHF_em_2[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
0146 hHF_em_2[i]->GetXaxis()->SetTitle(title);
0147 sprintf(name, "Em12%d", i);
0148 sprintf(title, "EM energy in depths 1,2 %s", sub);
0149 hHF_em_12[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
0150 hHF_em_12[i]->GetXaxis()->SetTitle(title);
0151 sprintf(name, "Had1%d", i);
0152 sprintf(title, "Had energy in depth 1 %s", sub);
0153 hHF_had_1[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 0.1);
0154 hHF_had_1[i]->GetXaxis()->SetTitle(title);
0155 sprintf(name, "Had2%d", i);
0156 sprintf(title, "Had energy in depth 2 %s", sub);
0157 hHF_had_2[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 0.1);
0158 hHF_had_2[i]->GetXaxis()->SetTitle(title);
0159 sprintf(name, "Had12%d", i);
0160 sprintf(title, "Had energy in depths 1,2 %s", sub);
0161 hHF_had_12[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 0.1);
0162 hHF_had_12[i]->GetXaxis()->SetTitle(title);
0163 }
0164
0165
0166 TFileDirectory HFPCaloHitTimeDir = fs->mkdir("HFPCaloHitTime");
0167 for (int i = 0; i < 3; i++) {
0168 if (i == 0)
0169 sprintf(sub, "(Absorber)");
0170 else if (i == 1)
0171 sprintf(sub, "(PMT)");
0172 else
0173 sprintf(sub, "(All)");
0174 sprintf(name, "Time1Ewt%d", i);
0175 sprintf(title, "Time (energy weighted) in depth 1 %s", sub);
0176 hHF1_time_Ewt[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
0177 hHF1_time_Ewt[i]->GetXaxis()->SetTitle(title);
0178 sprintf(name, "Time2Ewt%d", i);
0179 sprintf(title, "Time (energy weighted) in depth 2 %s", sub);
0180 hHF2_time_Ewt[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
0181 hHF2_time_Ewt[i]->GetXaxis()->SetTitle(title);
0182 sprintf(name, "Time12Ewt%d", i);
0183 sprintf(title, "Time (energy weighted) in depths 1,2 %s", sub);
0184 hHF12_time_Ewt[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
0185 hHF12_time_Ewt[i]->GetXaxis()->SetTitle(title);
0186 sprintf(name, "Time1%d", i);
0187 sprintf(title, "Time in depth 1 %s", sub);
0188 hHF1_time[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
0189 hHF1_time[i]->GetXaxis()->SetTitle(title);
0190 sprintf(name, "Time2%d", i);
0191 sprintf(title, "Time in depth 2 %s", sub);
0192 hHF2_time[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
0193 hHF2_time[i]->GetXaxis()->SetTitle(title);
0194 sprintf(name, "Time12%d", i);
0195 sprintf(title, "Time in depths 1,2 %s", sub);
0196 hHF12_time[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
0197 hHF12_time[i]->GetXaxis()->SetTitle(title);
0198 }
0199 }
0200
0201 void HFPMTHitAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0202 ++event_no;
0203 if (event_no % 500 == 0)
0204 edm::LogVerbatim("HcalSim") << "Event # " << event_no << " processed.";
0205
0206 const edm::Handle<edm::PCaloHitContainer> &hitsHcal = iEvent.getHandle(tok_calo_);
0207 const edm::Handle<edm::SimTrackContainer> &Tracks = iEvent.getHandle(tok_track_);
0208 const edm::Handle<edm::HepMCProduct> EvtHandle = iEvent.getHandle(tok_evt_);
0209 const HepMC::GenEvent *myGenEvent = EvtHandle->GetEvent();
0210
0211 float orig_energy = 0;
0212 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0213 ++p) {
0214 orig_energy = (*p)->momentum().e();
0215 hHF_MC_e->Fill(orig_energy);
0216 }
0217
0218 std::vector<PCaloHit> caloHits;
0219 caloHits.insert(caloHits.end(), hitsHcal->begin(), hitsHcal->end());
0220 analyzeHits(caloHits, *Tracks);
0221 }
0222
0223 void HFPMTHitAnalyzer::analyzeHits(std::vector<PCaloHit> &hits, const std::vector<SimTrack> &tracks1) {
0224 int nHit = hits.size();
0225 float energy1[3], energy2[3], energy12[3];
0226 float em1[3], had1[3], em2[3], had2[3], em12[3], had12[3];
0227 for (int i = 0; i < 3; i++) {
0228 energy1[i] = energy2[i] = energy12[i] = 0;
0229 em1[i] = em2[i] = em12[i] = 0;
0230 had1[i] = had2[i] = had12[i] = 0;
0231 }
0232 edm::LogVerbatim("HFShower") << "HFPMTHitAnalyser::Entry " << event_no << " with " << nHit << " hits";
0233 for (int i = 0; i < nHit; i++) {
0234 double energy = hits[i].energy();
0235 double em = hits[i].energyEM();
0236 double had = hits[i].energyHad();
0237 double time = hits[i].time();
0238 uint32_t id_ = hits[i].id();
0239 uint16_t pmtHit = hits[i].depth();
0240 uint16_t depthX = pmtHit;
0241 edm::LogVerbatim("HFShower") << "HFPMTHitAnalyser::Hit " << i << " ID " << std::hex << id_ << std::dec << " PMT "
0242 << pmtHit << " E (e|h|t) " << em << " " << had << " " << energy << " Time " << time;
0243
0244 HcalDetId id = HcalDetId(id_);
0245 int det = id.det();
0246 int subdet = id.subdet();
0247 int depth = id.depth();
0248 if (pmtHit != 0)
0249 pmtHit = 1;
0250
0251 if (det == 4) {
0252 if (subdet == static_cast<int>(HcalForward)) {
0253 h_HFDepHit->Fill(double(depth + 10 * depthX));
0254 if (depthX > 0) {
0255 int ieta = id.ietaAbs();
0256 int iphi = id.iphi();
0257 if (depth != 1) {
0258 ieta += 50;
0259 iphi += 50;
0260 }
0261 if (depthX <= 3) {
0262 h_HFEta[depthX - 1]->Fill(double(ieta));
0263 h_HFPhi[depthX - 1]->Fill(double(iphi));
0264 }
0265 }
0266 if (depth == 1) {
0267 energy1[pmtHit] += energy;
0268 energy12[pmtHit] += energy;
0269 em1[pmtHit] += em;
0270 em12[pmtHit] += em;
0271 had1[pmtHit] += had;
0272 had12[pmtHit] += had;
0273 energy1[2] += energy;
0274 energy12[2] += energy;
0275 em1[2] += em;
0276 em12[2] += em;
0277 had1[2] += had;
0278 had12[2] += had;
0279 hHF1_time[pmtHit]->Fill(time);
0280 hHF1_time[2]->Fill(time);
0281 hHF1_time_Ewt[pmtHit]->Fill(time, energy);
0282 hHF1_time_Ewt[2]->Fill(time, energy);
0283 }
0284 if (depth == 2) {
0285 energy2[pmtHit] += energy;
0286 energy12[pmtHit] += energy;
0287 em2[pmtHit] += em;
0288 em12[pmtHit] += em;
0289 had2[pmtHit] += had;
0290 had12[pmtHit] += had;
0291 energy2[2] += energy;
0292 energy12[2] += energy;
0293 em2[2] += em;
0294 em12[2] += em;
0295 had2[2] += had;
0296 had12[2] += had;
0297 hHF2_time[pmtHit]->Fill(time);
0298 hHF2_time[2]->Fill(time);
0299 hHF2_time_Ewt[pmtHit]->Fill(time, energy);
0300 hHF2_time_Ewt[2]->Fill(time, energy);
0301 }
0302 }
0303 }
0304 }
0305 for (int i = 0; i < 3; i++) {
0306 hHF_e_1[i]->Fill(energy1[i]);
0307 hHF_e_2[i]->Fill(energy2[i]);
0308 hHF_e_12[i]->Fill(energy12[i]);
0309 hHF_em_1[i]->Fill(em1[i]);
0310 hHF_em_2[i]->Fill(em2[i]);
0311 hHF_em_12[i]->Fill(em12[i]);
0312 hHF_had_1[i]->Fill(had1[i]);
0313 hHF_had_2[i]->Fill(had2[i]);
0314 hHF_had_12[i]->Fill(had12[i]);
0315 edm::LogVerbatim("HFShower") << "HFPMTHitAnalyser:: Type " << i << " Energy 1|2| " << energy1[i] << " "
0316 << energy2[i] << " " << energy12[i] << " EM Energy 1|2| " << em1[i] << " " << em2[i]
0317 << " " << em12[i] << " Had Energy 1|2| " << had1[i] << " " << had2[i] << " "
0318 << had12[i];
0319 }
0320 }
0321
0322 void HFPMTHitAnalyzer::endJob() {}
0323
0324 DEFINE_FWK_MODULE(HFPMTHitAnalyzer);
