EcalSimpleTBAnalyzer.cc

Go to the documentation of this file.

//
//
//

#include "RecoTBCalo/EcalSimpleTBAnalysis/interface/EcalSimpleTBAnalyzer.h"
#include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRecInfo.h"
#include "TBDataFormats/EcalTBObjects/interface/EcalTBTDCRecInfo.h"
#include "TBDataFormats/EcalTBObjects/interface/EcalTBEventHeader.h"

//#include<fstream>

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"

#include <iostream>
#include <string>
#include <stdexcept>
//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//

//========================================================================
EcalSimpleTBAnalyzer::EcalSimpleTBAnalyzer(const edm::ParameterSet& iConfig)
    : xtalInBeam_(0)
//========================================================================
{
  //now do what ever initialization is needed
  rootfile_ = iConfig.getUntrackedParameter<std::string>("rootfile", "ecalSimpleTBanalysis.root");
  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
  hitCollection_ = iConfig.getParameter<std::string>("hitCollection");
  hitProducer_ = iConfig.getParameter<std::string>("hitProducer");
  hodoRecInfoCollection_ = iConfig.getParameter<std::string>("hodoRecInfoCollection");
  hodoRecInfoProducer_ = iConfig.getParameter<std::string>("hodoRecInfoProducer");
  tdcRecInfoCollection_ = iConfig.getParameter<std::string>("tdcRecInfoCollection");
  tdcRecInfoProducer_ = iConfig.getParameter<std::string>("tdcRecInfoProducer");
  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");

  std::cout << "EcalSimpleTBAnalyzer: fetching hitCollection: " << hitCollection_.c_str() << " produced by "
            << hitProducer_.c_str() << std::endl;
}

//========================================================================
EcalSimpleTBAnalyzer::~EcalSimpleTBAnalyzer()
//========================================================================
{
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  // Amplitude vs TDC offset
  //   if (h_ampltdc)
  //   delete h_ampltdc;

  //   // Reconstructed energies
  //   delete h_e1x1;
  //   delete h_e3x3;
  //   delete h_e5x5;

  //   delete h_bprofx;
  //   delete h_bprofy;

  //   delete h_qualx;
  //   delete h_qualy;

  //   delete h_slopex;
  //   delete h_slopey;

  //   delete h_mapx;
  //   delete h_mapy;
}

