Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-18 08:23:52

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 //SimHits
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 //other stuff
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   //user parameters
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   //root file name and its objects
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   //energy Histograms
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   //Timing Histograms
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);