EcalSimpleTBAnalysis/src/EcalSimpleTBAnalyzer.cc

0001 /**\class EcalSimpleTBAnalyzer
0002
0003  Description: <one line class summary>
0004
0005  Implementation:
0006      <Notes on implementation>
0007 */
0008 //
0009 //
0010 //
0011
0012 #include "RecoTBCalo/EcalSimpleTBAnalysis/interface/EcalSimpleTBAnalyzer.h"
0013 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0014 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0015 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0016 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRecInfo.h"
0017 #include "TBDataFormats/EcalTBObjects/interface/EcalTBTDCRecInfo.h"
0018 #include "TBDataFormats/EcalTBObjects/interface/EcalTBEventHeader.h"
0019
0020 //#include<fstream>
0021
0022 #include "TFile.h"
0023 #include "TH1.h"
0024 #include "TH2.h"
0025 #include "TF1.h"
0026
0027 #include <iostream>
0028 #include <string>
0029 #include <stdexcept>
0030 //
0031 // constants, enums and typedefs
0032 //
0033
0034 //
0035 // static data member definitions
0036 //
0037
0038 //
0039 // constructors and destructor
0040 //
0041
0042 //========================================================================
0043 EcalSimpleTBAnalyzer::EcalSimpleTBAnalyzer(const edm::ParameterSet& iConfig)
0044     : xtalInBeam_(0)
0045 //========================================================================
0046 {
0047   //now do what ever initialization is needed
0048   rootfile_ = iConfig.getUntrackedParameter<std::string>("rootfile", "ecalSimpleTBanalysis.root");
0049   digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
0050   digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
0051   hitCollection_ = iConfig.getParameter<std::string>("hitCollection");
0052   hitProducer_ = iConfig.getParameter<std::string>("hitProducer");
0053   hodoRecInfoCollection_ = iConfig.getParameter<std::string>("hodoRecInfoCollection");
0054   hodoRecInfoProducer_ = iConfig.getParameter<std::string>("hodoRecInfoProducer");
0055   tdcRecInfoCollection_ = iConfig.getParameter<std::string>("tdcRecInfoCollection");
0056   tdcRecInfoProducer_ = iConfig.getParameter<std::string>("tdcRecInfoProducer");
0057   eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
0058   eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
0059
0060   std::cout << "EcalSimpleTBAnalyzer: fetching hitCollection: " << hitCollection_.c_str() << " produced by "
0061             << hitProducer_.c_str() << std::endl;
0062 }
0063
0064 //========================================================================
0065 EcalSimpleTBAnalyzer::~EcalSimpleTBAnalyzer()
0066 //========================================================================
0067 {
0068   // do anything here that needs to be done at desctruction time
0069   // (e.g. close files, deallocate resources etc.)
0070   // Amplitude vs TDC offset
0071   //   if (h_ampltdc)
0072   //   delete h_ampltdc;
0073
0074   //   // Reconstructed energies
0075   //   delete h_e1x1;
0076   //   delete h_e3x3;
0077   //   delete h_e5x5;
0078
0079   //   delete h_bprofx;
0080   //   delete h_bprofy;
0081
0082   //   delete h_qualx;
0083   //   delete h_qualy;
0084
0085   //   delete h_slopex;
0086   //   delete h_slopey;
0087
0088   //   delete h_mapx;
0089   //   delete h_mapy;
0090 }
0091
0092 //========================================================================
0093 void EcalSimpleTBAnalyzer::beginJob() {
0094   //========================================================================
0095
0096   // Amplitude vs TDC offset
0097   h_ampltdc = new TH2F("h_ampltdc", "Max Amplitude vs TDC offset", 100, 0., 1., 1000, 0., 4000.);
0098
0099   // Reconstructed energies
0100   h_tableIsMoving = new TH1F("h_tableIsMoving", "TableIsMoving", 100000, 0., 100000.);
0101
0102   h_e1x1 = new TH1F("h_e1x1", "E1x1 energy", 1000, 0., 4000.);
0103   h_e3x3 = new TH1F("h_e3x3", "E3x3 energy", 1000, 0., 4000.);
0104   h_e5x5 = new TH1F("h_e5x5", "E5x5 energy", 1000, 0., 4000.);
0105
0106   h_e1x1_center = new TH1F("h_e1x1_center", "E1x1 energy", 1000, 0., 4000.);
0107   h_e3x3_center = new TH1F("h_e3x3_center", "E3x3 energy", 1000, 0., 4000.);
0108   h_e5x5_center = new TH1F("h_e5x5_center", "E5x5 energy", 1000, 0., 4000.);
0109
0110   h_e1e9 = new TH1F("h_e1e9", "E1/E9 ratio", 600, 0., 1.2);
0111   h_e1e25 = new TH1F("h_e1e25", "E1/E25 ratio", 600, 0., 1.2);
0112   h_e9e25 = new TH1F("h_e9e25", "E9/E25 ratio", 600, 0., 1.2);
0113
0114   h_bprofx = new TH1F("h_bprofx", "Beam Profile X", 100, -20., 20.);
0115   h_bprofy = new TH1F("h_bprofy", "Beam Profile Y", 100, -20., 20.);
0116
0117   h_qualx = new TH1F("h_qualx", "Beam Quality X", 5000, 0., 5.);
0118   h_qualy = new TH1F("h_qualy", "Beam Quality X", 5000, 0., 5.);
0119
0120   h_slopex = new TH1F("h_slopex", "Beam Slope X", 500, -5e-4, 5e-4);
0121   h_slopey = new TH1F("h_slopey", "Beam Slope Y", 500, -5e-4, 5e-4);
0122
0123   char hname[50];
0124   char htitle[50];
0125   for (unsigned int icry = 0; icry < 25; icry++) {
0126     sprintf(hname, "h_mapx_%d", icry);
0127     sprintf(htitle, "Max Amplitude vs X %d", icry);
0128     h_mapx[icry] = new TH2F(hname, htitle, 80, -20, 20, 1000, 0., 4000.);
0129     sprintf(hname, "h_mapy_%d", icry);
0130     sprintf(htitle, "Max Amplitude vs Y %d", icry);
0131     h_mapy[icry] = new TH2F(hname, htitle, 80, -20, 20, 1000, 0., 4000.);
0132   }
0133
0134   h_e1e9_mapx = new TH2F("h_e1e9_mapx", "E1/E9 vs X", 80, -20, 20, 600, 0., 1.2);
0135   h_e1e9_mapy = new TH2F("h_e1e9_mapy", "E1/E9 vs Y", 80, -20, 20, 600, 0., 1.2);
0136
0137   h_e1e25_mapx = new TH2F("h_e1e25_mapx", "E1/E25 vs X", 80, -20, 20, 600, 0., 1.2);
0138   h_e1e25_mapy = new TH2F("h_e1e25_mapy", "E1/E25 vs Y", 80, -20, 20, 600, 0., 1.2);
0139
0140   h_e9e25_mapx = new TH2F("h_e9e25_mapx", "E9/E25 vs X", 80, -20, 20, 600, 0., 1.2);
0141   h_e9e25_mapy = new TH2F("h_e9e25_mapy", "E9/E25 vs Y", 80, -20, 20, 600, 0., 1.2);
0142
0143   h_Shape_ = new TH2F("h_Shape_", "Xtal in Beam Shape", 250, 0, 10, 350, 0, 3500);
0144 }
0145
0146 //========================================================================
0147 void EcalSimpleTBAnalyzer::endJob() {
0148   //========================================================================
0149
0150   TFile f(rootfile_.c_str(), "RECREATE");
0151
0152   // Amplitude vs TDC offset
0153   h_ampltdc->Write();
0154
0155   // Reconstructed energies
0156   h_e1x1->Write();
0157   h_e3x3->Write();
0158   h_e5x5->Write();
0159
0160   h_e1x1_center->Write();
0161   h_e3x3_center->Write();
0162   h_e5x5_center->Write();
0163
0164   h_e1e9->Write();
0165   h_e1e25->Write();
0166   h_e9e25->Write();
0167
0168   h_bprofx->Write();
0169   h_bprofy->Write();
0170
0171   h_qualx->Write();
0172   h_qualy->Write();
0173
0174   h_slopex->Write();
0175   h_slopey->Write();
0176
0177   h_Shape_->Write();
0178
0179   for (unsigned int icry = 0; icry < 25; icry++) {
0180     h_mapx[icry]->Write();
0181     h_mapy[icry]->Write();
0182   }
0183
0184   h_e1e9_mapx->Write();
0185   h_e1e9_mapy->Write();
0186
0187   h_e1e25_mapx->Write();
0188   h_e1e25_mapy->Write();
0189
0190   h_e9e25_mapx->Write();
0191   h_e9e25_mapy->Write();
0192
0193   h_tableIsMoving->Write();
0194
0195   f.Close();
0196 }
0197
0198 //
0199 // member functions
0200 //
0201
0202 //========================================================================
0203 void EcalSimpleTBAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0204   //========================================================================
0205
0206   using namespace edm;
0207   using namespace cms;
0208
0209   Handle<EBDigiCollection> pdigis;
0210   const EBDigiCollection* digis = nullptr;
0211   //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
0212   iEvent.getByLabel(digiProducer_, digiCollection_, pdigis);
0213   if (pdigis.isValid()) {
0214     digis = pdigis.product();  // get a ptr to the product
0215     //iEvent.getByLabel( hitProducer_, phits);
0216   } else {
0217     edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << digiCollection_;
0218   }
0219
0220   // fetch the digis and compute signal amplitude
0221   Handle<EBUncalibratedRecHitCollection> phits;
0222   const EBUncalibratedRecHitCollection* hits = nullptr;
0223   //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
0224   iEvent.getByLabel(hitProducer_, hitCollection_, phits);
0225   if (phits.isValid()) {
0226     hits = phits.product();  // get a ptr to the product
0227     //iEvent.getByLabel( hitProducer_, phits);
0228   } else {
0229     edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hitCollection_;
0230   }
0231
0232   Handle<EcalTBHodoscopeRecInfo> pHodo;
0233   const EcalTBHodoscopeRecInfo* recHodo = nullptr;
0234   //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
0235   iEvent.getByLabel(hodoRecInfoProducer_, hodoRecInfoCollection_, pHodo);
0236   if (pHodo.isValid()) {
0237     recHodo = pHodo.product();  // get a ptr to the product
0238   } else {
0239     edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hodoRecInfoCollection_;
0240   }
0241
0242   Handle<EcalTBTDCRecInfo> pTDC;
0243   const EcalTBTDCRecInfo* recTDC = nullptr;
0244   //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
0245   iEvent.getByLabel(tdcRecInfoProducer_, tdcRecInfoCollection_, pTDC);
0246   if (pTDC.isValid()) {
0247     recTDC = pTDC.product();  // get a ptr to the product
0248   } else {
0249     edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << tdcRecInfoCollection_;
0250   }
0251
0252   Handle<EcalTBEventHeader> pEventHeader;
0253   const EcalTBEventHeader* evtHeader = nullptr;
0254   //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
0255   iEvent.getByLabel(eventHeaderProducer_, pEventHeader);
0256   if (pEventHeader.isValid()) {
0257     evtHeader = pEventHeader.product();  // get a ptr to the product
0258   } else {
0259     edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << eventHeaderProducer_;
0260   }
0261
0262   if (!hits)
0263     return;
0264
0265   if (!recTDC)
0266     return;
0267
0268   if (!recHodo)
0269     return;
0270
0271   if (!evtHeader)
0272     return;
0273
0274   if (hits->empty())
0275     return;
0276
0277   if (evtHeader->tableIsMoving())
0278     h_tableIsMoving->Fill(evtHeader->eventNumber());
0279
0280   // Crystal hit by beam
0281   if (xtalInBeam_.null()) {
0282     xtalInBeam_ = EBDetId(1, evtHeader->crystalInBeam(), EBDetId::SMCRYSTALMODE);
0283     std::cout << "Xtal In Beam is " << xtalInBeam_.ic() << std::endl;
0284   } else if (xtalInBeam_ != EBDetId(1, evtHeader->crystalInBeam(), EBDetId::SMCRYSTALMODE))
0285     return;
0286
0287   if (evtHeader->tableIsMoving())
0288     return;
0289
0290   //    EBDetId maxHitId(0);
0291   //    float maxHit= -999999.;
0292
0293   //    for(EBUncalibratedRecHitCollection::const_iterator ithit = hits->begin(); ithit != hits->end(); ++ithit)
0294   //      {
0295   //        if (ithit->amplitude()>=maxHit)
0296   //     {
0297   //       maxHit=ithit->amplitude();
0298   //       maxHitId=ithit->id();
0299   //     }
0300
0301   //      }
0302
0303   //    if (maxHitId==EBDetId(0))
0304   //      return;
0305
0306   //   EBDetId maxHitId(1,704,EBDetId::SMCRYSTALMODE);
0307
0308   //Find EBDetId in a 5x5 Matrix (to be substituted by the Selector code)
0309   // Something like
0310   // EBFixedWindowSelector<EcalUncalibratedRecHit> Simple5x5Matrix(hits,maxHitId,5,5);
0311   // std::vector<EcalUncalibratedRecHit> Energies5x5 = Simple5x5Matrix.getHits();
0312
0313   EBDetId Xtals5x5[25];
0314   for (unsigned int icry = 0; icry < 25; icry++) {
0315     unsigned int row = icry / 5;
0316     unsigned int column = icry % 5;
0317     int ieta = xtalInBeam_.ieta() + column - 2;
0318     int iphi = xtalInBeam_.iphi() + row - 2;
0319     EBDetId tempId(ieta, iphi, EBDetId::ETAPHIMODE);
0320     if (tempId.ism() == 1)
0321       Xtals5x5[icry] = tempId;
0322     else
0323       Xtals5x5[icry] = EBDetId(0);
0324     ///       std::cout << "** Xtal in the matrix **** row " << row  << ", column " << column << ", xtal " << Xtals5x5[icry].ic() << std::endl;
0325   }
0326
0327   double samples_save[10];
0328   for (int i = 0; i < 10; ++i)
0329     samples_save[i] = 0.0;
0330
0331   // find the rechit corresponding digi and the max sample
0332   EBDigiCollection::const_iterator myDg = digis->find(xtalInBeam_);
0333   double eMax = 0.;
0334   if (myDg != digis->end()) {
0335     EBDataFrame myDigi = (*myDg);
0336     for (int sample = 0; sample < myDigi.size(); ++sample) {
0337       double analogSample = myDigi.sample(sample).adc();
0338       samples_save[sample] = analogSample;
0339       //  std::cout << analogSample << " ";
0340       if (eMax < analogSample) {
0341         eMax = analogSample;
0342       }
0343     }
0344     // std::cout << std::endl;
0345   }
0346
0347   for (int i = 0; i < 10; ++i) {
0348     h_Shape_->Fill(double(i) + recTDC->offset(), samples_save[i]);
0349   }
0350
0351   double amplitude[25];
0352
0353   double amplitude3x3 = 0;
0354   double amplitude5x5 = 0;
0355
0356   for (unsigned int icry = 0; icry < 25; icry++) {
0357     if (!Xtals5x5[icry].null()) {
0358       amplitude[icry] = (hits->find(Xtals5x5[icry]))->amplitude();
0359       amplitude5x5 += amplitude[icry];
0360       // Is in 3x3?
0361       if (icry == 6 || icry == 7 || icry == 8 || icry == 11 || icry == 12 || icry == 13 || icry == 16 || icry == 17 ||
0362           icry == 18) {
0363         amplitude3x3 += amplitude[icry];
0364       }
0365     }
0366   }
0367
0368   h_e1x1->Fill(amplitude[12]);
0369   h_e3x3->Fill(amplitude3x3);
0370   h_e5x5->Fill(amplitude5x5);
0371
0372   h_e1e9->Fill(amplitude[12] / amplitude3x3);
0373   h_e1e25->Fill(amplitude[12] / amplitude5x5);
0374   h_e9e25->Fill(amplitude3x3 / amplitude5x5);
0375
0376   if (recTDC)
0377     h_ampltdc->Fill(recTDC->offset(), amplitude[12]);
0378
0379   if (recHodo) {
0380     float x = recHodo->posX();
0381     float y = recHodo->posY();
0382     float xslope = recHodo->slopeX();
0383     float yslope = recHodo->slopeY();
0384     float xqual = recHodo->qualX();
0385     float yqual = recHodo->qualY();
0386
0387     //Filling beam profiles
0388     h_bprofx->Fill(x);
0389     h_bprofy->Fill(y);
0390     h_qualx->Fill(xqual);
0391     h_qualy->Fill(yqual);
0392     h_slopex->Fill(xslope);
0393     h_slopey->Fill(yslope);
0394
0395     //Fill central events
0396
0397     if (fabs(x + 2.5) < 2.5 && fabs(y + 0.5) < 2.5) {
0398       h_e1x1_center->Fill(amplitude[12]);
0399       h_e3x3_center->Fill(amplitude3x3);
0400       h_e5x5_center->Fill(amplitude5x5);
0401
0402       h_e1e9->Fill(amplitude[12] / amplitude3x3);
0403       h_e1e25->Fill(amplitude[12] / amplitude5x5);
0404       h_e9e25->Fill(amplitude3x3 / amplitude5x5);
0405     }
0406
0407     for (unsigned int icry = 0; icry < 25; icry++) {
0408       h_mapx[icry]->Fill(x, amplitude[icry]);
0409       h_mapy[icry]->Fill(y, amplitude[icry]);
0410     }
0411
0412     h_e1e9_mapx->Fill(x, amplitude[12] / amplitude3x3);
0413     h_e1e9_mapy->Fill(y, amplitude[12] / amplitude3x3);
0414
0415     h_e1e25_mapx->Fill(x, amplitude[12] / amplitude5x5);
0416     h_e1e25_mapy->Fill(y, amplitude[12] / amplitude5x5);
0417
0418     h_e9e25_mapx->Fill(x, amplitude3x3 / amplitude5x5);
0419     h_e9e25_mapy->Fill(y, amplitude3x3 / amplitude5x5);
0420   }
0421 }