Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-25 03:21:54 

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalLaserAnalyzerYousi
0004 // Class:      EcalLaserAnalyzerYousi
0005 //
0006 /**\class EcalLaserAnalyzerYousi EcalLaserAnalyzerYousi.cc CalibCalorimetry/EcalLaserAnalyzerYousi/src/EcalLaserAnalyzerYousi.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Yousi Ma
0015 //         Created:  Tue Jun 19 23:06:36 CEST 2007
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <iostream>
0022 #include <fstream>
0023 
0024 // user include files
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0027 
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 
0034 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0035 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0036 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0037 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0038 #include "DataFormats/EcalRawData/interface/EcalRawDataCollections.h"
0039 
0040 #include "TFile.h"
0041 #include "TH1F.h"
0042 #include "TH2F.h"
0043 #include "TROOT.h"
0044 #include "TF1.h"
0045 #include "TNtuple.h"
0046 #include "TDirectory.h"
0047 
0048 //
0049 // class decleration
0050 //
0051 
0052 class EcalLaserAnalyzerYousi : public edm::one::EDAnalyzer<> {
0053 public:
0054   explicit EcalLaserAnalyzerYousi(const edm::ParameterSet &);
0055   ~EcalLaserAnalyzerYousi() override = default;
0056 
0057 private:
0058   void beginJob() override;
0059   void analyze(const edm::Event &, const edm::EventSetup &) override;
0060   void endJob() override;
0061 
0062   // ----------member data ---------------------------
0063   // Declare histograms and ROOT trees, etc.
0064   TH1F *C_APD[1700];
0065   TH1F *C_APDPN[1700];
0066   TH1F *C_PN[1700];
0067   TH1F *C_J[1700];
0068   TH2F *C_APDPN_J[1700];
0069 
0070   TH1F *peakAPD[2];
0071   TH1F *peakAPDPN[2];
0072   TH1F *APD_LM[9];
0073   TH1F *APDPN_LM[9];
0074   TH2F *APDPN_J_LM[9];
0075   TH2F *APDPN_J_H[2];
0076 
0077   //fixme make declare and init separate
0078   TH2F *APD;
0079   TH2F *APD_RMS;
0080   TH2F *APDPN;
0081   TH2F *APDPN_RMS;
0082   TH2F *PN;
0083   TH2F *APDPN_J;
0084   TH2F *APDPN_C;
0085 
0086   TH1F *FitHist;
0087   TH2F *Count;
0088 
0089   TFile *fPN;
0090   TFile *fAPD;
0091   TFile *fROOT;
0092 
0093   TNtuple *C_Tree[1700];
0094 
0095   //parameters
0096 
0097   const std::string hitCollection_;
0098   const std::string hitProducer_;
0099   //  std::string PNFileName_ ;
0100   //  std::string ABFileName_ ;
0101   const std::string outFileName_;
0102   const std::string SM_;
0103   const std::string Run_;
0104   const std::string digiProducer_;
0105   const std::string PNdigiCollection_;
0106   const edm::EDGetTokenT<EcalRawDataCollection> rawDataToken_;
0107   const edm::EDGetTokenT<EBUncalibratedRecHitCollection> hitToken_;
0108   const edm::EDGetTokenT<EcalPnDiodeDigiCollection> pnDigiToken_;
0109 };
0110 
0111 //
0112 // constants, enums and typedefs
0113 //
0114 
0115 //
0116 // static data member definitions
0117 //
0118 
0119 //
0120 // constructors and destructor
0121 //
0122 EcalLaserAnalyzerYousi::EcalLaserAnalyzerYousi(const edm::ParameterSet &iConfig)
0123     : hitCollection_(iConfig.getUntrackedParameter<std::string>("hitCollection")),
0124       hitProducer_(iConfig.getUntrackedParameter<std::string>("hitProducer")),
0125       outFileName_(iConfig.getUntrackedParameter<std::string>("outFileName")),
0126       SM_(iConfig.getUntrackedParameter<std::string>("SM")),
0127       Run_(iConfig.getUntrackedParameter<std::string>("Run")),
0128       digiProducer_(iConfig.getUntrackedParameter<std::string>("digiProducer")),
0129       PNdigiCollection_(iConfig.getUntrackedParameter<std::string>("PNdigiCollection")),
0130       rawDataToken_(consumes<EcalRawDataCollection>(edm::InputTag(digiProducer_))),
0131       hitToken_(consumes<EBUncalibratedRecHitCollection>(edm::InputTag(hitProducer_, hitCollection_))),
0132       pnDigiToken_(consumes<EcalPnDiodeDigiCollection>(edm::InputTag(digiProducer_, PNdigiCollection_))) {
0133   //now do what ever initialization is needed
0134   //get the PN and AB file names
0135   //get the output file names, digi producers, etc
0136 
0137   //  PNFileName_ = iConfig.getUntrackedParameter<std::string>("PNFileName");
0138   //  ABFileName_ = iConfig.getUntrackedParameter<std::string>("ABFileName");
0139 }
0140 
0141 //
0142 // member functions
0143 //
0144 
0145 // ------------ method called to for each event  ------------
0146 void EcalLaserAnalyzerYousi::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0147   using namespace edm;
0148 
0149   //  if ( fPN->IsOpen() ) { edm::LogInfo("EcalLaserAnalyzerYousi") <<"fPN is open in analyze OKAAAAAAAAYYY \n\n"; }
0150 
0151   const edm::Handle<EcalRawDataCollection> &DCCHeaders = iEvent.getHandle(rawDataToken_);
0152 
0153   EcalDCCHeaderBlock::EcalDCCEventSettings settings = DCCHeaders->begin()->getEventSettings();
0154 
0155   int wavelength = settings.wavelength;
0156 
0157   //   std::cout<<"wavelength: "<<wavelength<<"\n\n";
0158 
0159   if (wavelength != 0) {
0160     return;
0161   }  //only process blue laser
0162 
0163   const edm::Handle<EBUncalibratedRecHitCollection> &hits = iEvent.getHandle(hitToken_);
0164 
0165   if (!hits.isValid()) {
0166     edm::LogError("EcalLaserAnalyzerYousi")
0167         << "Cannot get product:  EBRecHitCollection from: " << hitCollection_ << " - returning.\n\n";
0168     return;
0169   }
0170 
0171   const edm::Handle<EcalPnDiodeDigiCollection> &pndigis = iEvent.getHandle(pnDigiToken_);
0172   if (!pndigis.isValid()) {
0173     edm::LogError("EcalLaserAnalyzerYousi") << "Cannot get product:  EBdigiCollection from: getHandle - returning.\n\n";
0174     return;
0175   }
0176 
0177   Float_t PN_amp[5];
0178 
0179   //do some averaging over each pair of PNs
0180   for (int j = 0; j < 5; ++j) {
0181     PN_amp[j] = 0;
0182     for (int z = 0; z < 2; ++z) {
0183       FitHist->Reset();
0184       TF1 peakFit("peakFit", "[0] +[1]*x +[2]*x^2", 30, 50);
0185       TF1 pedFit("pedFit", "[0]", 0, 5);
0186 
0187       for (int k = 0; k < 50; k++) {
0188         FitHist->SetBinContent(k, (*pndigis)[j + z * 5].sample(k).adc());
0189       }
0190       pedFit.SetParameter(0, 750);
0191       FitHist->Fit(&pedFit, "RQI");
0192       Float_t ped = pedFit.GetParameter(0);
0193 
0194       Int_t maxbin = FitHist->GetMaximumBin();
0195       peakFit.SetRange(FitHist->GetBinCenter(maxbin) - 4 * FitHist->GetBinWidth(maxbin),
0196                        FitHist->GetBinCenter(maxbin) + 4 * FitHist->GetBinWidth(maxbin));
0197       peakFit.SetParameters(750, 4, -.05);
0198       FitHist->Fit(&peakFit, "RQI");
0199       Float_t max = peakFit.Eval(-peakFit.GetParameter(1) / (2 * peakFit.GetParameter(2)));
0200       if (ped != max) {
0201         PN_amp[j] = PN_amp[j] + max - ped;
0202       } else {
0203         PN_amp[j] = PN_amp[j] + max;
0204       }
0205 
0206     }  //end of z loop
0207     PN_amp[j] = PN_amp[j] / 2.0;
0208   }  //end of j loop
0209 
0210   //do some real PN, APD calculations
0211 
0212   //FIXME. previously used .info files to get time, what to do now?
0213 
0214   //   TNtuple *Time = new TNtuple("Time", "Time", "Time");
0215   //   Int_t iTime = Get_Time(Input_File);
0216   //   Time->Fill(iTime);
0217 
0218   Float_t fTree[7];
0219 
0220   //   b->GetEntry(EVT);
0221   EBDetId ID;
0222   Float_t theAPD;
0223   Float_t thePN;
0224   Float_t Jitter;
0225   Float_t Chi2;
0226   Int_t CN = hits->size();
0227   //            cout<<"num CN: "<<CN<<endl;
0228   for (int j = 0; j < CN; j++) {
0229     ID = (*hits)[j].id();
0230     theAPD = (*hits)[j].amplitude();
0231     Jitter = (*hits)[j].jitter();
0232     Chi2 = (*hits)[j].chi2();
0233     thePN = PN_amp[(ID.ic() + 299) / 400];
0234 
0235     //              cout<<"THE APD: "<<theAPD<<endl;
0236     //              cout<<"THE PN: "<<thePN<<endl;
0237 
0238     C_APD[ID.ic() - 1]->Fill(theAPD);
0239     C_APDPN[ID.ic() - 1]->Fill(theAPD / thePN);
0240     C_PN[ID.ic() - 1]->Fill(thePN);
0241     C_J[ID.ic() - 1]->Fill(Jitter);
0242     C_APDPN_J[ID.ic() - 1]->Fill(Jitter, theAPD / thePN);
0243     APDPN_J->Fill(Jitter, theAPD / thePN);
0244     APDPN_C->Fill(Chi2, theAPD / thePN);
0245     fTree[0] = theAPD;
0246     fTree[1] = thePN;
0247     fTree[2] = theAPD / thePN;
0248     fTree[3] = Jitter;
0249     fTree[4] = Chi2;
0250     fTree[5] = (*hits)[j].pedestal();
0251     fTree[6] = iEvent.id().event();
0252     C_Tree[ID.ic() - 1]->Fill(fTree);
0253     if (((ID.ic() - 1) % 20 > 9) || ((ID.ic() - 1) < 100)) {
0254       peakAPD[0]->Fill(theAPD);
0255       peakAPDPN[0]->Fill(theAPD / thePN);
0256       APDPN_J_H[0]->Fill(Jitter, theAPD / thePN);
0257     } else {
0258       peakAPD[1]->Fill(theAPD);
0259       peakAPDPN[1]->Fill(theAPD / thePN);
0260       APDPN_J_H[1]->Fill(Jitter, theAPD / thePN);
0261     }
0262     if ((ID.ic() - 1) < 100) {
0263       APD_LM[0]->Fill(theAPD);
0264       APDPN_LM[0]->Fill(theAPD / thePN);
0265       APDPN_J_LM[0]->Fill(Jitter, theAPD / thePN);
0266     } else {
0267       Int_t index;
0268       if (((ID.ic() - 1) % 20) < 10) {
0269         index = ((ID.ic() - 101) / 400) * 2 + 1;
0270         APD_LM[index]->Fill(theAPD);
0271         APDPN_LM[index]->Fill(theAPD / thePN);
0272         APDPN_J_LM[index]->Fill(Jitter, theAPD / thePN);
0273       } else {
0274         index = ((ID.ic() - 101) / 400) * 2 + 2;
0275         APD_LM[index]->Fill(theAPD);
0276         APDPN_LM[index]->Fill(theAPD / thePN);
0277         APDPN_J_LM[index]->Fill(Jitter, theAPD / thePN);
0278       }
0279     }
0280   }  //end of CN loop
0281 
0282   //now that you got the PN and APD's, make the ntuples. done
0283 
0284   //vec from ROOT version should correspond to hits_itr or something similar. done
0285 
0286   //check WL from PNdiodedigi, should be ==0, o.w (blue data only). don't process. done
0287 
0288   //get PN pulse, and do fitting of pulse. i.e. fill hist with PN.apd() or equivalent. done
0289 
0290   //fit to first 5 for PED, and 30-50 bins for pulse (poly2 for the moment). done
0291 }
0292 
0293 // ------------ method called once each job just before starting event loop  ------------
0294 void EcalLaserAnalyzerYousi::beginJob() {
0295   edm::LogInfo("EcalLaserAnalyzerYousi") << "running laser analyzer \n\n";
0296 
0297   fROOT = new TFile(outFileName_.c_str(), "RECREATE");
0298   fROOT->cd();
0299 
0300   //init all the histos and files?
0301   APD = new TH2F("APD", "APD", 85, 0., 85., 20, 0., 20.);
0302   APD_RMS = new TH2F("APD_RMS", "APD_RMS", 85, 0., 85., 20, 0., 20.);
0303   APDPN = new TH2F("APDPN", "APDPN", 85, 0., 85., 20, 0., 20.);
0304   APDPN_RMS = new TH2F("APDPN_RMS", "APDPN_RMS", 85, 0., 85., 20, 0., 20.);
0305   PN = new TH2F("PN", "PN", 85, 0., 85., 20, 0., 20.);
0306   APDPN_J = new TH2F("JittervAPDPN", "JittervAPDPN", 250, 3., 7., 250, 1., 2.);
0307   APDPN_C = new TH2F("Chi2vAPDPN", "Chi2vAPDPN", 250, 0., 50., 250, 0., 5.0);
0308   FitHist = new TH1F("FitHist", "FitHist", 50, 0, 50);
0309   Count = new TH2F("Count", "Count", 85, 0., 1., 20, 0., 1.);
0310 
0311   for (int i = 0; i < 1700; i++) {
0312     std::ostringstream name_1;
0313     std::ostringstream name_2;
0314     std::ostringstream name_3;
0315     std::ostringstream name_4;
0316     std::ostringstream name_5;
0317     name_1 << "C_APD_" << i + 1;
0318     name_2 << "C_APDPN_" << i + 1;
0319     name_3 << "C_PN_" << i + 1;
0320     name_4 << "C_J_" << i + 1;
0321     name_5 << "C_APDPN_J_" << i + 1;
0322     C_APD[i] = new TH1F(name_1.str().c_str(), name_1.str().c_str(), 2500, 0., 5000.);
0323     C_APDPN[i] = new TH1F(name_2.str().c_str(), name_2.str().c_str(), 20000, 0., 25.);
0324     C_PN[i] = new TH1F(name_3.str().c_str(), name_3.str().c_str(), 1000, 0., 4000.);
0325     C_J[i] = new TH1F(name_4.str().c_str(), name_4.str().c_str(), 250, 3.0, 7.);
0326     C_APDPN_J[i] = new TH2F(name_5.str().c_str(), name_5.str().c_str(), 250, 3.0, 6., 250, 1., 2.2);
0327   }
0328 
0329   for (int i = 0; i < 2; i++) {
0330     std::ostringstream aname_1;
0331     std::ostringstream aname_2;
0332     std::ostringstream aname_3;
0333     aname_1 << "peakAPD_" << i;
0334     aname_2 << "peakAPDPN_" << i;
0335     aname_3 << "JittervAPDPN_Half_" << i;
0336     peakAPD[i] = new TH1F(aname_1.str().c_str(), aname_1.str().c_str(), 1000, 0., 5000.);
0337     peakAPDPN[i] = new TH1F(aname_2.str().c_str(), aname_2.str().c_str(), 1000, 0., 8.);
0338     APDPN_J_H[i] = new TH2F(aname_3.str().c_str(), aname_3.str().c_str(), 250, 3., 7., 250, 1., 2.2);
0339   }
0340 
0341   for (int i = 0; i < 9; i++) {
0342     std::ostringstream bname_1;
0343     std::ostringstream bname_2;
0344     std::ostringstream bname_3;
0345     bname_1 << "APD_LM_" << i;
0346     bname_2 << "APDPN_LM_" << i;
0347     bname_3 << "APDPN_J_LM_" << i;
0348     APD_LM[i] = new TH1F(bname_1.str().c_str(), bname_1.str().c_str(), 500, 0., 5000.);
0349     APDPN_LM[i] = new TH1F(bname_2.str().c_str(), bname_2.str().c_str(), 500, 0., 8.);
0350     APDPN_J_LM[i] = new TH2F(bname_3.str().c_str(), bname_3.str().c_str(), 250, 3., 7., 250, 1., 2.2);
0351   }
0352 
0353   //get the PN file. or don't get, and read from event.
0354 
0355   //don't need to get AB, it will be read in via framework poolsource = ???
0356 
0357   //configure the final NTuple
0358   std::ostringstream varlist;
0359   varlist << "APD:PN:APDPN:Jitter:Chi2:ped:EVT";
0360   for (int i = 0; i < 1700; i++) {
0361     std::ostringstream name;
0362     name << "C_Tree_" << i + 1;
0363     C_Tree[i] = (TNtuple *)new TNtuple(name.str().c_str(), name.str().c_str(), varlist.str().c_str());
0364   }
0365 }
0366 
0367 // ------------ method called once each job just after ending the event loop  ------------
0368 void EcalLaserAnalyzerYousi::endJob() {
0369   //write the file (get ouput file name first).
0370   TFile *fROOT = (TFile *)new TFile(outFileName_.c_str(), "RECREATE");
0371 
0372   //  TDirectory *DIR = fROOT->Get(Run_.c_str());
0373   TDirectory *DIR;
0374   //  if(DIR == NULL){
0375   DIR = fROOT->mkdir(Run_.c_str());
0376   //  }
0377   DIR->cd();
0378   for (int j = 0; j < 1700; j++) {
0379     Float_t min_r, max_r;
0380     Float_t RMS, Sigma, K;
0381     Int_t iCount;
0382     TF1 *gs3;
0383 
0384     RMS = C_APD[j]->GetRMS();
0385     APD_RMS->SetBinContent(85 - (j / 20), 20 - (j % 20), RMS);
0386     Sigma = 999999;
0387     K = 2.5;
0388     iCount = 0;
0389     while (Sigma > RMS) {
0390       min_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) - K * RMS;
0391       max_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) + K * RMS;
0392       TF1 gs1("gs1", "gaus", min_r, max_r);
0393       C_APD[j]->Fit(&gs1, "RQI");
0394       Sigma = gs1.GetParameter(2);
0395       K = K * 1.5;
0396       iCount++;
0397       if (iCount > 2) {
0398         C_APD[j]->Fit("gaus", "QI");
0399         break;
0400       }
0401     }
0402 
0403     RMS = C_APDPN[j]->GetRMS();
0404     APDPN_RMS->SetBinContent(85 - (j / 20), 20 - (j % 20), RMS);
0405     Sigma = 999999;
0406     K = 2.5;
0407     iCount = 0;
0408     while (Sigma > RMS) {
0409       min_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) - K * RMS;
0410       max_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) + K * RMS;
0411       TF1 gs2("gs2", "gaus", min_r, max_r);
0412       C_APDPN[j]->Fit(&gs2, "RQI");
0413       Sigma = gs2.GetParameter(2);
0414       K = K * 1.5;
0415       iCount++;
0416       if (iCount > 2) {
0417         C_APDPN[j]->Fit("gaus", "QI");
0418         break;
0419       }
0420     }
0421 
0422     TF1 *newgs1;
0423     TF1 *newgs2;
0424 
0425     C_PN[j]->Fit("gaus", "Q");
0426     C_APD[j]->Fit("gaus", "QI");
0427     C_APDPN[j]->Fit("gaus", "QI");
0428     C_APD[j]->Write("", TObject::kOverwrite);
0429     C_APDPN[j]->Write("", TObject::kOverwrite);
0430     C_PN[j]->Write("", TObject::kOverwrite);
0431     C_J[j]->Write("", TObject::kOverwrite);
0432     C_APDPN_J[j]->Write("", TObject::kOverwrite);
0433     newgs1 = C_APD[j]->GetFunction("gaus");
0434     newgs2 = C_APDPN[j]->GetFunction("gaus");
0435     gs3 = C_PN[j]->GetFunction("gaus");
0436     Float_t theAPD = newgs1->GetParameter(1);
0437     APD->SetBinContent(85 - (j / 20), 20 - (j % 20), theAPD);
0438     Float_t theAPDPN = newgs2->GetParameter(1);
0439     APDPN->SetBinContent(85 - (j / 20), 20 - (j % 20), theAPDPN);
0440     Float_t thePN = gs3->GetParameter(1);
0441     //      cout<<"LOOK HERE thePN = "<< thePN<<endl;
0442     PN->SetBinContent(85 - (j / 20), 20 - (j % 20), thePN);
0443     C_Tree[j]->Write("", TObject::kOverwrite);
0444   }
0445 
0446   for (int i = 0; i < 9; i++) {
0447     APD_LM[i]->Write("", TObject::kOverwrite);
0448     APDPN_LM[i]->Write("", TObject::kOverwrite);
0449     APDPN_J_LM[i]->Write("", TObject::kOverwrite);
0450   }
0451 
0452   //  Time->Write("", TObject::kOverwrite);
0453   APD->Write("", TObject::kOverwrite);
0454   APD_RMS->Write("", TObject::kOverwrite);
0455   APDPN_RMS->Write("", TObject::kOverwrite);
0456   APDPN->Write("", TObject::kOverwrite);
0457   APDPN_J->Write("", TObject::kOverwrite);
0458   APDPN_C->Write("", TObject::kOverwrite);
0459   PN->Write("", TObject::kOverwrite);
0460   peakAPD[0]->Write("", TObject::kOverwrite);
0461   peakAPD[1]->Write("", TObject::kOverwrite);
0462   peakAPDPN[0]->Write("", TObject::kOverwrite);
0463   peakAPDPN[1]->Write("", TObject::kOverwrite);
0464   APDPN_J_H[0]->Write("", TObject::kOverwrite);
0465   APDPN_J_H[1]->Write("", TObject::kOverwrite);
0466 
0467   // don't Make plots
0468   //  fROOT->Close();
0469 
0470   //   fPN->Close();
0471   //   fAPD->Close();
0472 
0473   fROOT->Write();
0474   //   fROOT->Close();
0475 }
0476 
0477 //define this as a plug-in
0478 DEFINE_FWK_MODULE(EcalLaserAnalyzerYousi);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!