HcalSimAlgos/test/SiPMNonlinearityAnalyzer.cc

0001 // system include files
0002 #include <memory>
0003
0004 // user include files
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 //ROOT headers
0015 #include <TProfile.h>
0016 #include <TH1.h>
0017 #include <TF1.h>
0018
0019 //SiPM headers
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 //STL headers
0030 #include <vector>
0031 #include <utility>
0032 #include <sstream>
0033 #include <string>
0034 #include <algorithm>
0035 #include <cmath>
0036
0037 //
0038 // class declaration
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   // ----------member data ---------------------------
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 // constructors and destructor
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 // member functions
0088 //
0089
0090 // ------------ method called on each new Event  ------------
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   //instantiate SiPM classes
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   //shift to TS 4
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       //smear pes according to Y11 time distribution
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       //reset SiPM
0126       sipm.setNCells(pixels);
0127
0128       //evaluate SiPM response
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       //fill profiles
0142       profio->Fill(sumPE, sumHits);
0143       profoi->Fill(sumHits, sumPE);
0144     }
0145     npe += npeStep;
0146   }
0147
0148   //calculate correction factor
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   //fit with requested function
0155   corr->Fit(fitname.c_str(), "Q");
0156
0157   //print fit results
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 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0167 void SiPMNonlinearityAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0168   //The following says we do not know what parameters are allowed so do no validation
0169   // Please change this to state exactly what you do use, even if it is no parameters
0170   edm::ParameterSetDescription desc;
0171   desc.setUnknown();
0172   descriptions.addDefault(desc);
0173 }
0174 //define this as a plug-in
0175 DEFINE_FWK_MODULE(SiPMNonlinearityAnalyzer);