//========================================================================
void EcalSimpleTBAnalyzer::beginJob() {
  //========================================================================

  // Amplitude vs TDC offset
  h_ampltdc = new TH2F("h_ampltdc", "Max Amplitude vs TDC offset", 100, 0., 1., 1000, 0., 4000.);

  // Reconstructed energies
  h_tableIsMoving = new TH1F("h_tableIsMoving", "TableIsMoving", 100000, 0., 100000.);

  h_e1x1 = new TH1F("h_e1x1", "E1x1 energy", 1000, 0., 4000.);
  h_e3x3 = new TH1F("h_e3x3", "E3x3 energy", 1000, 0., 4000.);
  h_e5x5 = new TH1F("h_e5x5", "E5x5 energy", 1000, 0., 4000.);

  h_e1x1_center = new TH1F("h_e1x1_center", "E1x1 energy", 1000, 0., 4000.);
  h_e3x3_center = new TH1F("h_e3x3_center", "E3x3 energy", 1000, 0., 4000.);
  h_e5x5_center = new TH1F("h_e5x5_center", "E5x5 energy", 1000, 0., 4000.);

  h_e1e9 = new TH1F("h_e1e9", "E1/E9 ratio", 600, 0., 1.2);
  h_e1e25 = new TH1F("h_e1e25", "E1/E25 ratio", 600, 0., 1.2);
  h_e9e25 = new TH1F("h_e9e25", "E9/E25 ratio", 600, 0., 1.2);

  h_bprofx = new TH1F("h_bprofx", "Beam Profile X", 100, -20., 20.);
  h_bprofy = new TH1F("h_bprofy", "Beam Profile Y", 100, -20., 20.);

  h_qualx = new TH1F("h_qualx", "Beam Quality X", 5000, 0., 5.);
  h_qualy = new TH1F("h_qualy", "Beam Quality X", 5000, 0., 5.);

  h_slopex = new TH1F("h_slopex", "Beam Slope X", 500, -5e-4, 5e-4);
  h_slopey = new TH1F("h_slopey", "Beam Slope Y", 500, -5e-4, 5e-4);

  char hname[50];
  char htitle[50];
  for (unsigned int icry = 0; icry < 25; icry++) {
    sprintf(hname, "h_mapx_%d", icry);
    sprintf(htitle, "Max Amplitude vs X %d", icry);
    h_mapx[icry] = new TH2F(hname, htitle, 80, -20, 20, 1000, 0., 4000.);
    sprintf(hname, "h_mapy_%d", icry);
    sprintf(htitle, "Max Amplitude vs Y %d", icry);
    h_mapy[icry] = new TH2F(hname, htitle, 80, -20, 20, 1000, 0., 4000.);
  }

  h_e1e9_mapx = new TH2F("h_e1e9_mapx", "E1/E9 vs X", 80, -20, 20, 600, 0., 1.2);
  h_e1e9_mapy = new TH2F("h_e1e9_mapy", "E1/E9 vs Y", 80, -20, 20, 600, 0., 1.2);

  h_e1e25_mapx = new TH2F("h_e1e25_mapx", "E1/E25 vs X", 80, -20, 20, 600, 0., 1.2);
  h_e1e25_mapy = new TH2F("h_e1e25_mapy", "E1/E25 vs Y", 80, -20, 20, 600, 0., 1.2);

  h_e9e25_mapx = new TH2F("h_e9e25_mapx", "E9/E25 vs X", 80, -20, 20, 600, 0., 1.2);
  h_e9e25_mapy = new TH2F("h_e9e25_mapy", "E9/E25 vs Y", 80, -20, 20, 600, 0., 1.2);

  h_Shape_ = new TH2F("h_Shape_", "Xtal in Beam Shape", 250, 0, 10, 350, 0, 3500);
}

//========================================================================
void EcalSimpleTBAnalyzer::endJob() {
  //========================================================================

  TFile f(rootfile_.c_str(), "RECREATE");

  // Amplitude vs TDC offset
  h_ampltdc->Write();

  // Reconstructed energies
  h_e1x1->Write();
  h_e3x3->Write();
  h_e5x5->Write();

  h_e1x1_center->Write();
  h_e3x3_center->Write();
  h_e5x5_center->Write();

  h_e1e9->Write();
  h_e1e25->Write();
  h_e9e25->Write();

  h_bprofx->Write();
  h_bprofy->Write();

  h_qualx->Write();
  h_qualy->Write();

  h_slopex->Write();
  h_slopey->Write();

  h_Shape_->Write();

  for (unsigned int icry = 0; icry < 25; icry++) {
    h_mapx[icry]->Write();
    h_mapy[icry]->Write();
  }

  h_e1e9_mapx->Write();
  h_e1e9_mapy->Write();

  h_e1e25_mapx->Write();
  h_e1e25_mapy->Write();

  h_e9e25_mapx->Write();
  h_e9e25_mapy->Write();

  h_tableIsMoving->Write();

  f.Close();
}

//
// member functions
//

