File indexing completed on 2024-04-06 12:29:32
0001
0002 #include <memory>
0003
0004
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
0013 #include <TTree.h>
0014
0015
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
0038 #include <vector>
0039 #include <unordered_map>
0040 #include <string>
0041
0042
0043
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
0081 edm::Service<TFileService> fs;
0082 TTree* tree;
0083 const CaloGeometry* theGeometry;
0084 const HcalDDDRecConstants* theRecNumber;
0085 CaloHitResponse* theResponse;
0086 HcalSimParameterMap* theParameterMap;
0087
0088 CaloNtuple entry;
0089 std::unordered_map<unsigned, CaloNtuple> treemap;
0090
0091 bool TestNumbering;
0092
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
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
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
0155 if (geometry != theGeometry) {
0156 theGeometry = geometry;
0157 theRecNumber = pHRNDC;
0158 theResponse->setGeometry(theGeometry);
0159 }
0160 }
0161
0162
0163 void CaloSamplesAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0164 treemap.clear();
0165
0166
0167 auto conditions = &iSetup.getData(tokDB_);
0168 theParameterMap->setDbService(conditions);
0169
0170
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
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
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
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
0208 ntup.tof = theResponse->timeOfFlight(hid);
0209 }
0210
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
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
0239 ntup.energy.push_back(iSH.energy());
0240 ntup.time.push_back(iSH.time());
0241 }
0242
0243
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
0267 void CaloSamplesAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0268
0269
0270 edm::ParameterSetDescription desc;
0271 desc.setUnknown();
0272 descriptions.addDefault(desc);
0273 }
0274
0275 DEFINE_FWK_MODULE(CaloSamplesAnalyzer);