HcalRecAlgos/test/MahiDebugger.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoLocalCalo/HcalRecAlgos
0004 // Class:      MahiDebugger
0005 //
0006 /**\class MahiDebugger MahiDebugger.cc RecoLocalCalo/HcalRecAlgos/plugins/MahiDebugger.cc
0007
0008  Description: Tool to extract and store debugging information from the HBHE Reconstruction algorithm Mahi
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Jay Mathew Lawhorn
0015 //         Created:  Sat, 10 Feb 2018 10:02:38 GMT
0016 //
0017 //
0018
0019 // system include files
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 // user include files
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 //#include "RecoLocalCalo/HcalRecAlgos/test/SimpleHBHEPhase1AlgoDebug.h"
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 // class declaration
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   // Special HB- correction
0084   float hbminusCorrectionFactor(const HcalDetId& cell, int runnum, float energy, bool isRealData) const;
0085   void endJob() override;
0086
0087   // ----------member data ---------------------------
0088
0089   // Python-configurables
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   int timeAlgo_;
0100   float thEnergeticPulses_;
0101   float meanTime_;
0102   float timeSigmaHPD_;
0103   float timeSigmaSiPM_;
0104
0105   std::vector<int> activeBXs_;
0106
0107   int nMaxItersMin_;
0108   int nMaxItersNNLS_;
0109
0110   float deltaChiSqThresh_;
0111   float nnlsThresh_;
0112
0113   unsigned int bxSizeConf_;
0114   int bxOffsetConf_;
0115
0116   //for pulse shapes
0117   HcalPulseShapes theHcalPulseShapes_;
0118   int cntsetPulseShape_;
0119   std::unique_ptr<FitterFuncs::PulseShapeFunctor> psfPtr_;
0120   std::unique_ptr<ROOT::Math::Functor> pfunctor_;
0121
0122   std::unique_ptr<MahiFit> mahi_;
0123
0124   edm::EDGetTokenT<HBHEChannelInfoCollection> token_ChannelInfo_;
0125   const edm::ESGetToken<HcalTimeSlew, HcalTimeSlewRecord> tokDelay_;
0126
0127   const HcalTimeSlew* hcalTimeSlewDelay;
0128
0129   edm::Service<TFileService> FileService;
0130   TTree* outTree;
0131
0132   int run;
0133   int evt;
0134   int ls;
0135
0136   int nBxTrain;
0137
0138   int ieta;
0139   int iphi;
0140   int depth;
0141
0142   int nSamples;
0143   int soi;
0144
0145   bool use8;
0146
0147   float inTimeConst;
0148   float inDarkCurrent;
0149   float inPedAvg;
0150   float inGain;
0151
0152   float inNoiseADC[10];
0153   float inNoiseDC[10];
0154   float inNoisePhoto[10];
0155   float inPedestal[10];
0156
0157   float totalUCNoise[10];
0158
0159   float mahiEnergy;  //SOI charge
0160   float chiSq;
0161   float arrivalTime;
0162
0163   float ootEnergy[7];  //OOT charge
0164   float pedEnergy;     //pedestal charge
0165
0166   float count[10];        //TS value 0-9
0167   float inputTS[10];      //input TS samples
0168   int inputTDC[10];       //input TS samples
0169   float itPulse[10];      //SOI pulse shape
0170   float ootPulse[7][10];  //OOT pulse shape
0171 };
0172
0173 MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig)
0174     : dynamicPed_(iConfig.getParameter<bool>("dynamicPed")),
0175       ts4Thresh_(iConfig.getParameter<double>("ts4Thresh")),
0176       chiSqSwitch_(iConfig.getParameter<double>("chiSqSwitch")),
0177       applyTimeSlew_(iConfig.getParameter<bool>("applyTimeSlew")),
0178       calculateArrivalTime_(iConfig.getParameter<bool>("calculateArrivalTime")),
0179       timeAlgo_(iConfig.getParameter<int>("timeAlgo")),
0180       thEnergeticPulses_(iConfig.getParameter<double>("thEnergeticPulses")),
0181       meanTime_(iConfig.getParameter<double>("meanTime")),
0182       timeSigmaHPD_(iConfig.getParameter<double>("timeSigmaHPD")),
0183       timeSigmaSiPM_(iConfig.getParameter<double>("timeSigmaSiPM")),
0184       activeBXs_(iConfig.getParameter<std::vector<int>>("activeBXs")),
0185       nMaxItersMin_(iConfig.getParameter<int>("nMaxItersMin")),
0186       nMaxItersNNLS_(iConfig.getParameter<int>("nMaxItersNNLS")),
0187       deltaChiSqThresh_(iConfig.getParameter<double>("deltaChiSqThresh")),
0188       nnlsThresh_(iConfig.getParameter<double>("nnlsThresh")),
0189       tokDelay_(esConsumes<HcalTimeSlew, HcalTimeSlewRecord>(edm::ESInputTag("", "HBHE"))) {
0190   usesResource("TFileService");
0191
0192   mahi_ = std::make_unique<MahiFit>();
0193
0194   mahi_->setParameters(dynamicPed_,
0195                        ts4Thresh_,
0196                        chiSqSwitch_,
0197                        applyTimeSlew_,
0198                        HcalTimeSlew::Medium,
0199                        calculateArrivalTime_,
0200                        timeAlgo_,
0201                        thEnergeticPulses_,
0202                        meanTime_,
0203                        timeSigmaHPD_,
0204                        timeSigmaSiPM_,
0205                        activeBXs_,
0206                        nMaxItersMin_,
0207                        nMaxItersNNLS_,
0208                        deltaChiSqThresh_,
0209                        nnlsThresh_);
0210
0211   token_ChannelInfo_ = consumes<HBHEChannelInfoCollection>(iConfig.getParameter<edm::InputTag>("recoLabel"));
0212 }
0213
0214 MahiDebugger::~MahiDebugger() {}
0215
0216 //
0217 // member functions
0218 //
0219
0220 // ------------ method called for each event  ------------
0221 void MahiDebugger::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0222   using namespace edm;
0223
0224   hcalTimeSlewDelay = &iSetup.getData(tokDelay_);
0225
0226   run = iEvent.id().run();
0227   evt = iEvent.id().event();
0228   ls = iEvent.id().luminosityBlock();
0229
0230   edm::EventBase const& eventbase = iEvent;
0231   nBxTrain = int(eventbase.bunchCrossing());
0232
0233   Handle<HBHEChannelInfoCollection> hChannelInfo;
0234   iEvent.getByToken(token_ChannelInfo_, hChannelInfo);
0235
0236   for (HBHEChannelInfoCollection::const_iterator iter = hChannelInfo->begin(); iter != hChannelInfo->end(); iter++) {
0237     const HBHEChannelInfo& hci(*iter);
0238     const HcalDetId detid = hci.id();
0239
0240     ieta = detid.ieta();
0241     iphi = detid.iphi();
0242     depth = detid.depth();
0243
0244     const bool isRealData = true;
0245
0246     const MahiFit* mahi = mahi_.get();
0247     mahi_->setPulseShapeTemplate(
0248         hci.recoShape(), theHcalPulseShapes_, hci.hasTimeInfo(), hcalTimeSlewDelay, hci.nSamples(), hci.tsGain(0));
0249     MahiDebugInfo mdi;
0250     // initialize energies so that the values in the previous iteration are not stored
0251     mdi.mahiEnergy = 0;
0252     for (unsigned int ioot = 0; ioot < 7; ioot++)
0253       mdi.ootEnergy[ioot] = 0;
0254     mahi->phase1Debug(hci, mdi);
0255
0256     nSamples = mdi.nSamples;
0257     soi = mdi.soi;
0258
0259     inTimeConst = mdi.inTimeConst;
0260     inDarkCurrent = mdi.inDarkCurrent;
0261     inPedAvg = mdi.inPedAvg;
0262     inGain = mdi.inGain;
0263
0264     use8 = mdi.use3;
0265     chiSq = mdi.chiSq;
0266     arrivalTime = mdi.arrivalTime;
0267     mahiEnergy = mdi.mahiEnergy;
0268     mahiEnergy *= hbminusCorrectionFactor(detid, run, mahiEnergy, isRealData);
0269     for (unsigned int ioot = 0; ioot < 7; ioot++) {
0270       ootEnergy[ioot] = mdi.ootEnergy[ioot];
0271       ootEnergy[ioot] *= hbminusCorrectionFactor(detid, run, ootEnergy[ioot], isRealData);
0272     }
0273     pedEnergy = mdi.pedEnergy;
0274     pedEnergy *= hbminusCorrectionFactor(detid, run, pedEnergy, isRealData);
0275
0276     for (int i = 0; i < nSamples; i++) {
0277       count[i] = mdi.count[i];
0278       inputTS[i] = mdi.inputTS[i];
0279       inputTDC[i] = mdi.inputTDC[i];
0280       itPulse[i] = mdi.itPulse[i];
0281       for (unsigned int ioot = 0; ioot < 7; ioot++)
0282         ootPulse[ioot][i] = mdi.ootPulse[ioot][i];
0283
0284       inNoiseADC[i] = mdi.inNoiseADC[i];
0285       inNoiseDC[i] = mdi.inNoiseDC[i];
0286       inNoisePhoto[i] = mdi.inNoisePhoto[i];
0287       inPedestal[i] = mdi.inPedestal[i];
0288       totalUCNoise[i] = mdi.totalUCNoise[i];
0289     }
0290     if (nSamples == 8) {
0291       count[8] = 8;
0292       count[9] = 9;
0293     }
0294
0295     outTree->Fill();
0296   }
0297 }
0298
0299 float MahiDebugger::hbminusCorrectionFactor(const HcalDetId& cell,
0300                                             int runnum,
0301                                             const float energy,
0302                                             const bool isRealData) const {
0303   float corr = 1.f;
0304   if (isRealData && runnum > 0)
0305     if (cell.subdet() == HcalBarrel) {
0306       const int ieta = cell.ieta();
0307       const int iphi = cell.iphi();
0308       corr = hbminus_special_ecorr(ieta, iphi, energy, runnum);
0309     }
0310   return corr;
0311 }
0312
0313 // ------------ method called once each job just before starting event loop  ------------
0314 void MahiDebugger::beginJob() {
0315   outTree = FileService->make<TTree>("HcalTree", "HcalTree");
0316
0317   outTree->Branch("run", &run, "run/I");
0318   outTree->Branch("evt", &evt, "evt/I");
0319   outTree->Branch("ls", &ls, "ls/I");
0320   outTree->Branch("nBxTrain", &nBxTrain, "nBxTrain/I");
0321
0322   outTree->Branch("ieta", &ieta, "ieta/I");
0323   outTree->Branch("iphi", &iphi, "iphi/I");
0324   outTree->Branch("depth", &depth, "depth/I");
0325   outTree->Branch("nSamples", &nSamples, "nSamples/I");
0326   outTree->Branch("soi", &soi, "soi/I");
0327
0328   outTree->Branch("inTimeConst", &inTimeConst, "inTimeConst/F");
0329   outTree->Branch("inDarkCurrent", &inDarkCurrent, "inDarkCurrent/F");
0330   outTree->Branch("inPedAvg", &inPedAvg, "inPedAvg/F");
0331   outTree->Branch("inGain", &inGain, "inGain/F");
0332
0333   outTree->Branch("use8", &use8, "use8/B");
0334   outTree->Branch("mahiEnergy", &mahiEnergy, "mahiEnergy/F");
0335   outTree->Branch("chiSq", &chiSq, "chiSq/F");
0336   outTree->Branch("arrivalTime", &arrivalTime, "arrivalTime/F");
0337   outTree->Branch("ootEnergy", &ootEnergy, "ootEnergy[7]/F");
0338   outTree->Branch("pedEnergy", &pedEnergy, "pedEnergy/F");
0339   outTree->Branch("count", &count, "count[10]/F");
0340   outTree->Branch("inputTS", &inputTS, "inputTS[10]/F");
0341   outTree->Branch("inputTDC", &inputTDC, "inputTDC[10]/I");
0342   outTree->Branch("itPulse", &itPulse, "itPulse[10]/F");
0343   outTree->Branch("ootPulse", &ootPulse, "ootPulse[7][10]/F");
0344
0345   outTree->Branch("inNoiseADC", &inNoiseADC, "inNoiseADC[10]/F");
0346   outTree->Branch("inNoiseDC", &inNoiseDC, "inNoiseDC[10]/F");
0347   outTree->Branch("inNoisePhoto", &inNoisePhoto, "inNoisePhoto[10]/F");
0348   outTree->Branch("inPedestal", &inPedestal, "inPedestal[10]/F");
0349   outTree->Branch("totalUCNoise", &totalUCNoise, "totalUCNoise[10]/F");
0350 }
0351
0352 // ------------ method called once each job just after ending the event loop  ------------
0353 void MahiDebugger::endJob() {}
0354
0355 #define add_param_set(name)          \
0356   edm::ParameterSetDescription name; \
0357   name.setAllowAnything();           \
0358   desc.add<edm::ParameterSetDescription>(#name, name)
0359
0360 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0361 void MahiDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0362   edm::ParameterSetDescription desc;
0363
0364   desc.add<edm::InputTag>("recoLabel");
0365   desc.add<bool>("dynamicPed");
0366   desc.add<bool>("calculateArrivalTime");
0367   desc.add<int>("timeAlgo");
0368   desc.add<double>("thEnergeticPulse");
0369   desc.add<double>("ts4Thresh");
0370   desc.add<double>("chiSqSwitch");
0371   desc.add<bool>("applyTimeSlew");
0372   desc.add<double>("meanTime");
0373   desc.add<double>("timeSigmaHPD");
0374   desc.add<double>("timeSigmaSiPM");
0375   desc.add<std::vector<int>>("activeBXs");
0376   desc.add<int>("nMaxItersMin");
0377   desc.add<int>("nMaxItersNNLS");
0378   desc.add<double>("deltaChiSqThresh");
0379   desc.add<double>("nnlsThresh");
0380
0381   //desc.add<std::string>("algoConfigClass");
0382   //add_param_set(algorithm);
0383   descriptions.addDefault(desc);
0384 }
0385
0386 //define this as a plug-in
0387 DEFINE_FWK_MODULE(MahiDebugger);