Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-23 01:39:03

0001 // -*- C++ -*-
0002 //
0003 // Package:     HcalTestBeam
0004 // Class  :     HcalTB06Analysis
0005 //
0006 // Implementation:
0007 //     Main analysis class for Hcal Test Beam 2006 Analysis
0008 //
0009 // Original Author:
0010 //         Created:  19 November 2015
0011 //
0012 
0013 // user include files
0014 #include "SimG4CMS/HcalTestBeam/interface/HcalTB06Histo.h"
0015 #include "SimG4CMS/HcalTestBeam/interface/HcalTB06BeamSD.h"
0016 
0017 // to retreive hits
0018 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0019 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0020 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0021 #include "DataFormats/Common/interface/Handle.h"
0022 
0023 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0024 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0025 
0026 #include "FWCore/PluginManager/interface/ModuleDef.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/Run.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 
0037 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0038 
0039 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
0040 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0041 
0042 #include "Randomize.hh"
0043 #include "globals.hh"
0044 
0045 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0046 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0047 
0048 // system include files
0049 #include <memory>
0050 #include <string>
0051 #include <vector>
0052 
0053 //#define EDM_ML_DEBUG
0054 
0055 class HcalTB06Analysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0056 public:
0057   explicit HcalTB06Analysis(const edm::ParameterSet& p);
0058   ~HcalTB06Analysis() override = default;
0059 
0060   void beginJob() override;
0061   void endJob() override;
0062   void analyze(const edm::Event& e, const edm::EventSetup& c) override;
0063 
0064   HcalTB06Analysis(const HcalTB06Analysis&) = delete;
0065   const HcalTB06Analysis& operator=(const HcalTB06Analysis&) = delete;
0066 
0067 private:
0068   const bool m_ECAL;
0069   const double m_eta;
0070   const double m_phi;
0071   const double m_ener;
0072   const std::vector<int> m_PDG;
0073 
0074   edm::EDGetTokenT<edm::PCaloHitContainer> m_EcalToken;
0075   const edm::EDGetTokenT<edm::PCaloHitContainer> m_HcalToken;
0076   const edm::EDGetTokenT<edm::PCaloHitContainer> m_BeamToken;
0077 
0078   const edm::ParameterSet m_ptb;
0079   const double m_timeLimit;
0080   const double m_widthEcal;
0081   const double m_widthHcal;
0082   const double m_factEcal;
0083   const double m_factHcal;
0084   const double m_eMIP;
0085 
0086   int count;
0087   int m_idxetaEcal;
0088   int m_idxphiEcal;
0089   int m_idxetaHcal;
0090   int m_idxphiHcal;
0091 
0092   std::unique_ptr<HcalTB06Histo> m_histo;
0093 };
0094 
0095 HcalTB06Analysis::HcalTB06Analysis(const edm::ParameterSet& p)
0096     : m_ECAL(p.getParameter<bool>("ECAL")),
0097       m_eta(p.getParameter<double>("MinEta")),
0098       m_phi(p.getParameter<double>("MinPhi")),
0099       m_ener(p.getParameter<double>("MinE")),
0100       m_PDG(p.getParameter<std::vector<int> >("PartID")),
0101       m_HcalToken(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"))),
0102       m_BeamToken(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalTB06BeamHits"))),
0103       m_ptb(p.getParameter<edm::ParameterSet>("TestBeamAnalysis")),
0104       m_timeLimit(m_ptb.getParameter<double>("TimeLimit")),
0105       m_widthEcal(m_ptb.getParameter<double>("EcalWidth")),
0106       m_widthHcal(m_ptb.getParameter<double>("HcalWidth")),
0107       m_factEcal(m_ptb.getParameter<double>("EcalFactor")),
0108       m_factHcal(m_ptb.getParameter<double>("HcalFactor")),
0109       m_eMIP(m_ptb.getParameter<double>("MIP")),
0110       count(0) {
0111   usesResource("TFileService");
0112   if (m_ECAL)
0113     m_EcalToken = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
0114   double minEta = p.getParameter<double>("MinEta");
0115   double maxEta = p.getParameter<double>("MaxEta");
0116   double minPhi = p.getParameter<double>("MinPhi");
0117   double maxPhi = p.getParameter<double>("MaxPhi");
0118   double beamEta = (maxEta + minEta) * 0.5;
0119   double beamPhi = (maxPhi + minPhi) * 0.5;
0120   if (beamPhi < 0) {
0121     beamPhi += twopi;
0122   }
0123 
0124   m_idxetaEcal = 13;
0125   m_idxphiEcal = 13;
0126 
0127   m_idxetaHcal = static_cast<int>(beamEta / 0.087) + 1;
0128   m_idxphiHcal = static_cast<int>(beamPhi / 0.087) + 6;
0129   if (m_idxphiHcal > 72) {
0130     m_idxphiHcal -= 73;
0131   }
0132 
0133   edm::LogInfo("HcalTB06Analysis") << "Beam parameters: E(GeV)= " << m_ener << " pdgID= " << m_PDG[0]
0134                                    << "\n eta= " << m_eta << "  idx_etaEcal= " << m_idxetaEcal
0135                                    << "  idx_etaHcal= " << m_idxetaHcal << "  phi= " << m_phi
0136                                    << "  idx_phiEcal= " << m_idxphiEcal << "  idx_phiHcal= " << m_idxphiHcal
0137                                    << "\n        EcalFactor= " << m_factEcal << "  EcalWidth= " << m_widthEcal << " GeV"
0138                                    << "\n        HcalFactor= " << m_factHcal << "  HcalWidth= " << m_widthHcal << " GeV"
0139                                    << "  MIP=       " << m_eMIP << " GeV"
0140                                    << "\n        TimeLimit=  " << m_timeLimit << " ns"
0141                                    << "\n";
0142   m_histo = std::make_unique<HcalTB06Histo>(m_ptb);
0143 }
0144 
0145 void HcalTB06Analysis::beginJob() { edm::LogInfo("HcalTB06Analysis") << " =====> Begin of Run"; }
0146 
0147 void HcalTB06Analysis::endJob() {
0148   edm::LogInfo("HcalTB06Analysis") << " =====> End of Run; Total number of events: " << count;
0149 }
0150 
0151 void HcalTB06Analysis::analyze(const edm::Event& evt, const edm::EventSetup&) {
0152   ++count;
0153 
0154   //Beam Information
0155   m_histo->fillPrimary(m_ener, m_eta, m_phi);
0156 
0157   std::vector<double> eCalo(6, 0), eTrig(7, 0);
0158 
0159   const std::vector<PCaloHit>* EcalHits = nullptr;
0160   if (m_ECAL) {
0161     const edm::Handle<edm::PCaloHitContainer>& Ecal = evt.getHandle(m_EcalToken);
0162     EcalHits = Ecal.product();
0163   }
0164   const edm::Handle<edm::PCaloHitContainer>& Hcal = evt.getHandle(m_HcalToken);
0165   const std::vector<PCaloHit>* HcalHits = Hcal.product();
0166   const edm::Handle<edm::PCaloHitContainer>& Beam = evt.getHandle(m_BeamToken);
0167   const std::vector<PCaloHit>* BeamHits = Beam.product();
0168 
0169   // Total Energy
0170   double eecals = 0.;
0171   double ehcals = 0.;
0172 
0173   unsigned int ne = 0;
0174   unsigned int nh = 0;
0175   if (m_ECAL) {
0176     ne = EcalHits->size();
0177     for (unsigned int i = 0; i < ne; ++i) {
0178       EBDetId ecalid((*EcalHits)[i].id());
0179 #ifdef EDM_ML_DEBUG
0180       edm::LogVerbatim("HcalTBSim") << "EB " << i << " " << ecalid.ieta() << ":" << m_idxetaEcal << "   "
0181                                     << ecalid.iphi() << ":" << m_idxphiEcal << "   " << (*EcalHits)[i].time() << ":"
0182                                     << m_timeLimit << "   " << (*EcalHits)[i].energy();
0183 #endif
0184       // 7x7 crystal selection
0185       if (std::abs(m_idxetaEcal - ecalid.ieta()) <= 3 && std::abs(m_idxphiEcal - ecalid.iphi()) <= 3 &&
0186           (*EcalHits)[i].time() < m_timeLimit) {
0187         eCalo[0] += (*EcalHits)[i].energy();
0188       }
0189     }
0190     if (m_widthEcal > 0.0) {
0191       eCalo[1] = G4RandGauss::shoot(0.0, m_widthEcal);
0192     }
0193     eecals = m_factEcal * (eCalo[0] + eCalo[1]);
0194   }
0195   if (HcalHits) {
0196     nh = HcalHits->size();
0197     for (unsigned int i = 0; i < nh; ++i) {
0198       HcalDetId hcalid((*HcalHits)[i].id());
0199 #ifdef EDM_ML_DEBUG
0200       edm::LogVerbatim("HcalTBSim") << "HC " << i << " " << hcalid.subdet() << "  " << hcalid.ieta() << ":"
0201                                     << m_idxetaHcal << "   " << hcalid.iphi() << ":" << m_idxphiHcal << "   "
0202                                     << (*HcalHits)[i].time() << ":" << m_timeLimit << "   " << (*HcalHits)[i].energy();
0203 #endif
0204       // 3x3 towers selection
0205       if (std::abs(m_idxetaHcal - hcalid.ieta()) <= 1 && std::abs(m_idxphiHcal - hcalid.iphi()) <= 1 &&
0206           (*HcalHits)[i].time() < m_timeLimit) {
0207         if (hcalid.subdet() != HcalOuter) {
0208           eCalo[2] += (*HcalHits)[i].energy();
0209         } else {
0210           eCalo[4] += (*HcalHits)[i].energy();
0211         }
0212       }
0213     }
0214     if (m_widthHcal > 0.0) {
0215       eCalo[3] = G4RandGauss::shoot(0.0, m_widthHcal);
0216       eCalo[5] = G4RandGauss::shoot(0.0, m_widthHcal);
0217     }
0218     ehcals = m_factHcal * eCalo[2] + eCalo[3];
0219   }
0220   double etots = eecals + ehcals;
0221 
0222   edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Etot(MeV)= " << etots << "   E(Ecal)= " << eecals
0223                             << "   E(Hcal)= " << ehcals << "  Nhits(ECAL)= " << ne << "  Nhits(HCAL)= " << nh;
0224   m_histo->fillEdep(etots, eecals, ehcals);
0225 
0226   if (BeamHits) {
0227     for (unsigned int i = 0; i < BeamHits->size(); ++i) {
0228       unsigned int id = ((*BeamHits)[i].id());
0229       int det, lay, ix, iy;
0230       HcalTestBeamNumbering::unpackIndex(id, det, lay, ix, iy);
0231       if ((det == 1) && ((*BeamHits)[i].time() < m_timeLimit)) {
0232         if (lay > 0 && lay <= 4) {
0233           eTrig[lay - 1] += (*BeamHits)[i].energy();
0234         } else if (lay == 7 || lay == 8) {
0235           eTrig[lay - 2] += (*BeamHits)[i].energy();
0236         } else if (lay >= 11 && lay <= 14) {
0237           eTrig[4] += (*BeamHits)[i].energy();
0238         }
0239       }
0240     }
0241   }
0242 
0243   edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Trigger Info: " << eTrig[0] << ":" << eTrig[1] << ":" << eTrig[2]
0244                             << ":" << eTrig[3] << ":" << eTrig[4] << ":" << eTrig[5] << ":" << eTrig[6];
0245 
0246   m_histo->fillTree(eCalo, eTrig);
0247 }
0248 
0249 #include "FWCore/PluginManager/interface/ModuleDef.h"
0250 
0251 DEFINE_FWK_MODULE(HcalTB06Analysis);