EcalRecHits/src/EcalTBValidation.cc

0001 // Ecal H4 tesbeam analysis
0002 #include "DQMServices/Core/interface/DQMStore.h"
0003 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0004 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0005 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "TBDataFormats/EcalTBObjects/interface/EcalTBEventHeader.h"
0008 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRecInfo.h"
0009 #include "TBDataFormats/EcalTBObjects/interface/EcalTBTDCRecInfo.h"
0010 #include "Validation/EcalRecHits/interface/EcalTBValidation.h"
0011
0012 #include <iostream>
0013 #include <string>
0014
0015 EcalTBValidation::EcalTBValidation(const edm::ParameterSet &config) {
0016   data_ = config.getUntrackedParameter<int>("data", -1000);
0017   xtalInBeam_ = config.getUntrackedParameter<int>("xtalInBeam", -1000);
0018
0019   digiCollection_ = config.getParameter<std::string>("digiCollection");
0020   digiProducer_ = config.getParameter<std::string>("digiProducer");
0021   hitCollection_ = config.getParameter<std::string>("hitCollection");
0022   hitProducer_ = config.getParameter<std::string>("hitProducer");
0023   hodoRecInfoCollection_ = config.getParameter<std::string>("hodoRecInfoCollection");
0024   hodoRecInfoProducer_ = config.getParameter<std::string>("hodoRecInfoProducer");
0025   tdcRecInfoCollection_ = config.getParameter<std::string>("tdcRecInfoCollection");
0026   tdcRecInfoProducer_ = config.getParameter<std::string>("tdcRecInfoProducer");
0027   eventHeaderCollection_ = config.getParameter<std::string>("eventHeaderCollection");
0028   eventHeaderProducer_ = config.getParameter<std::string>("eventHeaderProducer");
0029
0030   digi_Token_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0031   hit_Token_ = consumes<EBUncalibratedRecHitCollection>(edm::InputTag(hitProducer_, hitCollection_));
0032   hodoRec_Token_ = consumes<EcalTBHodoscopeRecInfo>(edm::InputTag(hodoRecInfoProducer_, hodoRecInfoCollection_));
0033   tdcRec_Token_ = consumes<EcalTBTDCRecInfo>(edm::InputTag(tdcRecInfoProducer_, tdcRecInfoCollection_));
0034   eventHeader_Token_ = consumes<EcalTBEventHeader>(edm::InputTag(eventHeaderProducer_));
0035   // rootfile_              =
0036   // config.getUntrackedParameter<std::string>("rootfile","EcalTBValidation.root");
0037
0038   // verbosity...
0039   verbose_ = config.getUntrackedParameter<bool>("verbose", false);
0040
0041   meETBxib_ = nullptr;
0042   meETBampltdc_ = nullptr;
0043   meETBShape_ = nullptr;
0044   meETBhodoX_ = nullptr;
0045   meETBhodoY_ = nullptr;
0046   meETBe1x1_ = nullptr;
0047   meETBe3x3_ = nullptr;
0048   meETBe5x5_ = nullptr;
0049   meETBe1e9_ = nullptr;
0050   meETBe1e25_ = nullptr;
0051   meETBe9e25_ = nullptr;
0052   meETBe1x1_center_ = nullptr;
0053   meETBe3x3_center_ = nullptr;
0054   meETBe5x5_center_ = nullptr;
0055   meETBe1vsX_ = nullptr;
0056   meETBe1vsY_ = nullptr;
0057   meETBe1e9vsX_ = nullptr;
0058   meETBe1e9vsY_ = nullptr;
0059   meETBe1e25vsX_ = nullptr;
0060   meETBe1e25vsY_ = nullptr;
0061   meETBe9e25vsX_ = nullptr;
0062   meETBe9e25vsY_ = nullptr;
0063 }
0064
0065 EcalTBValidation::~EcalTBValidation() {}
0066
0067 void EcalTBValidation::bookHistograms(DQMStore::IBooker &ibooker, edm::Run const &, edm::EventSetup const &) {
0068   std::string hname;
0069   ibooker.setCurrentFolder("EcalRecHitsV/EcalTBValidationTask");
0070
0071   hname = "xtal in beam position";
0072   meETBxib_ = ibooker.book2D(hname, hname, 85, 0., 85., 20, 0., 20.);
0073   hname = "Max Amplitude vs TDC offset";
0074   meETBampltdc_ = ibooker.book2D(hname, hname, 100, 0., 1., 1000, 0., 4000.);
0075   hname = "Beam Profile X";
0076   meETBhodoX_ = ibooker.book1D(hname, hname, 100, -20., 20.);
0077   hname = "Beam Profile Y";
0078   meETBhodoY_ = ibooker.book1D(hname, hname, 100, -20., 20.);
0079   hname = "E1x1 energy";
0080   meETBe1x1_ = ibooker.book1D(hname, hname, 1000, 0., 4000.);
0081   hname = "E3x3 energy";
0082   meETBe3x3_ = ibooker.book1D(hname, hname, 1000, 0., 4000.);
0083   hname = "E5x5 energy";
0084   meETBe5x5_ = ibooker.book1D(hname, hname, 1000, 0., 4000.);
0085   hname = "E1x1 energy center";
0086   meETBe1x1_center_ = ibooker.book1D(hname, hname, 1000, 0., 4000.);
0087   hname = "E3x3 energy center";
0088   meETBe3x3_center_ = ibooker.book1D(hname, hname, 1000, 0., 4000.);
0089   hname = "E5x5 energy center";
0090   meETBe5x5_center_ = ibooker.book1D(hname, hname, 1000, 0., 4000.);
0091   hname = "E1 over E9 ratio";
0092   meETBe1e9_ = ibooker.book1D(hname, hname, 600, 0., 1.2);
0093   hname = "E1 over E25 ratio";
0094   meETBe1e25_ = ibooker.book1D(hname, hname, 600, 0., 1.2);
0095   hname = "E9 over E25 ratio";
0096   meETBe9e25_ = ibooker.book1D(hname, hname, 600, 0., 1.2);
0097   hname = "E1 vs X";
0098   meETBe1vsX_ = ibooker.book2D(hname, hname, 80, -20, 20, 1000, 0., 4000.);
0099   hname = "E1 vs Y";
0100   meETBe1vsY_ = ibooker.book2D(hname, hname, 80, -20, 20, 1000, 0., 4000.);
0101   hname = "E1 over E9 vs X";
0102   meETBe1e9vsX_ = ibooker.book2D(hname, hname, 80, -20, 20, 600, 0., 1.2);
0103   hname = "E1 over E9 vs Y";
0104   meETBe1e9vsY_ = ibooker.book2D(hname, hname, 80, -20, 20, 600, 0., 1.2);
0105   hname = "E1 over E25 vs X";
0106   meETBe1e25vsX_ = ibooker.book2D(hname, hname, 80, -20, 20, 600, 0., 1.2);
0107   hname = "E1 over E25 vs Y";
0108   meETBe1e25vsY_ = ibooker.book2D(hname, hname, 80, -20, 20, 600, 0., 1.2);
0109   hname = "E9 over E25 vs X";
0110   meETBe9e25vsX_ = ibooker.book2D(hname, hname, 80, -20, 20, 600, 0., 1.2);
0111   hname = "E9 over E25 vs Y";
0112   meETBe9e25vsY_ = ibooker.book2D(hname, hname, 80, -20, 20, 600, 0., 1.2);
0113   hname = "Xtal in Beam Shape";
0114   meETBShape_ = ibooker.book2D(hname, hname, 250, 0, 10, 350, 0, 3500);
0115 }
0116
0117 void EcalTBValidation::analyze(const edm::Event &event, const edm::EventSetup &setup) {
0118   using namespace edm;
0119   using namespace cms;
0120
0121   // digis
0122   const EBDigiCollection *theDigis = nullptr;
0123   Handle<EBDigiCollection> pdigis;
0124   event.getByToken(digi_Token_, pdigis);
0125   if (pdigis.isValid()) {
0126     theDigis = pdigis.product();
0127   } else {
0128     std::cerr << "Error! can't get the product " << digiCollection_.c_str() << std::endl;
0129     return;
0130   }
0131
0132   // rechits
0133   const EBUncalibratedRecHitCollection *theHits = nullptr;
0134   Handle<EBUncalibratedRecHitCollection> phits;
0135   event.getByToken(hit_Token_, phits);
0136   if (phits.isValid()) {
0137     theHits = phits.product();
0138   } else {
0139     std::cerr << "Error! can't get the product " << hitCollection_.c_str() << std::endl;
0140     return;
0141   }
0142
0143   // hodoscopes
0144   const EcalTBHodoscopeRecInfo *theHodo = nullptr;
0145   Handle<EcalTBHodoscopeRecInfo> pHodo;
0146   event.getByToken(hodoRec_Token_, pHodo);
0147   if (pHodo.isValid()) {
0148     theHodo = pHodo.product();
0149   } else {
0150     std::cerr << "Error! can't get the product " << hodoRecInfoCollection_.c_str() << std::endl;
0151     return;
0152   }
0153
0154   // tdc
0155   const EcalTBTDCRecInfo *theTDC = nullptr;
0156   Handle<EcalTBTDCRecInfo> pTDC;
0157   event.getByToken(tdcRec_Token_, pTDC);
0158   if (pTDC.isValid()) {
0159     theTDC = pTDC.product();
0160   } else {
0161     std::cerr << "Error! can't get the product " << tdcRecInfoCollection_.c_str() << std::endl;
0162     return;
0163   }
0164
0165   // event header
0166   const EcalTBEventHeader *evtHeader = nullptr;
0167   Handle<EcalTBEventHeader> pEventHeader;
0168   event.getByToken(eventHeader_Token_, pEventHeader);
0169   if (pEventHeader.isValid()) {
0170     evtHeader = pEventHeader.product();
0171   } else {
0172     std::cerr << "Error! can't get the product " << eventHeaderProducer_.c_str() << std::endl;
0173     return;
0174   }
0175
0176   // -----------------------------------------------------------------------
0177   // xtal-in-beam
0178   EBDetId xtalInBeamId(1, xtalInBeam_, EBDetId::SMCRYSTALMODE);
0179   if (xtalInBeamId == EBDetId(0)) {
0180     return;
0181   }
0182   int xibEta = xtalInBeamId.ieta();
0183   int xibPhi = xtalInBeamId.iphi();
0184
0185   // skipping events with moving table (in data)
0186   if (data_ && (evtHeader->tableIsMoving()))
0187     return;
0188
0189   // amplitudes
0190   EBDetId Xtals5x5[25];
0191   for (unsigned int icry = 0; icry < 25; icry++) {
0192     unsigned int row = icry / 5;
0193     unsigned int column = icry % 5;
0194     int ieta = xtalInBeamId.ieta() + column - 2;
0195     int iphi = xtalInBeamId.iphi() + row - 2;
0196     if (EBDetId::validDetId(ieta, iphi)) {
0197       EBDetId tempId(ieta, iphi, EBDetId::ETAPHIMODE);
0198       if (tempId.ism() == 1)
0199         Xtals5x5[icry] = tempId;
0200       else
0201         Xtals5x5[icry] = EBDetId(0);
0202     } else {
0203       Xtals5x5[icry] = EBDetId(0);
0204     }
0205   }
0206
0207   // matrices
0208   double ampl1x1 = 0.;
0209   double ampl3x3 = 0.;
0210   double ampl5x5 = 0.;
0211   for (unsigned int icry = 0; icry < 25; icry++) {
0212     if (!Xtals5x5[icry].null()) {
0213       double theAmpl = (theHits->find(Xtals5x5[icry]))->amplitude();
0214       ampl5x5 += theAmpl;
0215       if (icry == 12) {
0216         ampl1x1 = theAmpl;
0217       }
0218       if (icry == 6 || icry == 7 || icry == 8 || icry == 11 || icry == 12 || icry == 13 || icry == 16 || icry == 17 ||
0219           icry == 18) {
0220         ampl3x3 += theAmpl;
0221       }
0222     }
0223   }
0224
0225   // pulse shape
0226   double sampleSave[10];
0227   for (int ii = 0; ii < 10; ++ii) {
0228     sampleSave[ii] = 0.0;
0229   }
0230   EBDigiCollection::const_iterator thisDigi = theDigis->find(xtalInBeamId);
0231   // int sMax = -1; // UNUSED
0232   double eMax = 0.;
0233   if (thisDigi != theDigis->end()) {
0234     EBDataFrame myDigi = (*thisDigi);
0235     for (int sample = 0; sample < myDigi.size(); ++sample) {
0236       double analogSample = myDigi.sample(sample).adc();
0237       sampleSave[sample] = analogSample;
0238       if (eMax < analogSample) {
0239         eMax = analogSample;
0240         // sMax = sample; // UNUSED
0241       }
0242     }
0243   }
0244
0245   // beam profile
0246   double xBeam = theHodo->posX();
0247   double yBeam = theHodo->posY();
0248
0249   // filling histos
0250   meETBxib_->Fill(xibEta, xibPhi);
0251   meETBhodoX_->Fill(xBeam);
0252   meETBhodoY_->Fill(yBeam);
0253   meETBampltdc_->Fill(theTDC->offset(), ampl1x1);
0254   meETBe1x1_->Fill(ampl1x1);
0255   meETBe3x3_->Fill(ampl3x3);
0256   meETBe5x5_->Fill(ampl5x5);
0257   meETBe1e9_->Fill(ampl1x1 / ampl3x3);
0258   meETBe1e25_->Fill(ampl1x1 / ampl5x5);
0259   meETBe9e25_->Fill(ampl3x3 / ampl5x5);
0260   meETBe1vsX_->Fill(xBeam, ampl1x1);
0261   meETBe1vsY_->Fill(yBeam, ampl1x1);
0262   meETBe1e9vsX_->Fill(xBeam, ampl1x1 / ampl3x3);
0263   meETBe1e9vsY_->Fill(yBeam, ampl1x1 / ampl3x3);
0264   meETBe1e25vsX_->Fill(xBeam, ampl1x1 / ampl5x5);
0265   meETBe1e25vsY_->Fill(yBeam, ampl1x1 / ampl5x5);
0266   meETBe9e25vsX_->Fill(xBeam, ampl3x3 / ampl5x5);
0267   meETBe9e25vsY_->Fill(yBeam, ampl3x3 / ampl5x5);
0268
0269   for (int ii = 0; ii < 10; ++ii) {
0270     meETBShape_->Fill(double(ii) + theTDC->offset(), sampleSave[ii]);
0271   }
0272
0273   if ((fabs(xBeam) < 2.5) && (fabs(yBeam) < 2.5)) {
0274     meETBe1x1_center_->Fill(ampl1x1);
0275     meETBe3x3_center_->Fill(ampl3x3);
0276     meETBe5x5_center_->Fill(ampl5x5);
0277   }
0278 }