HGCalSimAlgos/test/HFNoseNoiseMapAnalyzer.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 HFNoseNoiseMapAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0043 public:
0044   explicit HFNoseNoiseMapAnalyzer(const edm::ParameterSet &);
0045   ~HFNoseNoiseMapAnalyzer() 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<HFNoseDetId>>> 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 HFNoseNoiseMapAnalyzer::HFNoseNoiseMapAnalyzer(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   std::string doseMapURL("SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb_fluka_HFNose_3.7.20.12_Eta2.4.txt");
0076   unsigned int doseMapAlgo(iConfig.getParameter<unsigned int>("doseMapAlgo"));
0077   double scaleByDoseFactor(iConfig.getParameter<double>("scaleByDoseFactor"));
0078   std::vector<double> ileakParam(
0079       iConfig.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
0080   std::vector<double> cceParamFine(
0081       iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"));
0082   std::vector<double> cceParamThin(
0083       iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"));
0084   std::vector<double> cceParamThick(
0085       iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));
0086
0087   noiseMaps_[DetId::Forward] = std::unique_ptr<HGCalSiNoiseMap<HFNoseDetId>>(new HGCalSiNoiseMap<HFNoseDetId>);
0088   noiseMaps_[DetId::Forward]->setDoseMap(doseMapURL, doseMapAlgo);
0089   noiseMaps_[DetId::Forward]->setFluenceScaleFactor(scaleByDoseFactor);
0090   noiseMaps_[DetId::Forward]->setIleakParam(ileakParam);
0091   noiseMaps_[DetId::Forward]->setCceParam(cceParamFine, cceParamThin, cceParamThick);
0092
0093   aimMIPtoADC_ = iConfig.getParameter<int>("aimMIPtoADC");
0094   ignoreGainSettings_ = iConfig.getParameter<bool>("ignoreGainSettings");
0095 }
0096
0097 //
0098 HFNoseNoiseMapAnalyzer::~HFNoseNoiseMapAnalyzer() {}
0099
0100 //
0101 void HFNoseNoiseMapAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &es) {
0102   //get geometry
0103   const auto &geom = es.getHandle(tokGeom_);
0104
0105   std::vector<DetId::Detector> dets = {DetId::Forward};
0106   for (const auto &d : dets) {
0107     noiseMaps_[d]->setGeometry(geom->getSubdetectorGeometry(d, ForwardSubdetector::HFNose));
0108     //sub-detector boundaries
0109     unsigned int nlay = noiseMaps_[d]->ddd()->layers(true);
0110     std::pair<double, double> ranZ = noiseMaps_[d]->ddd()->rangeZ(true);
0111     std::pair<double, double> ranRAtZ = noiseMaps_[d]->ddd()->rangeR(ranZ.first, true);
0112     std::pair<double, double> ranR(ranRAtZ.first - plotMargin_, ranRAtZ.second + plotMargin_);
0113
0114     const std::vector<DetId> &detIdVec = noiseMaps_[d]->geom()->getValidDetIds();
0115     /*
0116     cout << "Subdetector:" << d << " has " << detIdVec.size() << " valid cells" << endl
0117          << "\t" << ranR.first << "<r<" << ranR.second << "\t" << ranZ.first << "<z<" << ranZ.second << endl;
0118     */
0119
0120     //start histos
0121     TString baseName(Form("d%d_", d));
0122     TString title("HFNose");
0123     Int_t nbinsR(100);
0124     for (unsigned int ilay = 0; ilay < nlay; ilay++) {
0125       //this layer histos
0126       int layer(ilay + 1);
0127       std::pair<DetId::Detector, int> key(d, layer);
0128       TString layerBaseName(Form("%slayer%d_", baseName.Data(), layer));
0129       TString layerTitle(Form("%s %d", title.Data(), layer));
0130       layerTitle += ";Radius [cm];";
0131       layerN_[key] = fs_->make<TH1F>(layerBaseName + "ncells", layerTitle + "Cells", nbinsR, ranR.first, ranR.second);
0132       layerCCE_[key] = fs_->make<TH1F>(layerBaseName + "cce", layerTitle + "<CCE>", nbinsR, ranR.first, ranR.second);
0133       layerNoise_[key] =
0134           fs_->make<TH1F>(layerBaseName + "noise", layerTitle + "<Noise> [fC]", nbinsR, ranR.first, ranR.second);
0135       layerIleak_[key] =
0136           fs_->make<TH1F>(layerBaseName + "ileak", layerTitle + "<I_{leak}> [#muA]", nbinsR, ranR.first, ranR.second);
0137       layerSN_[key] = fs_->make<TH1F>(layerBaseName + "sn", layerTitle + "<S/N>", nbinsR, ranR.first, ranR.second);
0138       layerF_[key] = fs_->make<TH1F>(
0139           layerBaseName + "fluence", layerTitle + "<F> [n_{eq}/cm^{2}]", nbinsR, ranR.first, ranR.second);
0140       layerGain_[key] = fs_->make<TH1F>(layerBaseName + "gain", layerTitle + "<Gain>", nbinsR, ranR.first, ranR.second);
0141       layerMipPeak_[key] =
0142           fs_->make<TH1F>(layerBaseName + "mippeak", layerTitle + "<MIP peak> [ADC]", nbinsR, ranR.first, ranR.second);
0143     }
0144
0145     //cce vs fluence
0146     for (unsigned int wafertype = 0; wafertype < (noiseMaps_[d]->getMipEqfC()).size(); ++wafertype) {
0147       std::pair<DetId::Detector, int> key2(d, wafertype);
0148       detCCEVsFluence_[key2] = fs_->make<TProfile>(
0149           baseName + Form("wt%d_", wafertype) + "cceVsFluence", title + ";<F> [n_{eq}/cm^{2}];<CCE>", 1000, 1e14, 1e16);
0150     }
0151
0152     //sub-detector histos
0153     title += ";Layer;Radius [cm];";
0154     detN_[d] =
0155         fs_->make<TH2F>(baseName + "ncells", title + "Cells", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0156     detCCE_[d] = fs_->make<TH2F>(baseName + "cce", title + "<CCE>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0157     detNoise_[d] =
0158         fs_->make<TH2F>(baseName + "noise", title + "<Noise> [fC]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0159     detIleak_[d] = fs_->make<TH2F>(
0160         baseName + "ileak", title + "<I_{leak}> [#muA]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0161     detSN_[d] = fs_->make<TH2F>(baseName + "sn", title + "<S/N>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0162     detF_[d] = fs_->make<TH2F>(
0163         baseName + "fluence", title + "<F> [n_{eq}/cm^{2}]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0164     detGain_[d] =
0165         fs_->make<TH2F>(baseName + "gain", title + "<Gain>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0166     detMipPeak_[d] = fs_->make<TH2F>(
0167         baseName + "mippeak", title + "<MIP peak> [ADC]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
0168
0169     //fill histos
0170     for (const auto &cellId : detIdVec) {
0171       HFNoseDetId id(cellId.rawId());
0172       int layer = std::abs(id.layer());
0173       GlobalPoint pt = noiseMaps_[d]->geom()->getPosition(id);
0174       double r(pt.perp());
0175
0176       HGCalSiNoiseMap<HFNoseDetId>::GainRange_t gainToSet(HGCalSiNoiseMap<HFNoseDetId>::AUTO);
0177       if (ignoreGainSettings_)
0178         gainToSet = HGCalSiNoiseMap<HFNoseDetId>::q80fC;
0179       HGCalSiNoiseMap<HFNoseDetId>::SiCellOpCharacteristics siop =
0180           noiseMaps_[d]->getSiCellOpCharacteristics(id, gainToSet, aimMIPtoADC_);
0181
0182       //fill histos (layer,radius)
0183       detN_[d]->Fill(layer, r, 1);
0184       detCCE_[d]->Fill(layer, r, siop.core.cce);
0185       detNoise_[d]->Fill(layer, r, siop.core.noise);
0186       detSN_[d]->Fill(layer, r, siop.mipfC / siop.core.noise);
0187       detIleak_[d]->Fill(layer, r, siop.ileak);
0188       detF_[d]->Fill(layer, r, siop.fluence);
0189       detGain_[d]->Fill(layer, r, siop.core.gain + 1);
0190       detMipPeak_[d]->Fill(layer, r, siop.mipADC);
0191
0192       //per layer histograms
0193       std::pair<DetId::Detector, int> key(d, layer);
0194       layerN_[key]->Fill(r, 1);
0195       layerCCE_[key]->Fill(r, siop.core.cce);
0196       layerNoise_[key]->Fill(r, siop.core.noise);
0197       layerSN_[key]->Fill(r, siop.mipfC / siop.core.noise);
0198       layerIleak_[key]->Fill(r, siop.ileak);
0199       layerF_[key]->Fill(r, siop.fluence);
0200       layerGain_[key]->Fill(r, siop.core.gain + 1);
0201       layerMipPeak_[key]->Fill(r, siop.mipADC);
0202
0203       std::pair<DetId::Detector, int> key2(d, id.type());
0204       detCCEVsFluence_[key2]->Fill(siop.fluence, siop.core.cce);
0205     }
0206
0207     //normalize histos per cell counts
0208     detF_[d]->Divide(detN_[d]);
0209     detCCE_[d]->Divide(detN_[d]);
0210     detNoise_[d]->Divide(detN_[d]);
0211     detSN_[d]->Divide(detN_[d]);
0212     detIleak_[d]->Divide(detN_[d]);
0213     detGain_[d]->Divide(detN_[d]);
0214     detMipPeak_[d]->Divide(detN_[d]);
0215     for (unsigned int ilay = 0; ilay < nlay; ilay++) {
0216       int layer(ilay + 1);
0217       std::pair<DetId::Detector, int> key(d, layer);
0218       layerF_[key]->Divide(layerN_[key]);
0219       layerCCE_[key]->Divide(layerN_[key]);
0220       layerNoise_[key]->Divide(layerN_[key]);
0221       layerSN_[key]->Divide(layerN_[key]);
0222       layerIleak_[key]->Divide(layerN_[key]);
0223       layerGain_[key]->Divide(layerN_[key]);
0224       layerMipPeak_[key]->Divide(layerN_[key]);
0225     }
0226   }
0227 }
0228
0229 //define this as a plug-in
0230 DEFINE_FWK_MODULE(HFNoseNoiseMapAnalyzer);