Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-08 08:15:44 

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;
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   std::string hitCollection_;
0098   std::string hitProducer_;
0099   //  std::string PNFileName_ ;
0100   //  std::string ABFileName_ ;
0101   std::string outFileName_;
0102   std::string SM_;
0103   std::string Run_;
0104   std::string digiProducer_;
0105   std::string PNdigiCollection_;
0106 };
0107 
0108 //
0109 // constants, enums and typedefs
0110 //
0111 
0112 //
0113 // static data member definitions
0114 //
0115 
0116 //
0117 // constructors and destructor
0118 //
0119 EcalLaserAnalyzerYousi::EcalLaserAnalyzerYousi(const edm::ParameterSet &iConfig) {
0120   //now do what ever initialization is needed
0121   //get the PN and AB file names
0122   //get the output file names, digi producers, etc
0123 
0124   hitCollection_ = iConfig.getUntrackedParameter<std::string>("hitCollection");
0125   hitProducer_ = iConfig.getUntrackedParameter<std::string>("hitProducer");
0126   //  PNFileName_ = iConfig.getUntrackedParameter<std::string>("PNFileName");
0127   //  ABFileName_ = iConfig.getUntrackedParameter<std::string>("ABFileName");
0128   outFileName_ = iConfig.getUntrackedParameter<std::string>("outFileName");
0129   SM_ = iConfig.getUntrackedParameter<std::string>("SM");
0130   Run_ = iConfig.getUntrackedParameter<std::string>("Run");
0131   digiProducer_ = iConfig.getUntrackedParameter<std::string>("digiProducer");
0132   PNdigiCollection_ = iConfig.getUntrackedParameter<std::string>("PNdigiCollection");
0133 }
0134 
0135 EcalLaserAnalyzerYousi::~EcalLaserAnalyzerYousi() {
0136   // do anything here that needs to be done at desctruction time
0137   // (e.g. close files, deallocate resources etc.)
0138 }
0139 
0140 //
0141 // member functions
0142 //
0143 
0144 // ------------ method called to for each event  ------------
0145 void EcalLaserAnalyzerYousi::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0146   using namespace edm;
0147 
0148   //  if ( fPN->IsOpen() ) { edm::LogInfo("EcalLaserAnalyzerYousi") <<"fPN is open in analyze OKAAAAAAAAYYY \n\n"; }
0149 
0150   edm::Handle<EcalRawDataCollection> DCCHeaders;
0151   iEvent.getByLabel(digiProducer_, DCCHeaders);
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   edm::Handle<EBUncalibratedRecHitCollection> hits;
0164 
0165   try {
0166     iEvent.getByLabel(hitProducer_, hitCollection_, hits);
0167     //     iEvent.getByType(hits);
0168   } catch (std::exception &ex) {
0169     LogError("EcalLaserAnalyzerYousi") << "Cannot get product:  EBRecHitCollection from: " << hitCollection_
0170                                        << " - returning.\n\n";
0171     //    return;
0172   }
0173 
0174   edm::Handle<EcalPnDiodeDigiCollection> pndigis;
0175 
0176   try {
0177     //     iEvent.getByLabel(hitProducer_, hits);
0178     iEvent.getByLabel(digiProducer_, PNdigiCollection_, pndigis);
0179     //iEvent.getByType( pndigis );
0180   } catch (std::exception &ex) {
0181     LogError("EcalLaserAnalyzerYousi") << "Cannot get product:  EBdigiCollection from: "
0182                                        << "getbytype"
0183                                        << " - returning.\n\n";
0184     //    return;
0185   }
0186 
0187   Float_t PN_amp[5];
0188 
0189   //do some averaging over each pair of PNs
0190   for (int j = 0; j < 5; ++j) {
0191     PN_amp[j] = 0;
0192     for (int z = 0; z < 2; ++z) {
0193       FitHist->Reset();
0194       TF1 peakFit("peakFit", "[0] +[1]*x +[2]*x^2", 30, 50);
0195       TF1 pedFit("pedFit", "[0]", 0, 5);
0196 
0197       for (int k = 0; k < 50; k++) {
0198         FitHist->SetBinContent(k, (*pndigis)[j + z * 5].sample(k).adc());
0199       }
0200       pedFit.SetParameter(0, 750);
0201       FitHist->Fit(&pedFit, "RQI");
0202       Float_t ped = pedFit.GetParameter(0);
0203 
0204       Int_t maxbin = FitHist->GetMaximumBin();
0205       peakFit.SetRange(FitHist->GetBinCenter(maxbin) - 4 * FitHist->GetBinWidth(maxbin),
0206                        FitHist->GetBinCenter(maxbin) + 4 * FitHist->GetBinWidth(maxbin));
0207       peakFit.SetParameters(750, 4, -.05);
0208       FitHist->Fit(&peakFit, "RQI");
0209       Float_t max = peakFit.Eval(-peakFit.GetParameter(1) / (2 * peakFit.GetParameter(2)));
0210       if (ped != max) {
0211         PN_amp[j] = PN_amp[j] + max - ped;
0212       } else {
0213         PN_amp[j] = PN_amp[j] + max;
0214       }
0215 
0216     }  //end of z loop
0217     PN_amp[j] = PN_amp[j] / 2.0;
0218   }  //end of j loop
0219 
0220   //do some real PN, APD calculations
0221 
0222   //FIXME. previously used .info files to get time, what to do now?
0223 
0224   //   TNtuple *Time = new TNtuple("Time", "Time", "Time");
0225   //   Int_t iTime = Get_Time(Input_File);
0226   //   Time->Fill(iTime);
0227 
0228   Float_t fTree[7];
0229 
0230   //   b->GetEntry(EVT);
0231   EBDetId ID;
0232   Float_t theAPD;
0233   Float_t thePN;
0234   Float_t Jitter;
0235   Float_t Chi2;
0236   Int_t CN = hits->size();
0237   //            cout<<"num CN: "<<CN<<endl;
0238   for (int j = 0; j < CN; j++) {
0239     ID = (*hits)[j].id();
0240     theAPD = (*hits)[j].amplitude();
0241     Jitter = (*hits)[j].jitter();
0242     Chi2 = (*hits)[j].chi2();
0243     thePN = PN_amp[(ID.ic() + 299) / 400];
0244 
0245     //              cout<<"THE APD: "<<theAPD<<endl;
0246     //              cout<<"THE PN: "<<thePN<<endl;
0247 
0248     C_APD[ID.ic() - 1]->Fill(theAPD);
0249     C_APDPN[ID.ic() - 1]->Fill(theAPD / thePN);
0250     C_PN[ID.ic() - 1]->Fill(thePN);
0251     C_J[ID.ic() - 1]->Fill(Jitter);
0252     C_APDPN_J[ID.ic() - 1]->Fill(Jitter, theAPD / thePN);
0253     APDPN_J->Fill(Jitter, theAPD / thePN);
0254     APDPN_C->Fill(Chi2, theAPD / thePN);
0255     fTree[0] = theAPD;
0256     fTree[1] = thePN;
0257     fTree[2] = theAPD / thePN;
0258     fTree[3] = Jitter;
0259     fTree[4] = Chi2;
0260     fTree[5] = (*hits)[j].pedestal();
0261     fTree[6] = iEvent.id().event();
0262     C_Tree[ID.ic() - 1]->Fill(fTree);
0263     if (((ID.ic() - 1) % 20 > 9) || ((ID.ic() - 1) < 100)) {
0264       peakAPD[0]->Fill(theAPD);
0265       peakAPDPN[0]->Fill(theAPD / thePN);
0266       APDPN_J_H[0]->Fill(Jitter, theAPD / thePN);
0267     } else {
0268       peakAPD[1]->Fill(theAPD);
0269       peakAPDPN[1]->Fill(theAPD / thePN);
0270       APDPN_J_H[1]->Fill(Jitter, theAPD / thePN);
0271     }
0272     if ((ID.ic() - 1) < 100) {
0273       APD_LM[0]->Fill(theAPD);
0274       APDPN_LM[0]->Fill(theAPD / thePN);
0275       APDPN_J_LM[0]->Fill(Jitter, theAPD / thePN);
0276     } else {
0277       Int_t index;
0278       if (((ID.ic() - 1) % 20) < 10) {
0279         index = ((ID.ic() - 101) / 400) * 2 + 1;
0280         APD_LM[index]->Fill(theAPD);
0281         APDPN_LM[index]->Fill(theAPD / thePN);
0282         APDPN_J_LM[index]->Fill(Jitter, theAPD / thePN);
0283       } else {
0284         index = ((ID.ic() - 101) / 400) * 2 + 2;
0285         APD_LM[index]->Fill(theAPD);
0286         APDPN_LM[index]->Fill(theAPD / thePN);
0287         APDPN_J_LM[index]->Fill(Jitter, theAPD / thePN);
0288       }
0289     }
0290   }  //end of CN loop
0291 
0292   //now that you got the PN and APD's, make the ntuples. done
0293 
0294   //vec from ROOT version should correspond to hits_itr or something similar. done
0295 
0296   //check WL from PNdiodedigi, should be ==0, o.w (blue data only). don't process. done
0297 
0298   //get PN pulse, and do fitting of pulse. i.e. fill hist with PN.apd() or equivalent. done
0299 
0300   //fit to first 5 for PED, and 30-50 bins for pulse (poly2 for the moment). done
0301 }
0302 
0303 // ------------ method called once each job just before starting event loop  ------------
0304 void EcalLaserAnalyzerYousi::beginJob() {
0305   edm::LogInfo("EcalLaserAnalyzerYousi") << "running laser analyzer \n\n";
0306 
0307   fROOT = new TFile(outFileName_.c_str(), "RECREATE");
0308   fROOT->cd();
0309 
0310   //init all the histos and files?
0311   APD = new TH2F("APD", "APD", 85, 0., 85., 20, 0., 20.);
0312   APD_RMS = new TH2F("APD_RMS", "APD_RMS", 85, 0., 85., 20, 0., 20.);
0313   APDPN = new TH2F("APDPN", "APDPN", 85, 0., 85., 20, 0., 20.);
0314   APDPN_RMS = new TH2F("APDPN_RMS", "APDPN_RMS", 85, 0., 85., 20, 0., 20.);
0315   PN = new TH2F("PN", "PN", 85, 0., 85., 20, 0., 20.);
0316   APDPN_J = new TH2F("JittervAPDPN", "JittervAPDPN", 250, 3., 7., 250, 1., 2.);
0317   APDPN_C = new TH2F("Chi2vAPDPN", "Chi2vAPDPN", 250, 0., 50., 250, 0., 5.0);
0318   FitHist = new TH1F("FitHist", "FitHist", 50, 0, 50);
0319   Count = new TH2F("Count", "Count", 85, 0., 1., 20, 0., 1.);
0320 
0321   for (int i = 0; i < 1700; i++) {
0322     std::ostringstream name_1;
0323     std::ostringstream name_2;
0324     std::ostringstream name_3;
0325     std::ostringstream name_4;
0326     std::ostringstream name_5;
0327     name_1 << "C_APD_" << i + 1;
0328     name_2 << "C_APDPN_" << i + 1;
0329     name_3 << "C_PN_" << i + 1;
0330     name_4 << "C_J_" << i + 1;
0331     name_5 << "C_APDPN_J_" << i + 1;
0332     C_APD[i] = new TH1F(name_1.str().c_str(), name_1.str().c_str(), 2500, 0., 5000.);
0333     C_APDPN[i] = new TH1F(name_2.str().c_str(), name_2.str().c_str(), 20000, 0., 25.);
0334     C_PN[i] = new TH1F(name_3.str().c_str(), name_3.str().c_str(), 1000, 0., 4000.);
0335     C_J[i] = new TH1F(name_4.str().c_str(), name_4.str().c_str(), 250, 3.0, 7.);
0336     C_APDPN_J[i] = new TH2F(name_5.str().c_str(), name_5.str().c_str(), 250, 3.0, 6., 250, 1., 2.2);
0337   }
0338 
0339   for (int i = 0; i < 2; i++) {
0340     std::ostringstream aname_1;
0341     std::ostringstream aname_2;
0342     std::ostringstream aname_3;
0343     aname_1 << "peakAPD_" << i;
0344     aname_2 << "peakAPDPN_" << i;
0345     aname_3 << "JittervAPDPN_Half_" << i;
0346     peakAPD[i] = new TH1F(aname_1.str().c_str(), aname_1.str().c_str(), 1000, 0., 5000.);
0347     peakAPDPN[i] = new TH1F(aname_2.str().c_str(), aname_2.str().c_str(), 1000, 0., 8.);
0348     APDPN_J_H[i] = new TH2F(aname_3.str().c_str(), aname_3.str().c_str(), 250, 3., 7., 250, 1., 2.2);
0349   }
0350 
0351   for (int i = 0; i < 9; i++) {
0352     std::ostringstream bname_1;
0353     std::ostringstream bname_2;
0354     std::ostringstream bname_3;
0355     bname_1 << "APD_LM_" << i;
0356     bname_2 << "APDPN_LM_" << i;
0357     bname_3 << "APDPN_J_LM_" << i;
0358     APD_LM[i] = new TH1F(bname_1.str().c_str(), bname_1.str().c_str(), 500, 0., 5000.);
0359     APDPN_LM[i] = new TH1F(bname_2.str().c_str(), bname_2.str().c_str(), 500, 0., 8.);
0360     APDPN_J_LM[i] = new TH2F(bname_3.str().c_str(), bname_3.str().c_str(), 250, 3., 7., 250, 1., 2.2);
0361   }
0362 
0363   //get the PN file. or don't get, and read from event.
0364 
0365   //don't need to get AB, it will be read in via framework poolsource = ???
0366 
0367   //configure the final NTuple
0368   std::ostringstream varlist;
0369   varlist << "APD:PN:APDPN:Jitter:Chi2:ped:EVT";
0370   for (int i = 0; i < 1700; i++) {
0371     std::ostringstream name;
0372     name << "C_Tree_" << i + 1;
0373     C_Tree[i] = (TNtuple *)new TNtuple(name.str().c_str(), name.str().c_str(), varlist.str().c_str());
0374   }
0375 }
0376 
0377 // ------------ method called once each job just after ending the event loop  ------------
0378 void EcalLaserAnalyzerYousi::endJob() {
0379   //write the file (get ouput file name first).
0380   TFile *fROOT = (TFile *)new TFile(outFileName_.c_str(), "RECREATE");
0381 
0382   //  TDirectory *DIR = fROOT->Get(Run_.c_str());
0383   TDirectory *DIR;
0384   //  if(DIR == NULL){
0385   DIR = fROOT->mkdir(Run_.c_str());
0386   //  }
0387   DIR->cd();
0388   for (int j = 0; j < 1700; j++) {
0389     Float_t min_r, max_r;
0390     Float_t RMS, Sigma, K;
0391     Int_t iCount;
0392     TF1 *gs3;
0393 
0394     RMS = C_APD[j]->GetRMS();
0395     APD_RMS->SetBinContent(85 - (j / 20), 20 - (j % 20), RMS);
0396     Sigma = 999999;
0397     K = 2.5;
0398     iCount = 0;
0399     while (Sigma > RMS) {
0400       min_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) - K * RMS;
0401       max_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) + K * RMS;
0402       TF1 gs1("gs1", "gaus", min_r, max_r);
0403       C_APD[j]->Fit(&gs1, "RQI");
0404       Sigma = gs1.GetParameter(2);
0405       K = K * 1.5;
0406       iCount++;
0407       if (iCount > 2) {
0408         C_APD[j]->Fit("gaus", "QI");
0409         break;
0410       }
0411     }
0412 
0413     RMS = C_APDPN[j]->GetRMS();
0414     APDPN_RMS->SetBinContent(85 - (j / 20), 20 - (j % 20), RMS);
0415     Sigma = 999999;
0416     K = 2.5;
0417     iCount = 0;
0418     while (Sigma > RMS) {
0419       min_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) - K * RMS;
0420       max_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) + K * RMS;
0421       TF1 gs2("gs2", "gaus", min_r, max_r);
0422       C_APDPN[j]->Fit(&gs2, "RQI");
0423       Sigma = gs2.GetParameter(2);
0424       K = K * 1.5;
0425       iCount++;
0426       if (iCount > 2) {
0427         C_APDPN[j]->Fit("gaus", "QI");
0428         break;
0429       }
0430     }
0431 
0432     TF1 *newgs1;
0433     TF1 *newgs2;
0434 
0435     C_PN[j]->Fit("gaus", "Q");
0436     C_APD[j]->Fit("gaus", "QI");
0437     C_APDPN[j]->Fit("gaus", "QI");
0438     C_APD[j]->Write("", TObject::kOverwrite);
0439     C_APDPN[j]->Write("", TObject::kOverwrite);
0440     C_PN[j]->Write("", TObject::kOverwrite);
0441     C_J[j]->Write("", TObject::kOverwrite);
0442     C_APDPN_J[j]->Write("", TObject::kOverwrite);
0443     newgs1 = C_APD[j]->GetFunction("gaus");
0444     newgs2 = C_APDPN[j]->GetFunction("gaus");
0445     gs3 = C_PN[j]->GetFunction("gaus");
0446     Float_t theAPD = newgs1->GetParameter(1);
0447     APD->SetBinContent(85 - (j / 20), 20 - (j % 20), theAPD);
0448     Float_t theAPDPN = newgs2->GetParameter(1);
0449     APDPN->SetBinContent(85 - (j / 20), 20 - (j % 20), theAPDPN);
0450     Float_t thePN = gs3->GetParameter(1);
0451     //      cout<<"LOOK HERE thePN = "<< thePN<<endl;
0452     PN->SetBinContent(85 - (j / 20), 20 - (j % 20), thePN);
0453     C_Tree[j]->Write("", TObject::kOverwrite);
0454   }
0455 
0456   for (int i = 0; i < 9; i++) {
0457     APD_LM[i]->Write("", TObject::kOverwrite);
0458     APDPN_LM[i]->Write("", TObject::kOverwrite);
0459     APDPN_J_LM[i]->Write("", TObject::kOverwrite);
0460   }
0461 
0462   //  Time->Write("", TObject::kOverwrite);
0463   APD->Write("", TObject::kOverwrite);
0464   APD_RMS->Write("", TObject::kOverwrite);
0465   APDPN_RMS->Write("", TObject::kOverwrite);
0466   APDPN->Write("", TObject::kOverwrite);
0467   APDPN_J->Write("", TObject::kOverwrite);
0468   APDPN_C->Write("", TObject::kOverwrite);
0469   PN->Write("", TObject::kOverwrite);
0470   peakAPD[0]->Write("", TObject::kOverwrite);
0471   peakAPD[1]->Write("", TObject::kOverwrite);
0472   peakAPDPN[0]->Write("", TObject::kOverwrite);
0473   peakAPDPN[1]->Write("", TObject::kOverwrite);
0474   APDPN_J_H[0]->Write("", TObject::kOverwrite);
0475   APDPN_J_H[1]->Write("", TObject::kOverwrite);
0476 
0477   // don't Make plots
0478   //  fROOT->Close();
0479 
0480   //   fPN->Close();
0481   //   fAPD->Close();
0482 
0483   fROOT->Write();
0484   //   fROOT->Close();
0485 }
0486 
0487 //define this as a plug-in
0488 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!