File indexing completed on 2021-11-26 23:20:34
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <utility>
0021 #include <algorithm>
0022 #include <memory>
0023 #include <vector>
0024 #include <iostream>
0025 #include <fstream>
0026
0027 #include "TTree.h"
0028 #include "TFile.h"
0029 #include "TH1D.h"
0030
0031
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include "FWCore/Framework/interface/MakerMacros.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037
0038 #include "FWCore/ServiceRegistry/interface/Service.h"
0039 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0040
0041 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0042 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0043
0044 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0045 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0046 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0047 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0048 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
0049
0050 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0051 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0052 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0053
0054 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0055 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0056 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
0057
0058
0059 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentManager.h"
0060 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCorrectionFunctions.h"
0061
0062 #include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFitOOTPileupCorrection.h"
0063 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalDeterministicFit.h"
0064 #include "RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h"
0065 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0066
0067 #include "RecoLocalCalo/HcalRecAlgos/interface/parseHBHEPhase1AlgoDescription.h"
0068
0069
0070
0071
0072
0073 class MahiDebugger : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0074 public:
0075 explicit MahiDebugger(const edm::ParameterSet&);
0076 ~MahiDebugger() override;
0077
0078 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0079
0080 private:
0081 void beginJob() override;
0082 void analyze(const edm::Event&, const edm::EventSetup&) override;
0083
0084 float hbminusCorrectionFactor(const HcalDetId& cell, int runnum, float energy, bool isRealData) const;
0085 void endJob() override;
0086
0087
0088
0089
0090 bool dynamicPed_;
0091 float ts4Thresh_;
0092 float chiSqSwitch_;
0093
0094 bool applyTimeSlew_;
0095 HcalTimeSlew::BiasSetting slewFlavor_;
0096 double tsDelay1GeV_ = 0;
0097
0098 bool calculateArrivalTime_;
0099 float meanTime_;
0100 float timeSigmaHPD_;
0101 float timeSigmaSiPM_;
0102
0103 std::vector<int> activeBXs_;
0104
0105 int nMaxItersMin_;
0106 int nMaxItersNNLS_;
0107
0108 float deltaChiSqThresh_;
0109 float nnlsThresh_;
0110
0111 unsigned int bxSizeConf_;
0112 int bxOffsetConf_;
0113
0114
0115 HcalPulseShapes theHcalPulseShapes_;
0116 int cntsetPulseShape_;
0117 std::unique_ptr<FitterFuncs::PulseShapeFunctor> psfPtr_;
0118 std::unique_ptr<ROOT::Math::Functor> pfunctor_;
0119
0120 std::unique_ptr<MahiFit> mahi_;
0121
0122 edm::EDGetTokenT<HBHEChannelInfoCollection> token_ChannelInfo_;
0123 const edm::ESGetToken<HcalTimeSlew, HcalTimeSlewRecord> tokDelay_;
0124
0125 const HcalTimeSlew* hcalTimeSlewDelay;
0126
0127 edm::Service<TFileService> FileService;
0128 TTree* outTree;
0129
0130 int run;
0131 int evt;
0132 int ls;
0133
0134 int nBxTrain;
0135
0136 int ieta;
0137 int iphi;
0138 int depth;
0139
0140 int nSamples;
0141 int soi;
0142
0143 bool use8;
0144
0145 float inTimeConst;
0146 float inDarkCurrent;
0147 float inPedAvg;
0148 float inGain;
0149
0150 float inNoiseADC[10];
0151 float inNoiseDC[10];
0152 float inNoisePhoto[10];
0153 float inPedestal[10];
0154
0155 float totalUCNoise[10];
0156
0157 float mahiEnergy;
0158 float chiSq;
0159 float arrivalTime;
0160
0161 float ootEnergy[7];
0162 float pedEnergy;
0163
0164 float count[10];
0165 float inputTS[10];
0166 int inputTDC[10];
0167 float itPulse[10];
0168 float ootPulse[7][10];
0169 };
0170
0171 MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig)
0172 : dynamicPed_(iConfig.getParameter<bool>("dynamicPed")),
0173 ts4Thresh_(iConfig.getParameter<double>("ts4Thresh")),
0174 chiSqSwitch_(iConfig.getParameter<double>("chiSqSwitch")),
0175 applyTimeSlew_(iConfig.getParameter<bool>("applyTimeSlew")),
0176 calculateArrivalTime_(iConfig.getParameter<bool>("calculateArrivalTime")),
0177 meanTime_(iConfig.getParameter<double>("meanTime")),
0178 timeSigmaHPD_(iConfig.getParameter<double>("timeSigmaHPD")),
0179 timeSigmaSiPM_(iConfig.getParameter<double>("timeSigmaSiPM")),
0180 activeBXs_(iConfig.getParameter<std::vector<int>>("activeBXs")),
0181 nMaxItersMin_(iConfig.getParameter<int>("nMaxItersMin")),
0182 nMaxItersNNLS_(iConfig.getParameter<int>("nMaxItersNNLS")),
0183 deltaChiSqThresh_(iConfig.getParameter<double>("deltaChiSqThresh")),
0184 nnlsThresh_(iConfig.getParameter<double>("nnlsThresh")),
0185 tokDelay_(esConsumes<HcalTimeSlew, HcalTimeSlewRecord>(edm::ESInputTag("", "HBHE"))) {
0186 usesResource("TFileService");
0187
0188 mahi_ = std::make_unique<MahiFit>();
0189
0190 mahi_->setParameters(dynamicPed_,
0191 ts4Thresh_,
0192 chiSqSwitch_,
0193 applyTimeSlew_,
0194 HcalTimeSlew::Medium,
0195 calculateArrivalTime_,
0196 meanTime_,
0197 timeSigmaHPD_,
0198 timeSigmaSiPM_,
0199 activeBXs_,
0200 nMaxItersMin_,
0201 nMaxItersNNLS_,
0202 deltaChiSqThresh_,
0203 nnlsThresh_);
0204
0205 token_ChannelInfo_ = consumes<HBHEChannelInfoCollection>(iConfig.getParameter<edm::InputTag>("recoLabel"));
0206 }
0207
0208 MahiDebugger::~MahiDebugger() {}
0209
0210
0211
0212
0213
0214
0215 void MahiDebugger::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0216 using namespace edm;
0217
0218 hcalTimeSlewDelay = &iSetup.getData(tokDelay_);
0219
0220 run = iEvent.id().run();
0221 evt = iEvent.id().event();
0222 ls = iEvent.id().luminosityBlock();
0223
0224 edm::EventBase const& eventbase = iEvent;
0225 nBxTrain = int(eventbase.bunchCrossing());
0226
0227 Handle<HBHEChannelInfoCollection> hChannelInfo;
0228 iEvent.getByToken(token_ChannelInfo_, hChannelInfo);
0229
0230 for (HBHEChannelInfoCollection::const_iterator iter = hChannelInfo->begin(); iter != hChannelInfo->end(); iter++) {
0231 const HBHEChannelInfo& hci(*iter);
0232 const HcalDetId detid = hci.id();
0233
0234 ieta = detid.ieta();
0235 iphi = detid.iphi();
0236 depth = detid.depth();
0237
0238 const bool isRealData = true;
0239
0240 const MahiFit* mahi = mahi_.get();
0241 mahi_->setPulseShapeTemplate(
0242 hci.recoShape(), theHcalPulseShapes_, hci.hasTimeInfo(), hcalTimeSlewDelay, hci.nSamples());
0243 MahiDebugInfo mdi;
0244
0245 mdi.mahiEnergy = 0;
0246 for (unsigned int ioot = 0; ioot < 7; ioot++)
0247 mdi.ootEnergy[ioot] = 0;
0248 mahi->phase1Debug(hci, mdi);
0249
0250 nSamples = mdi.nSamples;
0251 soi = mdi.soi;
0252
0253 inTimeConst = mdi.inTimeConst;
0254 inDarkCurrent = mdi.inDarkCurrent;
0255 inPedAvg = mdi.inPedAvg;
0256 inGain = mdi.inGain;
0257
0258 use8 = mdi.use3;
0259 chiSq = mdi.chiSq;
0260 arrivalTime = mdi.arrivalTime;
0261 mahiEnergy = mdi.mahiEnergy;
0262 mahiEnergy *= hbminusCorrectionFactor(detid, run, mahiEnergy, isRealData);
0263 for (unsigned int ioot = 0; ioot < 7; ioot++) {
0264 ootEnergy[ioot] = mdi.ootEnergy[ioot];
0265 ootEnergy[ioot] *= hbminusCorrectionFactor(detid, run, ootEnergy[ioot], isRealData);
0266 }
0267 pedEnergy = mdi.pedEnergy;
0268 pedEnergy *= hbminusCorrectionFactor(detid, run, pedEnergy, isRealData);
0269
0270 for (int i = 0; i < nSamples; i++) {
0271 count[i] = mdi.count[i];
0272 inputTS[i] = mdi.inputTS[i];
0273 inputTDC[i] = mdi.inputTDC[i];
0274 itPulse[i] = mdi.itPulse[i];
0275 for (unsigned int ioot = 0; ioot < 7; ioot++)
0276 ootPulse[ioot][i] = mdi.ootPulse[ioot][i];
0277
0278 inNoiseADC[i] = mdi.inNoiseADC[i];
0279 inNoiseDC[i] = mdi.inNoiseDC[i];
0280 inNoisePhoto[i] = mdi.inNoisePhoto[i];
0281 inPedestal[i] = mdi.inPedestal[i];
0282 totalUCNoise[i] = mdi.totalUCNoise[i];
0283 }
0284 if (nSamples == 8) {
0285 count[8] = 8;
0286 count[9] = 9;
0287 }
0288
0289 outTree->Fill();
0290 }
0291 }
0292
0293 float MahiDebugger::hbminusCorrectionFactor(const HcalDetId& cell,
0294 int runnum,
0295 const float energy,
0296 const bool isRealData) const {
0297 float corr = 1.f;
0298 if (isRealData && runnum > 0)
0299 if (cell.subdet() == HcalBarrel) {
0300 const int ieta = cell.ieta();
0301 const int iphi = cell.iphi();
0302 corr = hbminus_special_ecorr(ieta, iphi, energy, runnum);
0303 }
0304 return corr;
0305 }
0306
0307
0308 void MahiDebugger::beginJob() {
0309 outTree = FileService->make<TTree>("HcalTree", "HcalTree");
0310
0311 outTree->Branch("run", &run, "run/I");
0312 outTree->Branch("evt", &evt, "evt/I");
0313 outTree->Branch("ls", &ls, "ls/I");
0314 outTree->Branch("nBxTrain", &nBxTrain, "nBxTrain/I");
0315
0316 outTree->Branch("ieta", &ieta, "ieta/I");
0317 outTree->Branch("iphi", &iphi, "iphi/I");
0318 outTree->Branch("depth", &depth, "depth/I");
0319 outTree->Branch("nSamples", &nSamples, "nSamples/I");
0320 outTree->Branch("soi", &soi, "soi/I");
0321
0322 outTree->Branch("inTimeConst", &inTimeConst, "inTimeConst/F");
0323 outTree->Branch("inDarkCurrent", &inDarkCurrent, "inDarkCurrent/F");
0324 outTree->Branch("inPedAvg", &inPedAvg, "inPedAvg/F");
0325 outTree->Branch("inGain", &inGain, "inGain/F");
0326
0327 outTree->Branch("use8", &use8, "use8/B");
0328 outTree->Branch("mahiEnergy", &mahiEnergy, "mahiEnergy/F");
0329 outTree->Branch("chiSq", &chiSq, "chiSq/F");
0330 outTree->Branch("arrivalTime", &arrivalTime, "arrivalTime/F");
0331 outTree->Branch("ootEnergy", &ootEnergy, "ootEnergy[7]/F");
0332 outTree->Branch("pedEnergy", &pedEnergy, "pedEnergy/F");
0333 outTree->Branch("count", &count, "count[10]/F");
0334 outTree->Branch("inputTS", &inputTS, "inputTS[10]/F");
0335 outTree->Branch("inputTDC", &inputTDC, "inputTDC[10]/I");
0336 outTree->Branch("itPulse", &itPulse, "itPulse[10]/F");
0337 outTree->Branch("ootPulse", &ootPulse, "ootPulse[7][10]/F");
0338
0339 outTree->Branch("inNoiseADC", &inNoiseADC, "inNoiseADC[10]/F");
0340 outTree->Branch("inNoiseDC", &inNoiseDC, "inNoiseDC[10]/F");
0341 outTree->Branch("inNoisePhoto", &inNoisePhoto, "inNoisePhoto[10]/F");
0342 outTree->Branch("inPedestal", &inPedestal, "inPedestal[10]/F");
0343 outTree->Branch("totalUCNoise", &totalUCNoise, "totalUCNoise[10]/F");
0344 }
0345
0346
0347 void MahiDebugger::endJob() {}
0348
0349 #define add_param_set(name) \
0350 edm::ParameterSetDescription name; \
0351 name.setAllowAnything(); \
0352 desc.add<edm::ParameterSetDescription>(#name, name)
0353
0354
0355 void MahiDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0356 edm::ParameterSetDescription desc;
0357
0358 desc.add<edm::InputTag>("recoLabel");
0359 desc.add<bool>("dynamicPed");
0360 desc.add<bool>("calculateArrivalTime");
0361 desc.add<double>("ts4Thresh");
0362 desc.add<double>("chiSqSwitch");
0363 desc.add<bool>("applyTimeSlew");
0364 desc.add<double>("meanTime");
0365 desc.add<double>("timeSigmaHPD");
0366 desc.add<double>("timeSigmaSiPM");
0367 desc.add<std::vector<int>>("activeBXs");
0368 desc.add<int>("nMaxItersMin");
0369 desc.add<int>("nMaxItersNNLS");
0370 desc.add<double>("deltaChiSqThresh");
0371 desc.add<double>("nnlsThresh");
0372
0373
0374
0375 descriptions.addDefault(desc);
0376 }
0377
0378
0379 DEFINE_FWK_MODULE(MahiDebugger);