EcalEBTrigPrimAnalyzer.cc

Go to the documentation of this file.

// system include files
#include <memory>
#include <utility>

// user include files
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"

#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"

#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

//#include "CalibCalorimetry/EcalTPGTools/interface/EcalEBTPGScale.h"
#include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"

#include "EcalEBTrigPrimAnalyzer.h"

#include <TMath.h>

using namespace edm;
class CaloSubdetectorGeometry;

EcalEBTrigPrimAnalyzer::EcalEBTrigPrimAnalyzer(const edm::ParameterSet&  iConfig)

{
  ecal_parts_.push_back("Barrel");
 
  histfile_=new TFile("histos.root","RECREATE");
  tree_ = new TTree("TPGtree","TPGtree");
  tree_->Branch("tpIphi",&tpIphi_,"tpIphi/I");
  tree_->Branch("tpIeta",&tpIeta_,"tpIeta/I");
  tree_->Branch("rhIphi",&rhIphi_,"rhIphi/I");
  tree_->Branch("rhIeta",&rhIeta_,"rhIeta/I");
  tree_->Branch("eRec",&eRec_,"eRec/F");
  tree_->Branch("tpgADC",&tpgADC_,"tpgADC/I");
  tree_->Branch("tpgGeV",&tpgGeV_,"tpgGeV/F");
  tree_->Branch("ttf",&ttf_,"ttf/I");
  tree_->Branch("fg",&fg_,"fg/I");
  for (unsigned int i=0;i<ecal_parts_.size();++i) {
    char title[30];
    sprintf(title,"%s_Et",ecal_parts_[i].c_str());
    ecal_et_[i]=new TH1I(title,"Et",255,0,255);
    sprintf(title,"%s_ttf",ecal_parts_[i].c_str());
    ecal_tt_[i]=new TH1I(title,"TTF",10,0,10);
    sprintf(title,"%s_fgvb",ecal_parts_[i].c_str());
    ecal_fgvb_[i]=new TH1I(title,"FGVB",10,0,10);
  }

  recHits_= iConfig.getParameter<bool>("AnalyzeRecHits");
  debug_= iConfig.getParameter<bool>("Debug");
  rechits_labelEB_=consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("inputRecHitsEB"));
  primToken_=consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("inputTP"));
  tokenEBdigi_=consumes<EBDigiCollection>(iConfig.getParameter<edm::InputTag>("barrelEcalDigis"));
  nEvents_=0;

  hTPvsTow_eta_= new TH2F("TP_vs_Tow_eta","TP vs Tow eta ; #eta(tow); #eta(tp)",50,-2.5,2.5,50,-2.5,2.5);
  if (recHits_) {
    hTPvsTow_ieta_= new TH2F("TP_vs_Tow_ieta","TP vs Tow ieta ; i#eta(tow); i#eta(tp)",200,-100,100,200,-100,100);

    hTPvsRechit_= new TH2F("TP_vs_RecHit","TP vs rechit Et;E_{T}(rh) (GeV);E_{T}(tp) (GeV)",100,0,50,100,0,50);
    hTPoverRechit_= new TH1F("TP_over_RecHit","Et(TP/rechit); E_{T}(tp)/E_{T}(rh); Counts",500,0,4);
    hRechitEt_= new TH1F("RecHitEt","E_{T};E_{T}(rh) (GeV);Counts",100,0,50);
    hTPEt_= new TH1F("TPEt","E_{T}{tp);E_{T}(rh) (GeV);Count",100,0,50);
    hRatioEt_ = new TH1F("RatioTPoverRH","Et",100,0,50);
    hAllRechitEt_= new TH1F("AllRecHit","Et",100,0,50);

    hTP_iphiVsieta_= new TH2F("TP_iphiVsieta","TP i#phi vs i#eta ; i#eta(tp); i#phi(tp)",10,70,80,10,340,350);
    hRH_iphiVsieta_= new TH2F("RH_iphiVsieta","RH i#phi vs i#eta ; i#eta(rh); i#phi(rh)",10,70,80,10,340,350);
    hTP_iphiVsieta_fullrange_= new TH2F("TP_iphiVsieta_fullrange","TP i#phi vs i#eta ; i#eta(tp); i#phi(tp)",200,-100,100,350,0,350);
    hRH_iphiVsieta_fullrange_= new TH2F("RH_iphiVsieta_fullrange","RH i#phi vs i#eta ; i#eta(rh); i#phi(rh)",200,-100,100,350,0,350);


  }
}


