HGCalSimAlgos/test/HGCHEbackSignalScalerAnalyzer.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/HGCalSciNoiseMap.h"
0019
0020 //ROOT headers
0021 #include <TH2F.h>
0022
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/PluginManager/interface/ModuleDef.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027
0028 //STL headers
0029 #include <vector>
0030 #include <sstream>
0031 #include <string>
0032 #include <iomanip>
0033 #include <set>
0034
0035 #include "vdt/vdtMath.h"
0036
0037 //
0038 // class declaration
0039 //
0040
0041 class HGCHEbackSignalScalerAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0042 public:
0043   explicit HGCHEbackSignalScalerAnalyzer(const edm::ParameterSet&);
0044   ~HGCHEbackSignalScalerAnalyzer() override;
0045
0046 private:
0047   void beginJob() override {}
0048   void analyze(const edm::Event&, const edm::EventSetup&) override;
0049   void endJob() override {}
0050
0051   void createBinning(const std::vector<DetId>&);
0052   void printBoundaries();
0053
0054   // ----------member data ---------------------------
0055   edm::Service<TFileService> fs;
0056
0057   std::string doseMap_;
0058   uint32_t doseMapAlgo_;
0059   std::string sipmMap_;
0060   double refIdark_;
0061
0062   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0063   const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> dddToken_;
0064
0065   std::set<int> allLayers_, allIphi_, allIeta_;
0066
0067   const HGCalGeometry* gHGCal_;
0068   const HGCalDDDConstants* hgcCons_;
0069
0070   std::map<int, std::map<TString, TH1F*> > layerHistos_;
0071   std::map<TString, TH2F*> histos_;
0072 };
0073
0074 //
0075 // constructors and destructor
0076 //
0077 HGCHEbackSignalScalerAnalyzer::HGCHEbackSignalScalerAnalyzer(const edm::ParameterSet& iConfig)
0078     : doseMap_(iConfig.getParameter<std::string>("doseMap")),
0079       doseMapAlgo_(iConfig.getParameter<uint32_t>("doseMapAlgo")),
0080       sipmMap_(iConfig.getParameter<std::string>("sipmMap")),
0081       refIdark_(iConfig.getParameter<double>("referenceIdark")),
0082       geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag("", "HGCalHEScintillatorSensitive"))),
0083       dddToken_(
0084           esConsumes<HGCalDDDConstants, IdealGeometryRecord>(edm::ESInputTag("", "HGCalHEScintillatorSensitive"))) {
0085   usesResource("TFileService");
0086   fs->file().cd();
0087 }
0088
0089 HGCHEbackSignalScalerAnalyzer::~HGCHEbackSignalScalerAnalyzer() {}
0090
0091 //
0092 // member functions
0093 //
0094
0095 // ------------ method called on each new Event  ------------
0096 void HGCHEbackSignalScalerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0097   //get geometry
0098   const auto& geomhandle = iSetup.getHandle(geomToken_);
0099   if (!geomhandle.isValid()) {
0100     edm::LogError("HGCHEbackSignalScalerAnalyzer") << "Cannot get valid HGCalGeometry Object";
0101     return;
0102   }
0103   gHGCal_ = geomhandle.product();
0104   const std::vector<DetId>& detIdVec = gHGCal_->getValidDetIds();
0105   LogDebug("HGCHEbackSignalScalerAnalyzer") << "Total number of DetIDs: " << detIdVec.size();
0106
0107   //get ddd constants
0108   const auto& dddhandle = iSetup.getHandle(dddToken_);
0109   if (!dddhandle.isValid()) {
0110     edm::LogError("HGCHEbackSignalScalerAnalyzer") << "Cannot initiate HGCalDDDConstants";
0111     return;
0112   }
0113   hgcCons_ = dddhandle.product();
0114
0115   //setup maps
0116   createBinning(detIdVec);
0117   int minLayer(*std::min_element(allLayers_.begin(), allLayers_.end()));
0118   int maxLayer(*std::max_element(allLayers_.begin(), allLayers_.end()));
0119   int nLayers(maxLayer - minLayer + 1);
0120   int minIeta(*std::min_element(allIeta_.begin(), allIeta_.end()));
0121   int maxIeta(*std::max_element(allIeta_.begin(), allIeta_.end()));
0122   int nIeta(maxIeta - minIeta + 1);
0123   int minIphi(*std::min_element(allIphi_.begin(), allIphi_.end()));
0124   int maxIphi(*std::max_element(allIphi_.begin(), allIphi_.end()));
0125   int nIphi(maxIphi - minIphi + 1);
0126
0127   TString hnames[] = {"count", "radius", "dose", "fluence", "s", "n", "sn", "lysf", "gain", "thr"};
0128   TString htitles[] = {"tiles",
0129                        "#rho [cm]",
0130                        "<Dose> [krad]",
0131                        "<Fluence> [1 MeV n_{eq}/cm^{2}]",
0132                        "<S> [pe]",
0133                        "<#sigma_{N}> [pe]",
0134                        "<S/#sigma_{N}>",
0135                        "LY SF",
0136                        "<Gain>",
0137                        "<Threshold> [ADC]"};
0138   for (size_t i = 0; i < sizeof(hnames) / sizeof(TString); i++) {
0139     histos_[hnames[i] + "_ieta"] = fs->make<TH2F>(
0140         hnames[i] + "_ieta", ";Layer;i-#eta;" + htitles[i], nLayers, minLayer, maxLayer + 1, nIeta, minIeta, maxIeta);
0141     histos_[hnames[i] + "_iphi"] = fs->make<TH2F>(
0142         hnames[i] + "_iphi", ";Layer;i-#phi;" + htitles[i], nLayers, minLayer, maxLayer + 1, nIphi, minIphi, maxIphi);
0143
0144     //add per-layer profiles
0145     for (int ilay = minLayer; ilay <= maxLayer; ilay++) {
0146       if (layerHistos_.count(ilay) == 0) {
0147         std::map<TString, TH1F*> templ_map;
0148         layerHistos_[ilay] = templ_map;
0149       }
0150       TString hnameLay(Form("%s_lay%d_", hnames[i].Data(), ilay));
0151       layerHistos_[ilay][hnames[i] + "_ieta"] =
0152           fs->make<TH1F>(hnameLay + "ieta", ";i-#eta;" + htitles[i], nIeta, minIeta, maxIeta);
0153       layerHistos_[ilay][hnames[i] + "_iphi"] =
0154           fs->make<TH1F>(hnameLay + "iphi", ";i-#phi;" + htitles[i], nIphi, minIphi, maxIphi);
0155     }
0156   }
0157
0158   //instantiate scaler
0159   HGCalSciNoiseMap scal;
0160   scal.setDoseMap(doseMap_, doseMapAlgo_);
0161   scal.setReferenceDarkCurrent(refIdark_);
0162   scal.setSipmMap(sipmMap_);
0163   scal.setGeometry(gHGCal_);
0164
0165   //loop over valid detId from the HGCHEback
0166   for (auto myId : detIdVec) {
0167     HGCScintillatorDetId scId(myId.rawId());
0168
0169     int layer = std::abs(scId.layer());
0170     std::pair<int, int> ietaphi = scId.ietaphi();
0171     int ieta = std::abs(ietaphi.first);
0172     int iphi = std::abs(ietaphi.second);
0173
0174     //characteristics of this tile
0175     GlobalPoint global = gHGCal_->getPosition(scId);
0176     double radius = std::sqrt(std::pow(global.x(), 2) + std::pow(global.y(), 2));
0177     double dose = scal.getDoseValue(DetId::HGCalHSc, layer, radius);
0178     double fluence = scal.getFluenceValue(DetId::HGCalHSc, layer, radius);
0179     auto opChar = scal.scaleByDose(scId, radius);
0180     float s = opChar.s;
0181     float n = opChar.n;
0182     float sn = n > 0 ? s / n : -1.;
0183     float lysf = opChar.lySF;
0184     int gain = opChar.gain;
0185     int thrADC = opChar.thrADC;
0186
0187     histos_["count_ieta"]->Fill(layer, ieta, 1);
0188     histos_["count_iphi"]->Fill(layer, iphi, 1);
0189     histos_["radius_ieta"]->Fill(layer, ieta, radius);
0190     histos_["radius_iphi"]->Fill(layer, iphi, radius);
0191     histos_["dose_ieta"]->Fill(layer, ieta, dose);
0192     histos_["dose_iphi"]->Fill(layer, iphi, dose);
0193     histos_["fluence_ieta"]->Fill(layer, ieta, fluence);
0194     histos_["fluence_iphi"]->Fill(layer, iphi, fluence);
0195     histos_["s_ieta"]->Fill(layer, ieta, s);
0196     histos_["s_iphi"]->Fill(layer, iphi, s);
0197     histos_["n_ieta"]->Fill(layer, ieta, n);
0198     histos_["n_iphi"]->Fill(layer, iphi, n);
0199     histos_["sn_ieta"]->Fill(layer, ieta, sn);
0200     histos_["sn_iphi"]->Fill(layer, iphi, sn);
0201     histos_["lysf_ieta"]->Fill(layer, ieta, lysf);
0202     histos_["lysf_iphi"]->Fill(layer, iphi, lysf);
0203     histos_["gain_ieta"]->Fill(layer, ieta, gain);
0204     histos_["gain_iphi"]->Fill(layer, iphi, gain);
0205     histos_["thr_ieta"]->Fill(layer, ieta, thrADC);
0206     histos_["thr_iphi"]->Fill(layer, iphi, thrADC);
0207
0208     if (layerHistos_.count(layer) == 0) {
0209       continue;
0210     }
0211
0212     //per layer
0213     layerHistos_[layer]["count_ieta"]->Fill(ieta, 1);
0214     layerHistos_[layer]["count_iphi"]->Fill(iphi, 1);
0215     layerHistos_[layer]["radius_ieta"]->Fill(ieta, radius);
0216     layerHistos_[layer]["radius_iphi"]->Fill(iphi, radius);
0217     layerHistos_[layer]["dose_ieta"]->Fill(ieta, dose);
0218     layerHistos_[layer]["dose_iphi"]->Fill(iphi, dose);
0219     layerHistos_[layer]["fluence_ieta"]->Fill(ieta, fluence);
0220     layerHistos_[layer]["fluence_iphi"]->Fill(iphi, fluence);
0221     layerHistos_[layer]["s_ieta"]->Fill(ieta, s);
0222     layerHistos_[layer]["s_iphi"]->Fill(iphi, s);
0223     layerHistos_[layer]["n_ieta"]->Fill(ieta, n);
0224     layerHistos_[layer]["n_iphi"]->Fill(iphi, n);
0225     layerHistos_[layer]["sn_ieta"]->Fill(ieta, sn);
0226     layerHistos_[layer]["sn_iphi"]->Fill(iphi, sn);
0227     layerHistos_[layer]["lysf_ieta"]->Fill(ieta, lysf);
0228     layerHistos_[layer]["lysf_iphi"]->Fill(iphi, lysf);
0229     layerHistos_[layer]["gain_ieta"]->Fill(ieta, gain);
0230     layerHistos_[layer]["gain_iphi"]->Fill(iphi, gain);
0231     layerHistos_[layer]["thr_ieta"]->Fill(ieta, thrADC);
0232     layerHistos_[layer]["thr_iphi"]->Fill(iphi, thrADC);
0233   }
0234
0235   //normalize by the number of tiles to have the mean stored
0236   for (auto h : histos_) {
0237     TString name(h.first);
0238     if (name.Contains("count_"))
0239       continue;
0240     TString pfix(name.Contains("ieta") ? "ieta" : "iphi");
0241     h.second->Divide(histos_["count_" + pfix]);
0242   }
0243
0244   for (auto h : layerHistos_) {
0245     int lay(h.first);
0246     for (auto hh : layerHistos_[lay]) {
0247       TString name(hh.first);
0248       if (name.Contains("count_"))
0249         continue;
0250       TString pfix(name.Contains("ieta") ? "ieta" : "iphi");
0251       hh.second->Divide(layerHistos_[lay]["count_" + pfix]);
0252     }
0253   }
0254 }
0255
0256 //
0257 void HGCHEbackSignalScalerAnalyzer::createBinning(const std::vector<DetId>& detIdVec) {
0258   //loop over detids to find possible values
0259   for (auto myId : detIdVec) {
0260     HGCScintillatorDetId scId(myId.rawId());
0261     int layer = std::abs(scId.layer());
0262     std::pair<int, int> ietaphi = scId.ietaphi();
0263
0264     allLayers_.insert(std::abs(layer));
0265     allIeta_.insert(std::abs(ietaphi.first));
0266     allIphi_.insert(std::abs(ietaphi.second));
0267   }
0268 }
0269
0270 //define this as a plug-in
0271 DEFINE_FWK_MODULE(HGCHEbackSignalScalerAnalyzer);