EcalSimpleTBAnalysis/src/EcalSimple2007H4TBAnalyzer.cc

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