File indexing completed on 2024-04-06 12:29:32
0001
0002 #include <memory>
0003
0004
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/StreamID.h"
0013
0014
0015 #include <TProfile.h>
0016 #include <TH1.h>
0017 #include <TF1.h>
0018
0019
0020 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPM.h"
0021 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.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 "FWCore/Utilities/interface/RandomNumberGenerator.h"
0027 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0028
0029
0030 #include <vector>
0031 #include <utility>
0032 #include <sstream>
0033 #include <string>
0034 #include <algorithm>
0035 #include <cmath>
0036
0037
0038
0039
0040
0041 class SiPMNonlinearityAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0042 public:
0043 explicit SiPMNonlinearityAnalyzer(const edm::ParameterSet&);
0044 ~SiPMNonlinearityAnalyzer() override;
0045
0046 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048 private:
0049 void beginJob() override {}
0050 void analyze(const edm::Event&, const edm::EventSetup&) override;
0051 void endJob() override {}
0052
0053
0054 edm::Service<TFileService> fs;
0055 unsigned pixels;
0056 unsigned npeMin, npeMax, npeStep, nReps;
0057 unsigned nBins, binMin, binMax;
0058 double tau, dt;
0059 unsigned nPreciseBins;
0060 std::string fitname;
0061 unsigned signalShape;
0062 };
0063
0064
0065
0066
0067 SiPMNonlinearityAnalyzer::SiPMNonlinearityAnalyzer(const edm::ParameterSet& iConfig)
0068 : pixels(iConfig.getParameter<unsigned>("pixels")),
0069 npeMin(iConfig.getParameter<unsigned>("npeMin")),
0070 npeMax(iConfig.getParameter<unsigned>("npeMax")),
0071 npeStep(iConfig.getParameter<unsigned>("npeStep")),
0072 nReps(iConfig.getParameter<unsigned>("nReps")),
0073 nBins(iConfig.getParameter<unsigned>("nBins")),
0074 binMin(iConfig.getParameter<unsigned>("binMin")),
0075 binMax(iConfig.getParameter<unsigned>("binMax")),
0076 tau(iConfig.getParameter<double>("tau")),
0077 dt(iConfig.getParameter<double>("dt")),
0078 nPreciseBins(iConfig.getParameter<unsigned>("nPreciseBins")),
0079 fitname(iConfig.getParameter<std::string>("fitname")),
0080 signalShape(iConfig.getParameter<unsigned>("signalShape")) {
0081 usesResource("TFileService");
0082 }
0083
0084 SiPMNonlinearityAnalyzer::~SiPMNonlinearityAnalyzer() {}
0085
0086
0087
0088
0089
0090
0091 void SiPMNonlinearityAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0092 edm::Service<edm::RandomNumberGenerator> rng;
0093 CLHEP::HepRandomEngine* engine = &rng->getEngine(iEvent.streamID());
0094
0095
0096 HcalSiPM sipm(pixels, tau);
0097
0098 fs->file().cd();
0099 TProfile* profio = fs->make<TProfile>("input_vs_output", "", nBins, binMin, binMax);
0100 profio->GetXaxis()->SetTitle("input [pe]");
0101 profio->GetYaxis()->SetTitle("output [pe]");
0102 TProfile* profoi = fs->make<TProfile>("output_vs_input", "", nBins, binMin, binMax);
0103 profoi->GetXaxis()->SetTitle("output [pe]");
0104 profoi->GetYaxis()->SetTitle("input [pe]");
0105 TH1F* corr = fs->make<TH1F>("correction", "", nBins, binMin, binMax);
0106 corr->GetXaxis()->SetTitle("output [pe]");
0107 corr->GetYaxis()->SetTitle("correction factor");
0108
0109
0110 static const int tsOffset = 100;
0111 unsigned npe = npeMin;
0112 while (npe <= npeMax) {
0113 if (npe % 1000 == 0)
0114 std::cout << "npe = " << npe << std::endl;
0115 for (unsigned rep = 0; rep < nReps; ++rep) {
0116
0117 std::vector<unsigned> photonHist(nPreciseBins, 0);
0118 for (unsigned pe = 0; pe < npe; ++pe) {
0119 double t_pe = HcalPulseShapes::generatePhotonTime(engine, signalShape);
0120 int t_bin = int(t_pe + tsOffset + 0.5);
0121 if (t_bin > 0 and (static_cast<unsigned>(t_bin) < photonHist.size()))
0122 photonHist[t_bin] += 1;
0123 }
0124
0125
0126 sipm.setNCells(pixels);
0127
0128
0129 double elapsedTime = 0.;
0130 unsigned sumPE = 0;
0131 double sumHits = 0.;
0132 for (unsigned tbin = 0; tbin < photonHist.size(); ++tbin) {
0133 unsigned pe = photonHist[tbin];
0134 if (pe > 0) {
0135 sumPE += pe;
0136 sumHits += sipm.hitCells(engine, pe, 0., elapsedTime);
0137 }
0138 elapsedTime += dt;
0139 }
0140
0141
0142 profio->Fill(sumPE, sumHits);
0143 profoi->Fill(sumHits, sumPE);
0144 }
0145 npe += npeStep;
0146 }
0147
0148
0149 for (int b = 0; b < profoi->GetNbinsX(); ++b) {
0150 if (profoi->GetBinContent(b) == 0)
0151 continue;
0152 corr->SetBinContent(b, profoi->GetBinContent(b) / profoi->GetBinCenter(b));
0153 }
0154
0155 corr->Fit(fitname.c_str(), "Q");
0156
0157
0158 TF1* fit = corr->GetFunction(fitname.c_str());
0159 std::cout << "chisquare / ndf = " << fit->GetChisquare() << " / " << fit->GetNDF() << std::endl;
0160 std::cout << "parameters = ";
0161 std::copy(
0162 fit->GetParameters(), fit->GetParameters() + fit->GetNpar(), std::ostream_iterator<double>(std::cout, ", "));
0163 std::cout << std::endl;
0164 }
0165
0166
0167 void SiPMNonlinearityAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0168
0169
0170 edm::ParameterSetDescription desc;
0171 desc.setUnknown();
0172 descriptions.addDefault(desc);
0173 }
0174
0175 DEFINE_FWK_MODULE(SiPMNonlinearityAnalyzer);