HGCalSimAlgos/test/HGCSiNoiseMapAnalyzer.cc

0001 // user include files
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/StreamID.h"
0010
0011 // Geometry
0012 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0013 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0016 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0017
0018 #include "SimCalorimetry/HGCalSimAlgos/interface/HGCalSiNoiseMap.h"
0019
0020 #include <TH1F.h>
0021 #include <TH2F.h>
0022 #include <TProfile.h>
0023
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/PluginManager/interface/ModuleDef.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0028
0029 //STL headers
0030 #include <vector>
0031 #include <sstream>
0032 #include <string>
0033 #include <iomanip>
0034
0035 #include "vdt/vdtMath.h"
0036
0037 using namespace std;
0038
0039 //
0040 // class declaration
0041 //
0042 class HGCSiNoiseMapAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0043 public:
0044   explicit HGCSiNoiseMapAnalyzer(const edm::ParameterSet &);
0045   ~HGCSiNoiseMapAnalyzer() override;
0046
0047 private:
0048   void beginJob() override {}
0049   void analyze(const edm::Event &, const edm::EventSetup &) override;
0050   void endJob() override {}
0051
0052   // ----------member data ---------------------------
0053   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokGeom_;
0054   edm::Service<TFileService> fs_;
0055   std::map<DetId::Detector, std::unique_ptr<HGCalSiNoiseMap<HGCSiliconDetId>>> noiseMaps_;
0056   std::map<std::pair<DetId::Detector, int>, TH1F *> layerN_, layerCCE_, layerNoise_, layerIleak_, layerSN_, layerF_,
0057       layerGain_, layerMipPeak_;
0058   std::map<DetId::Detector, TH2F *> detN_, detCCE_, detNoise_, detIleak_, detSN_, detF_, detGain_, detMipPeak_;
0059   std::map<std::pair<DetId::Detector, int>, TProfile *> detCCEVsFluence_;
0060
0061   int aimMIPtoADC_;
0062   bool ignoreGainSettings_;
0063
0064   const int plotMargin_ = 20;
0065 };
0066
0067 //
0068 HGCSiNoiseMapAnalyzer::HGCSiNoiseMapAnalyzer(const edm::ParameterSet &iConfig)
0069     : tokGeom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0070   usesResource("TFileService");
0071   fs_->file().cd();
0072
0073   //configure the dose map
0074   std::string doseMapURL(iConfig.getParameter<std::string>("doseMap"));
0075   unsigned int doseMapAlgo(iConfig.getParameter<unsigned int>("doseMapAlgo"));
0076   double scaleByDoseFactor(iConfig.getParameter<double>("scaleByDoseFactor"));
0077   std::vector<double> ileakParam(
0078       iConfig.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
0079   std::vector<double> cceParamFine(
0080       iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"));
0081   std::vector<double> cceParamThin(
0082       iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"));
0083   std::vector<double> cceParamThick(
0084       iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));
0085
0086   noiseMaps_[DetId::HGCalEE] = std::unique_ptr<HGCalSiNoiseMap<HGCSiliconDetId>>(new HGCalSiNoiseMap<HGCSiliconDetId>);
0087   noiseMaps_[DetId::HGCalEE]->setDoseMap(doseMapURL, doseMapAlgo);
0088   noiseMaps_[DetId::HGCalEE]->setFluenceScaleFactor(scaleByDoseFactor);
0089
0090   noiseMaps_[DetId::HGCalEE]->setIleakParam(ileakParam);
0091   noiseMaps_[DetId::HGCalEE]->setCceParam(cceParamFine, cceParamThin, cceParamThick);
0092
0093   noiseMaps_[DetId::HGCalHSi] = std::unique_ptr<HGCalSiNoiseMap<HGCSiliconDetId>>(new HGCalSiNoiseMap<HGCSiliconDetId>);
0094   noiseMaps_[DetId::HGCalHSi]->setDoseMap(doseMapURL, doseMapAlgo);
0095   noiseMaps_[DetId::HGCalHSi]->setFluenceScaleFactor(scaleByDoseFactor);
0096   noiseMaps_[DetId::HGCalHSi]->setIleakParam(ileakParam);
0097   noiseMaps_[DetId::HGCalHSi]->setCceParam(cceParamFine, cceParamThin, cceParamThick);
0098
0099   aimMIPtoADC_ = iConfig.getParameter<int>("aimMIPtoADC");
0100   ignoreGainSettings_ = iConfig.getParameter<bool>("ignoreGainSettings");
0101 }
0102
0103 //
0104 HGCSiNoiseMapAnalyzer::~HGCSiNoiseMapAnalyzer() {}
0105
0106 //
0107 void HGCSiNoiseMapAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &es) {
0108   //get geometry
0109   const auto &geom = es.getHandle(tokGeom_);
0110
0111   std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi};
0112   for (const auto &d : dets) {
0113     noiseMaps_[d]->setGeometry(geom->getSubdetectorGeometry(d, ForwardSubdetector::ForwardEmpty));
0114     //sub-detector boundaries
0115     unsigned int nlay = noiseMaps_[d]->ddd()->layers(true);
0116     std::pair<double, double> ranZ = noiseMaps_[d]->ddd()->rangeZ(true);
0117     std::pair<double, double> ranRAtZ = noiseMaps_[d]->ddd()->rangeR(ranZ.first, true);
0118     std::pair<double, double> ranR(ranRAtZ.first - plotMargin_, ranRAtZ.second + plotMargin_);
0119
0120     const std::vector<DetId> &detIdVec = noiseMaps_[d]->geom()->getValidDetIds();
0121     cout << "Subdetector:" << d << " has " << detIdVec.size() << " valid cells" << endl
0122          << "\t" << ranR.first << "<r<" << ranR.second << "\t" << ranZ.first << "<z<" << ranZ.second << endl;
0123
0124     //start histos
0125     TString baseName(Form("d%d_", d));
0126     TString title(d == DetId::HGCalEE ? "CEE" : "CEH_{Si}");
0127     Int_t nbinsR(100);
0128     for (unsigned int ilay = 0; ilay < nlay; ilay++) {
0129       //this layer histos
0130       int layer(ilay + 1);
0131       std::pair<DetId::Detector, int> key(d, layer);
0132       TString layerBaseName(Form("%slayer%d_", baseName.Data(), layer));
0133       TString layerTitle(Form("%s %d", title.Data(), layer));
0134       layerTitle += ";Radius [cm];";
0135       layerN_[key] = fs_->make<TH1F>(layerBaseName + "ncells", layerTitle + "Cells", nbinsR, ranR.first, ranR.second);
0136       layerCCE_[key] = fs_->make<TH1F>(layerBaseName + "cce", layerTitle + "<CCE>", nbinsR, ranR.first, ranR.second);
0137       layerNoise_[key] =
0138           fs_->make<TH1F>(layerBaseName + "noise", layerTitle + "<Noise> [fC]", nbinsR, ranR.first, ranR.second);
0139       layerIleak_[key] =
0140           fs_->make<TH1F>(layerBaseName + "ileak", layerTitle + "<I_{leak}> [#muA]", nbinsR, ranR.first, ranR.second);
0141       layerSN_[key] = fs_->make<TH1F>(layerBaseName + "sn", layerTitle + "<S/N>", nbinsR, ranR.first, ranR.second);
0142       layerF_[key] = fs_->make<TH1F>(
0143           layerBaseName + "fluence", layerTitle + "<F> [n_{eq}/cm^{2}]", nbinsR, ranR.first, ranR.second);
0144       layerGain_[key] = fs_->make<TH1F>(layerBaseName + "gain", layerTitle + "<Gain>", nbinsR, ranR.first, ranR.second);
0145       layerMipPeak_[key] =
0146           fs_->make<TH1F>(layerBaseName + "mippeak", layerTitle + "<MIP peak> [ADC]", nbinsR, ranR.first, ranR.second);
0147     }
0148
0149     //cce vs fluence
0150     for (unsigned int wafertype = 0; wafertype < (noiseMaps_[d]->getMipEqfC()).size(); ++wafertype) {
0151       std::pair<DetId::Detector, int> key2(d, wafertype);
0152       detCCEVsFluence_[key2] = fs_->make<TProfile>(
0153           baseName + Form("wt%d_", wafertype) + "cceVsFluence", title + ";<F> [n_{eq}/cm^{2}];<CCE>", 1000, 1e14, 1e16);
0154     }
0155
0156     //sub-detector histos
0157     title += ";Layer;Radius [cm];";
0158     detN_[d] =
0159         fs_->make<TH2F>(baseName + "ncells", title + "Cells", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0160     detCCE_[d] = fs_->make<TH2F>(baseName + "cce", title + "<CCE>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0161     detNoise_[d] =
0162         fs_->make<TH2F>(baseName + "noise", title + "<Noise> [fC]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0163     detIleak_[d] = fs_->make<TH2F>(
0164         baseName + "ileak", title + "<I_{leak}> [#muA]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0165     detSN_[d] = fs_->make<TH2F>(baseName + "sn", title + "<S/N>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0166     detF_[d] = fs_->make<TH2F>(
0167         baseName + "fluence", title + "<F> [n_{eq}/cm^{2}]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0168     detGain_[d] =
0169         fs_->make<TH2F>(baseName + "gain", title + "<Gain>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0170     detMipPeak_[d] = fs_->make<TH2F>(
0171         baseName + "mippeak", title + "<MIP peak> [ADC]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0172
0173     //fill histos
0174     for (const auto &cellId : detIdVec) {
0175       HGCSiliconDetId id(cellId.rawId());
0176       int layer = std::abs(id.layer());
0177       GlobalPoint pt = noiseMaps_[d]->geom()->getPosition(id);
0178       double r(pt.perp());
0179
0180       HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t gainToSet(HGCalSiNoiseMap<HGCSiliconDetId>::AUTO);
0181       if (ignoreGainSettings_)
0182         gainToSet = HGCalSiNoiseMap<HGCSiliconDetId>::q80fC;
0183       HGCalSiNoiseMap<HGCSiliconDetId>::SiCellOpCharacteristics siop =
0184           noiseMaps_[d]->getSiCellOpCharacteristics(id, gainToSet, aimMIPtoADC_);
0185
0186       //fill histos (layer,radius)
0187       detN_[d]->Fill(layer, r, 1);
0188       detCCE_[d]->Fill(layer, r, siop.core.cce);
0189       detNoise_[d]->Fill(layer, r, siop.core.noise);
0190       detSN_[d]->Fill(layer, r, siop.mipfC / siop.core.noise);
0191       detIleak_[d]->Fill(layer, r, siop.ileak);
0192       detF_[d]->Fill(layer, r, siop.fluence);
0193       detGain_[d]->Fill(layer, r, siop.core.gain + 1);
0194       detMipPeak_[d]->Fill(layer, r, siop.mipADC);
0195
0196       //per layer histograms
0197       std::pair<DetId::Detector, int> key(d, layer);
0198       layerN_[key]->Fill(r, 1);
0199       layerCCE_[key]->Fill(r, siop.core.cce);
0200       layerNoise_[key]->Fill(r, siop.core.noise);
0201       layerSN_[key]->Fill(r, siop.mipfC / siop.core.noise);
0202       layerIleak_[key]->Fill(r, siop.ileak);
0203       layerF_[key]->Fill(r, siop.fluence);
0204       layerGain_[key]->Fill(r, siop.core.gain + 1);
0205       layerMipPeak_[key]->Fill(r, siop.mipADC);
0206
0207       std::pair<DetId::Detector, int> key2(d, id.type());
0208       detCCEVsFluence_[key2]->Fill(siop.fluence, siop.core.cce);
0209     }
0210
0211     //normalize histos per cell counts
0212     detF_[d]->Divide(detN_[d]);
0213     detCCE_[d]->Divide(detN_[d]);
0214     detNoise_[d]->Divide(detN_[d]);
0215     detSN_[d]->Divide(detN_[d]);
0216     detIleak_[d]->Divide(detN_[d]);
0217     detGain_[d]->Divide(detN_[d]);
0218     detMipPeak_[d]->Divide(detN_[d]);
0219     for (unsigned int ilay = 0; ilay < nlay; ilay++) {
0220       int layer(ilay + 1);
0221       std::pair<DetId::Detector, int> key(d, layer);
0222       layerF_[key]->Divide(layerN_[key]);
0223       layerCCE_[key]->Divide(layerN_[key]);
0224       layerNoise_[key]->Divide(layerN_[key]);
0225       layerSN_[key]->Divide(layerN_[key]);
0226       layerIleak_[key]->Divide(layerN_[key]);
0227       layerGain_[key]->Divide(layerN_[key]);
0228       layerMipPeak_[key]->Divide(layerN_[key]);
0229     }
0230   }
0231 }
0232
0233 //define this as a plug-in
0234 DEFINE_FWK_MODULE(HGCSiNoiseMapAnalyzer);