HcalSimAlgos/test/CaloSamplesAnalyzer.cc

0001 // system include files
0002 #include <memory>
0003
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Utilities/interface/StreamID.h"
0011
0012 //ROOT headers
0013 #include <TTree.h>
0014
0015 //calo headers
0016 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
0017 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
0018 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0019 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0020 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0021 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0022 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0023 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0024 #include "DataFormats/DetId/interface/DetId.h"
0025 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0026 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0027 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0028 #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
0029 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0030 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0031
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/PluginManager/interface/ModuleDef.h"
0034 #include "FWCore/ServiceRegistry/interface/Service.h"
0035 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0036
0037 //STL headers
0038 #include <vector>
0039 #include <unordered_map>
0040 #include <string>
0041
0042 //
0043 // class declaration
0044 //
0045
0046 class CaloSamplesAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0047 public:
0048   explicit CaloSamplesAnalyzer(const edm::ParameterSet&);
0049   ~CaloSamplesAnalyzer() override;
0050
0051   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0052
0053   struct CaloNtuple {
0054     unsigned run = 0;
0055     unsigned lumiblock = 0;
0056     unsigned long long event = 0;
0057     unsigned id = 0;
0058     int subdet = 0;
0059     int ieta = 0;
0060     int iphi = 0;
0061     int depth = 0;
0062     double samplingFactor = 0.;
0063     double fCtoGeV = 0.;
0064     double photoelectronsToAnalog = 0.;
0065     double simHitToPhotoelectrons = 0.;
0066     double tof = 0.;
0067     std::vector<double> energy;
0068     std::vector<double> time;
0069     std::vector<double> signalTot;
0070     std::vector<double> signalTotPrecise;
0071   };
0072
0073 private:
0074   void beginJob() override;
0075   void doBeginRun_(const edm::Run&, const edm::EventSetup&) override;
0076   void analyze(const edm::Event&, const edm::EventSetup&) override;
0077   void doEndRun_(const edm::Run&, const edm::EventSetup&) override {}
0078   void endJob() override {}
0079
0080   // ----------member data ---------------------------
0081   edm::Service<TFileService> fs;
0082   TTree* tree;
0083   const CaloGeometry* theGeometry;
0084   const HcalDDDRecConstants* theRecNumber;
0085   CaloHitResponse* theResponse;
0086   HcalSimParameterMap* theParameterMap;
0087   //for tree branches
0088   CaloNtuple entry;
0089   std::unordered_map<unsigned, CaloNtuple> treemap;
0090   //misc. params
0091   bool TestNumbering;
0092   //tokens
0093   edm::EDGetTokenT<std::vector<PCaloHit>> tok_sim;
0094   edm::EDGetTokenT<std::vector<CaloSamples>> tok_calo;
0095   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokGeom_;
0096   edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tokRecCns_;
0097   edm::ESGetToken<HcalDbService, HcalDbRecord> tokDB_;
0098 };
0099
0100 //
0101 // constructors and destructor
0102 //
0103 CaloSamplesAnalyzer::CaloSamplesAnalyzer(const edm::ParameterSet& iConfig)
0104     : tree(nullptr),
0105       theGeometry(nullptr),
0106       theRecNumber(nullptr),
0107       theResponse(new CaloHitResponse(nullptr, (CaloShapes*)nullptr)),
0108       theParameterMap(new HcalSimParameterMap(iConfig)),
0109       TestNumbering(iConfig.getParameter<bool>("TestNumbering")),
0110       tok_sim(consumes<std::vector<PCaloHit>>(
0111           edm::InputTag(iConfig.getParameter<std::string>("hitsProducer"), "HcalHits"))),
0112       tok_calo(consumes<std::vector<CaloSamples>>(iConfig.getParameter<edm::InputTag>("CaloSamplesTag"))),
0113       tokGeom_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0114       tokRecCns_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0115       tokDB_(esConsumes<HcalDbService, HcalDbRecord>()) {
0116   usesResource("TFileService");
0117 }
0118
0119 CaloSamplesAnalyzer::~CaloSamplesAnalyzer() {
0120   delete theResponse;
0121   delete theParameterMap;
0122 }
0123
0124 //
0125 // member functions
0126 //
0127
0128 void CaloSamplesAnalyzer::beginJob() {
0129   tree = fs->make<TTree>("tree", "tree");
0130
0131   tree->Branch("run", &entry.run, "run/i");
0132   tree->Branch("lumiblock", &entry.lumiblock, "lumiblock/i");
0133   tree->Branch("event", &entry.event, "event/l");
0134   tree->Branch("id", &entry.id, "id/i");
0135   tree->Branch("subdet", &entry.subdet, "subdet/I");
0136   tree->Branch("ieta", &entry.ieta, "ieta/I");
0137   tree->Branch("iphi", &entry.iphi, "iphi/I");
0138   tree->Branch("depth", &entry.depth, "depth/I");
0139   tree->Branch("samplingFactor", &entry.samplingFactor, "samplingFactor/D");
0140   tree->Branch("fCtoGeV", &entry.fCtoGeV, "fCtoGeV/D");
0141   tree->Branch("photoelectronsToAnalog", &entry.photoelectronsToAnalog, "photoelectronsToAnalog/D");
0142   tree->Branch("simHitToPhotoelectrons", &entry.simHitToPhotoelectrons, "simHitToPhotoelectrons/D");
0143   tree->Branch("energy", "vector<double>", &entry.energy, 32000, 0);
0144   tree->Branch("time", "vector<double>", &entry.time, 32000, 0);
0145   tree->Branch("tof", &entry.tof, "tof/D");
0146   tree->Branch("signalTot", "vector<double>", &entry.signalTot, 32000, 0);
0147   tree->Branch("signalTotPrecise", "vector<double>", &entry.signalTotPrecise, 32000, 0);
0148 }
0149
0150 void CaloSamplesAnalyzer::doBeginRun_(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0151   auto geometry = &iSetup.getData(tokGeom_);
0152   auto pHRNDC = &iSetup.getData(tokRecCns_);
0153
0154   // See if it's been updated
0155   if (geometry != theGeometry) {
0156     theGeometry = geometry;
0157     theRecNumber = pHRNDC;
0158     theResponse->setGeometry(theGeometry);
0159   }
0160 }
0161
0162 // ------------ method called on each new Event  ------------
0163 void CaloSamplesAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0164   treemap.clear();
0165
0166   //conditions
0167   auto conditions = &iSetup.getData(tokDB_);
0168   theParameterMap->setDbService(conditions);
0169
0170   // Event information
0171   const edm::EventAuxiliary& aux = iEvent.eventAuxiliary();
0172   unsigned run = aux.run();
0173   unsigned lumiblock = aux.luminosityBlock();
0174   unsigned long long event = aux.event();
0175
0176   edm::Handle<std::vector<CaloSamples>> h_CS;
0177   iEvent.getByToken(tok_calo, h_CS);
0178
0179   //loop over CaloSamples, extract info per channel
0180   for (const auto& iCS : *(h_CS.product())) {
0181     DetId did(iCS.id());
0182     if (did.det() != DetId::Hcal)
0183       continue;
0184     HcalDetId hid(did);
0185     CaloNtuple& ntup = treemap[hid.rawId()];
0186     //check for existence of entry
0187     if (ntup.id == 0) {
0188       ntup.id = hid.rawId();
0189       ntup.subdet = hid.subdet();
0190       ntup.ieta = hid.ieta();
0191       ntup.iphi = hid.iphi();
0192       ntup.depth = hid.depth();
0193       //get sim parameters
0194       if (hid.subdet() == HcalForward) {
0195         const HFSimParameters& pars = dynamic_cast<const HFSimParameters&>(theParameterMap->simParameters(hid));
0196         ntup.samplingFactor = pars.samplingFactor();
0197         ntup.fCtoGeV = pars.fCtoGeV(hid);
0198         ntup.photoelectronsToAnalog = pars.photoelectronsToAnalog(hid);
0199         ntup.simHitToPhotoelectrons = pars.simHitToPhotoelectrons(hid);
0200       } else {
0201         const HcalSimParameters& pars = dynamic_cast<const HcalSimParameters&>(theParameterMap->simParameters(hid));
0202         ntup.samplingFactor = pars.samplingFactor(hid);
0203         ntup.fCtoGeV = pars.fCtoGeV(hid);
0204         ntup.photoelectronsToAnalog = pars.photoelectronsToAnalog(hid);
0205         ntup.simHitToPhotoelectrons = pars.simHitToPhotoelectrons(hid);
0206       }
0207       //get time of flight
0208       ntup.tof = theResponse->timeOfFlight(hid);
0209     }
0210     //get pulse info every time
0211     if (ntup.signalTot.empty())
0212       ntup.signalTot.resize(iCS.size(), 0.0);
0213     for (int i = 0; i < iCS.size(); ++i) {
0214       ntup.signalTot[i] += iCS[i];
0215     }
0216     if (ntup.signalTotPrecise.empty())
0217       ntup.signalTotPrecise.resize(iCS.preciseSize(), 0.0);
0218     for (int i = 0; i < iCS.preciseSize(); ++i) {
0219       ntup.signalTotPrecise[i] += iCS.preciseAt(i);
0220     }
0221   }
0222
0223   edm::Handle<std::vector<PCaloHit>> h_SH;
0224   iEvent.getByToken(tok_sim, h_SH);
0225
0226   //loop over simhits, extract info per channel
0227   for (const auto& iSH : *(h_SH.product())) {
0228     HcalDetId hid;
0229     unsigned int id = iSH.id();
0230     if (TestNumbering)
0231       hid = HcalDetId(HcalHitRelabeller::relabel(id, theRecNumber));
0232     else
0233       hid = HcalDetId(id);
0234     auto ntupIt = treemap.find(hid.rawId());
0235     if (ntupIt == treemap.end())
0236       continue;
0237     CaloNtuple& ntup = ntupIt->second;
0238     //append simhit info
0239     ntup.energy.push_back(iSH.energy());
0240     ntup.time.push_back(iSH.time());
0241   }
0242
0243   //one tree entry per map entry
0244   entry.run = run;
0245   entry.lumiblock = lumiblock;
0246   entry.event = event;
0247   for (const auto& ntup : treemap) {
0248     entry.id = ntup.second.id;
0249     entry.subdet = ntup.second.subdet;
0250     entry.ieta = ntup.second.ieta;
0251     entry.iphi = ntup.second.iphi;
0252     entry.depth = ntup.second.depth;
0253     entry.samplingFactor = ntup.second.samplingFactor;
0254     entry.fCtoGeV = ntup.second.fCtoGeV;
0255     entry.photoelectronsToAnalog = ntup.second.photoelectronsToAnalog;
0256     entry.simHitToPhotoelectrons = ntup.second.simHitToPhotoelectrons;
0257     entry.tof = ntup.second.tof;
0258     entry.energy = ntup.second.energy;
0259     entry.time = ntup.second.time;
0260     entry.signalTot = ntup.second.signalTot;
0261     entry.signalTotPrecise = ntup.second.signalTotPrecise;
0262     tree->Fill();
0263   }
0264 }
0265
0266 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0267 void CaloSamplesAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0268   //The following says we do not know what parameters are allowed so do no validation
0269   // Please change this to state exactly what you do use, even if it is no parameters
0270   edm::ParameterSetDescription desc;
0271   desc.setUnknown();
0272   descriptions.addDefault(desc);
0273 }
0274 //define this as a plug-in
0275 DEFINE_FWK_MODULE(CaloSamplesAnalyzer);