//========================================================================
void EcalSimpleTBAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  //========================================================================

  using namespace edm;
  using namespace cms;

  Handle<EBDigiCollection> pdigis;
  const EBDigiCollection* digis = nullptr;
  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
  iEvent.getByLabel(digiProducer_, digiCollection_, pdigis);
  if (pdigis.isValid()) {
    digis = pdigis.product();  // get a ptr to the product
    //iEvent.getByLabel( hitProducer_, phits);
  } else {
    edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << digiCollection_;
  }

  // fetch the digis and compute signal amplitude
  Handle<EBUncalibratedRecHitCollection> phits;
  const EBUncalibratedRecHitCollection* hits = nullptr;
  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
  iEvent.getByLabel(hitProducer_, hitCollection_, phits);
  if (phits.isValid()) {
    hits = phits.product();  // get a ptr to the product
    //iEvent.getByLabel( hitProducer_, phits);
  } else {
    edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hitCollection_;
  }

  Handle<EcalTBHodoscopeRecInfo> pHodo;
  const EcalTBHodoscopeRecInfo* recHodo = nullptr;
  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
  iEvent.getByLabel(hodoRecInfoProducer_, hodoRecInfoCollection_, pHodo);
  if (pHodo.isValid()) {
    recHodo = pHodo.product();  // get a ptr to the product
  } else {
    edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hodoRecInfoCollection_;
  }

  Handle<EcalTBTDCRecInfo> pTDC;
  const EcalTBTDCRecInfo* recTDC = nullptr;
  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
  iEvent.getByLabel(tdcRecInfoProducer_, tdcRecInfoCollection_, pTDC);
  if (pTDC.isValid()) {
    recTDC = pTDC.product();  // get a ptr to the product
  } else {
    edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << tdcRecInfoCollection_;
  }

  Handle<EcalTBEventHeader> pEventHeader;
  const EcalTBEventHeader* evtHeader = nullptr;
  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
  iEvent.getByLabel(eventHeaderProducer_, pEventHeader);
  if (pEventHeader.isValid()) {
    evtHeader = pEventHeader.product();  // get a ptr to the product
  } else {
    edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << eventHeaderProducer_;
  }

  if (!hits)
    return;

  if (!recTDC)
    return;

  if (!recHodo)
    return;

  if (!evtHeader)
    return;

  if (hits->empty())
    return;

  if (evtHeader->tableIsMoving())
    h_tableIsMoving->Fill(evtHeader->eventNumber());

  // Crystal hit by beam
  if (xtalInBeam_.null()) {
    xtalInBeam_ = EBDetId(1, evtHeader->crystalInBeam(), EBDetId::SMCRYSTALMODE);
    std::cout << "Xtal In Beam is " << xtalInBeam_.ic() << std::endl;
  } else if (xtalInBeam_ != EBDetId(1, evtHeader->crystalInBeam(), EBDetId::SMCRYSTALMODE))
    return;

  if (evtHeader->tableIsMoving())
    return;

  //    EBDetId maxHitId(0);
  //    float maxHit= -999999.;

  //    for(EBUncalibratedRecHitCollection::const_iterator ithit = hits->begin(); ithit != hits->end(); ++ithit)
  //      {
  //        if (ithit->amplitude()>=maxHit)
  //     {
  //       maxHit=ithit->amplitude();
  //       maxHitId=ithit->id();
  //     }

  //      }

  //    if (maxHitId==EBDetId(0))
  //      return;

  //   EBDetId maxHitId(1,704,EBDetId::SMCRYSTALMODE);

  //Find EBDetId in a 5x5 Matrix (to be substituted by the Selector code)
  // Something like
  // EBFixedWindowSelector<EcalUncalibratedRecHit> Simple5x5Matrix(hits,maxHitId,5,5);
  // std::vector<EcalUncalibratedRecHit> Energies5x5 = Simple5x5Matrix.getHits();

  EBDetId Xtals5x5[25];
  for (unsigned int icry = 0; icry < 25; icry++) {
    unsigned int row = icry / 5;
    unsigned int column = icry % 5;
    int ieta = xtalInBeam_.ieta() + column - 2;
    int iphi = xtalInBeam_.iphi() + row - 2;
    EBDetId tempId(ieta, iphi, EBDetId::ETAPHIMODE);
    if (tempId.ism() == 1)
      Xtals5x5[icry] = tempId;
    else
      Xtals5x5[icry] = EBDetId(0);
  }

  double samples_save[10];
  for (int i = 0; i < 10; ++i)
    samples_save[i] = 0.0;

  // find the rechit corresponding digi and the max sample
  EBDigiCollection::const_iterator myDg = digis->find(xtalInBeam_);
  double eMax = 0.;
  if (myDg != digis->end()) {
    EBDataFrame myDigi = (*myDg);
    for (int sample = 0; sample < myDigi.size(); ++sample) {
      double analogSample = myDigi.sample(sample).adc();
      samples_save[sample] = analogSample;
      //  std::cout << analogSample << " ";
      if (eMax < analogSample) {
        eMax = analogSample;
      }
    }
    // std::cout << std::endl;
  }

  for (int i = 0; i < 10; ++i) {
    h_Shape_->Fill(double(i) + recTDC->offset(), samples_save[i]);
  }

  double amplitude[25];

  double amplitude3x3 = 0;
  double amplitude5x5 = 0;

  for (unsigned int icry = 0; icry < 25; icry++) {
    if (!Xtals5x5[icry].null()) {
      amplitude[icry] = (hits->find(Xtals5x5[icry]))->amplitude();
      amplitude5x5 += amplitude[icry];
      // Is in 3x3?
      if (icry == 6 || icry == 7 || icry == 8 || icry == 11 || icry == 12 || icry == 13 || icry == 16 || icry == 17 ||
          icry == 18) {
        amplitude3x3 += amplitude[icry];
      }
    }
  }

  h_e1x1->Fill(amplitude[12]);
  h_e3x3->Fill(amplitude3x3);
  h_e5x5->Fill(amplitude5x5);

  h_e1e9->Fill(amplitude[12] / amplitude3x3);
  h_e1e25->Fill(amplitude[12] / amplitude5x5);
  h_e9e25->Fill(amplitude3x3 / amplitude5x5);

  if (recTDC)
    h_ampltdc->Fill(recTDC->offset(), amplitude[12]);

  if (recHodo) {
    float x = recHodo->posX();
    float y = recHodo->posY();
    float xslope = recHodo->slopeX();
    float yslope = recHodo->slopeY();
    float xqual = recHodo->qualX();
    float yqual = recHodo->qualY();

    //Filling beam profiles
    h_bprofx->Fill(x);
    h_bprofy->Fill(y);
    h_qualx->Fill(xqual);
    h_qualy->Fill(yqual);
    h_slopex->Fill(xslope);
    h_slopey->Fill(yslope);

    //Fill central events

    if (fabs(x + 2.5) < 2.5 && fabs(y + 0.5) < 2.5) {
      h_e1x1_center->Fill(amplitude[12]);
      h_e3x3_center->Fill(amplitude3x3);
      h_e5x5_center->Fill(amplitude5x5);

      h_e1e9->Fill(amplitude[12] / amplitude3x3);
      h_e1e25->Fill(amplitude[12] / amplitude5x5);
      h_e9e25->Fill(amplitude3x3 / amplitude5x5);
    }

    for (unsigned int icry = 0; icry < 25; icry++) {
      h_mapx[icry]->Fill(x, amplitude[icry]);
      h_mapy[icry]->Fill(y, amplitude[icry]);
    }

    h_e1e9_mapx->Fill(x, amplitude[12] / amplitude3x3);
    h_e1e9_mapy->Fill(y, amplitude[12] / amplitude3x3);

    h_e1e25_mapx->Fill(x, amplitude[12] / amplitude5x5);
    h_e1e25_mapy->Fill(y, amplitude[12] / amplitude5x5);

    h_e9e25_mapx->Fill(x, amplitude3x3 / amplitude5x5);
    h_e9e25_mapy->Fill(y, amplitude3x3 / amplitude5x5);
  }
}