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   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   //for pulse shapes
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;  //SOI charge
0158   float chiSq;
0159   float arrivalTime;
0160
0161   float ootEnergy[7];  //OOT charge
0162   float pedEnergy;     //pedestal charge
0163
0164   float count[10];        //TS value 0-9
0165   float inputTS[10];      //input TS samples
0166   int inputTDC[10];       //input TS samples
0167   float itPulse[10];      //SOI pulse shape
0168   float ootPulse[7][10];  //OOT pulse shape
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 // member functions
0212 //
0213
0214 // ------------ method called for each event  ------------
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     // initialize energies so that the values in the previous iteration are not stored
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 // ------------ method called once each job just before starting event loop  ------------
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 // ------------ method called once each job just after ending the event loop  ------------
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 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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   //desc.add<std::string>("algoConfigClass");
0374   //add_param_set(algorithm);
0375   descriptions.addDefault(desc);
0376 }
0377
0378 //define this as a plug-in
0379 DEFINE_FWK_MODULE(MahiDebugger);