EcalEBTrigPrimAnalyzer::~EcalEBTrigPrimAnalyzer()
{

   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

  histfile_->Write();
  histfile_->Close();

}

void EcalEBTrigPrimAnalyzer::init(const edm::EventSetup & iSetup) {
  iSetup.get<IdealGeometryRecord>().get(eTTmap_);
}

//
// member functions
//

// ------------ method called to analyze the data  ------------
void
EcalEBTrigPrimAnalyzer::analyze(const edm::Event& iEvent, const  edm::EventSetup & iSetup)
{
  using namespace edm;
  using namespace std;
  nEvents_++;

  if ( nEvents_==1) this->init(iSetup);

  // Get input
  edm::Handle<EcalEBTrigPrimDigiCollection> tp;
  iEvent.getByToken(primToken_,tp);
  //
  /*
  edm::Handle<EBDigiCollection> barrelDigiHandle;
  const EBDigiCollection *ebdigi=NULL;
  iEvent.getByToken(tokenEBdigi_,barrelDigiHandle);
  ebdigi=barrelDigiHandle.product();
  */

  

  for (unsigned int i=0;i<tp.product()->size();i++) {
    EcalEBTriggerPrimitiveDigi d=(*(tp.product()))[i];
    int subdet=0;
    if (subdet==0) {
      ecal_et_[subdet]->Fill(d.compressedEt());
    }
  }
  if (!recHits_) return;

  // comparison with RecHits
  edm::Handle<EcalRecHitCollection> rechit_EB_col;
  iEvent.getByToken(rechits_labelEB_,rechit_EB_col);
 

  edm::ESHandle<CaloSubdetectorGeometry> theBarrelGeometry_handle;
  iSetup.get<EcalBarrelGeometryRecord>().get("EcalBarrel",theBarrelGeometry_handle);
  const CaloSubdetectorGeometry *theBarrelGeometry;
  theBarrelGeometry = &(*theBarrelGeometry_handle);

  if (debug_) std::cout << " TP analyzer  =================> Treating event  "<<iEvent.id()<< ", Number of EB rechits "<<  rechit_EB_col.product()->size() <<  " Number of TPs " <<  tp.product()->size() <<  std::endl;

  //if ( iEvent.id().event() != 648) return;

  //EcalEBTPGScale ecalScale ;
  EcalTPGScale ecalScale ;
  ecalScale.setEventSetup(iSetup) ;


  //  for(unsigned int iDigi = 0; iDigi < ebdigi->size() ; ++iDigi) {
  // EBDataFrame myFrame((*ebdigi)[iDigi]);  
  // const EBDetId & myId = myFrame.id();
     
   
    for (unsigned int i=0;i<tp.product()->size();i++) {
      EcalEBTriggerPrimitiveDigi d=(*(tp.product()))[i];
      const EBDetId TPid= d.id();
      // if ( myId != TPid ) continue;
   

      /*    
      int index=getIndex(ebdigi,coarser);
      std::cout << " Same xTal " << myId << " " << TPid << " coarser " << coarser << " index " << index << std::endl;
      double Et = ecalScale.getTPGInGeV(d.compressedEt(), coarser) ; 
      */
      float Et=d.compressedEt()*0.5;
      if ( Et==0) continue;

      //      EcalTrigTowerDetId coarser=(*eTTmap_).towerOf(myId);
      // does not work float etaTow =  theBarrelGeometry->getGeometry(coarser)->getPosition().theta();
      // float etaTP =  theBarrelGeometry->getGeometry(TPid)->getPosition().eta();
      // does not work hTPvsTow_eta_->Fill ( etaTow,  etaTP );
      //      hTPvsTow_ieta_->Fill ( coarser.ieta(),  TPid.ieta() );


      tpIphi_ = TPid.iphi() ;
      tpIeta_ = TPid.ieta() ;
      tpgADC_ = d.compressedEt();
      tpgGeV_ = Et ;

      hTP_iphiVsieta_->Fill(TPid.ieta(), TPid.iphi(), Et);
      hTP_iphiVsieta_fullrange_->Fill(TPid.ieta(), TPid.iphi(), Et);
      
      for (unsigned int j=0;j<rechit_EB_col.product()->size();j++) {
    const EBDetId & myid1=(*rechit_EB_col.product())[j].id();
    float theta =  theBarrelGeometry->getGeometry(myid1)->getPosition().theta();
    float rhEt=((*rechit_EB_col.product())[j].energy())*sin(theta);
    if ( myid1 == TPid ) {
      if (debug_) std::cout << " Analyzer same cristal " << myid1 << " " << TPid << std::endl;
      //      if ( rhEt < 1.5 && Et > 10 )  {
      // std::cout << " TP analyzer  =================> Treating event  "<<iEvent.id()<< ", Number of EB rechits "<<  rechit_EB_col.product()->size() <<  " Number of TPs " <<  tp.product()->size() <<  std::endl;
      //std::cout << " TP compressed et " << d.compressedEt()  << " Et in GeV  " <<  Et << " RH Et " << rhEt << " Et/rhEt " << Et/rhEt << std::endl;          
      //} 

      //std::cout << " TP out " <<  d << std::endl;

      //      for (int isam=0;isam< d.size();++isam) {
      // std::cout << " d[isam].raw() "  <<  d[isam].raw() << std::endl;
      //}

      rhIphi_ =  myid1.iphi() ;
      rhIeta_ =  myid1.ieta() ;
      hRH_iphiVsieta_->Fill(myid1.ieta(), myid1.iphi(), rhEt);
      hRH_iphiVsieta_fullrange_->Fill(myid1.ieta(), myid1.iphi(), rhEt);
      
      hTPvsRechit_->Fill(rhEt,Et);
      hTPoverRechit_->Fill(Et/rhEt);
      if (debug_) std::cout << " TP compressed et " << d.compressedEt()  << " Et in GeV  " <<  Et << " RH Et " << rhEt << " Et/rhEt " << Et/rhEt << std::endl;          
      hRechitEt_->Fill(rhEt);
      hTPEt_->Fill(Et);
      if ( rhEt < 1000000) eRec_ = rhEt;
      tree_->Fill() ;
    }

      }
      
    }   // end loop over TP collection
    
    //  } // end loop over digi collection
  
  hRatioEt_->Divide( hTPEt_, hRechitEt_);

  for (unsigned int j=0;j<rechit_EB_col.product()->size();j++) {
    const EBDetId & myid1=(*rechit_EB_col.product())[j].id();
    float theta =  theBarrelGeometry->getGeometry(myid1)->getPosition().theta();
    float rhEt=((*rechit_EB_col.product())[j].energy())*sin(theta);
    if ( rhEt >0 ) 
      hAllRechitEt_ ->Fill(rhEt);
  }



}

void
EcalEBTrigPrimAnalyzer::endJob(){
  for (unsigned int i=0;i<ecal_parts_.size();++i) {
    ecal_et_[i]->Write();
    ecal_tt_[i]->Write();
    ecal_fgvb_[i]->Write();
  }


  if (recHits_) {
    hTPvsTow_ieta_->Write();
    hTPvsTow_eta_->Write();    
    hTPvsRechit_->Write();
    hTPoverRechit_->Write();
    hAllRechitEt_->Write();
    hRechitEt_->Write();
    hTPEt_->Write();
    hRatioEt_->Write();
    hTP_iphiVsieta_->Write();
    hRH_iphiVsieta_->Write();
    hTP_iphiVsieta_fullrange_->Write();
    hRH_iphiVsieta_fullrange_->Write();

  }
}