ZeePlots.cc

Go to the documentation of this file.

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
#include "Calibration/Tools/interface/calibXMLwriter.h"
#include "Calibration/Tools/interface/CalibrationCluster.h"
#include "Calibration/Tools/interface/CalibElectron.h"
#include "Calibration/Tools/interface/HouseholderDecomposition.h"
#include "Calibration/Tools/interface/MinL3Algorithm.h"
#include "Calibration/EcalCalibAlgos/interface/ZeePlots.h"
#include "Calibration/EcalCalibAlgos/interface/ZeeKinematicTools.h"
#include "DataFormats/EgammaCandidates/interface/Electron.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TRandom.h"


#include <iostream>
#include <string>
#include <stdexcept>
#include <vector>

#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"

ZeePlots::ZeePlots( const char* fileName )
{

  fileName_ = fileName;
  file_ = new TFile(fileName_, "RECREATE");
}


ZeePlots::~ZeePlots()
{
  
  file_->Close();
  
  delete file_;

}

//========================================================================

void ZeePlots::openFile(){


  file_ -> cd();

}
//========================================================================

void ZeePlots::bookZMCHistograms(){

  file_ -> cd();

  h1_gen_ZMass_ = new TH1F("gen_ZMass","Generated Z mass",200,0.,150.);
  h1_gen_ZMass_->SetXTitle("gen_ZMass (GeV)");
  h1_gen_ZMass_->SetYTitle("events");

  h1_gen_ZEta_ = new TH1F("gen_ZEta","Eta of gen Z",200,-6.,6.);
  h1_gen_ZEta_->SetXTitle("#eta");
  h1_gen_ZEta_->SetYTitle("events");

  h1_gen_ZPhi_ = new TH1F("gen_ZPhi","Phi of gen Z",200,-4.,4.);
  h1_gen_ZPhi_->SetXTitle("#phi");
  h1_gen_ZPhi_->SetYTitle("events");

  h1_gen_ZRapidity_ = new TH1F("gen_ZRapidity","Rapidity of gen Z",200,-6.,6.);
  h1_gen_ZRapidity_->SetXTitle("Y");
  h1_gen_ZRapidity_->SetYTitle("events");

  h1_gen_ZPt_ = new TH1F("gen_ZPt","Pt of gen Z",200, 0.,100.);
  h1_gen_ZPt_->SetXTitle("p_{T} (GeV/c)");
  h1_gen_ZPt_->SetYTitle("events");


}

void ZeePlots::bookZHistograms(){

  file_ -> cd();

  h1_reco_ZEta_ = new TH1F("reco_ZEta","Eta of reco Z",200,-6.,6.);
  h1_reco_ZEta_->SetXTitle("#eta");
  h1_reco_ZEta_->SetYTitle("events");
  
  h1_reco_ZTheta_ = new TH1F("reco_ZTheta","Theta of reco Z",200, 0., 4.);
  h1_reco_ZTheta_->SetXTitle("#theta");
  h1_reco_ZTheta_->SetYTitle("events");
  
  h1_reco_ZRapidity_ = new TH1F("reco_ZRapidity","Rapidity of reco Z",200,-6.,6.);
  h1_reco_ZRapidity_->SetXTitle("Y");
  h1_reco_ZRapidity_->SetYTitle("events");
  
  h1_reco_ZPhi_ = new TH1F("reco_ZPhi","Phi of reco Z",100,-4.,4.);
  h1_reco_ZPhi_->SetXTitle("#phi");
  h1_reco_ZPhi_->SetYTitle("events");
  
  h1_reco_ZPt_ = new TH1F("reco_ZPt","Pt of reco Z",200,0.,100.);
  h1_reco_ZPt_->SetXTitle("p_{T} (GeV/c)");
  h1_reco_ZPt_->SetYTitle("events");
  
  
}

//========================================================================

void ZeePlots::fillZInfo( std::pair<calib::CalibElectron*,calib::CalibElectron*> myZeeCandidate ) {

  h1_reco_ZEta_->Fill( ZeeKinematicTools::calculateZEta(myZeeCandidate) );
  h1_reco_ZTheta_->Fill( ZeeKinematicTools::calculateZTheta(myZeeCandidate) );
  h1_reco_ZRapidity_->Fill( ZeeKinematicTools::calculateZRapidity(myZeeCandidate) );
  h1_reco_ZPhi_->Fill( ZeeKinematicTools::calculateZPhi(myZeeCandidate) );
  h1_reco_ZPt_->Fill( ZeeKinematicTools::calculateZPt(myZeeCandidate) );

}



//========================================================================

void ZeePlots::writeZHistograms() {

  file_->cd();

  h1_reco_ZEta_->Write();
  h1_reco_ZTheta_->Write();
  h1_reco_ZRapidity_->Write();
  h1_reco_ZPhi_->Write();
  h1_reco_ZPt_->Write();
  
}

//========================================================================

void ZeePlots::writeMCZHistograms() {

  file_->cd();
  
  h1_gen_ZRapidity_->Write();
  h1_gen_ZPt_->Write();
  h1_gen_ZPhi_->Write();

}

//========================================================================

void ZeePlots::fillZMCInfo( const HepMC::GenEvent* myGenEvent ) {

  file_->cd();

  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
    p != myGenEvent->particles_end(); ++p ) {//loop over MC particles
 
    if ( (*p)->pdg_id() == 23 && (*p)->status() == 2 ){
      
      h1_gen_ZMass_->Fill( (*p)->momentum().m() );
      h1_gen_ZEta_->Fill( (*p)->momentum().eta() );
      
      float genZ_Y = 0.5 * log ( ( (*p)->momentum().e() + (*p)->momentum().pz() ) /  ( (*p)->momentum().e() - (*p)->momentum().pz() ) )   ;
      
      h1_gen_ZRapidity_->Fill( genZ_Y );
      h1_gen_ZPt_->Fill((*p)->momentum().perp());
      h1_gen_ZPhi_->Fill((*p)->momentum().phi());

      

    }
  }//end loop over MC particles

  return;  
}

//========================================================================

void ZeePlots::bookEleMCHistograms(){

  file_->cd();

  h1_mcEle_Energy_ = new TH1F("mcEleEnergy","mc EleEnergy",300,0.,300.);
  h1_mcEle_Energy_->SetXTitle("E (GeV)");
  h1_mcEle_Energy_->SetYTitle("events");

  h1_mcElePt_ = new TH1F("mcElePt","p_{T} of MC electrons",300,0.,300.);
  h1_mcElePt_->SetXTitle("p_{T}(GeV/c)");
  h1_mcElePt_->SetYTitle("events");
  
  h1_mcEleEta_ = new TH1F("mcEleEta","Eta of MC electrons",100,-4.,4.);
  h1_mcEleEta_->SetXTitle("#eta");
  h1_mcEleEta_->SetYTitle("events");

  h1_mcElePhi_ = new TH1F("mcElePhi","Phi of MC electrons",100,-4.,4.);
  h1_mcElePhi_->SetXTitle("#phi");
  h1_mcElePhi_->SetYTitle("events");

}

//========================================================================

void ZeePlots::fillEleMCInfo( const HepMC::GenEvent* myGenEvent ) {

  file_->cd();
  
  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
    p != myGenEvent->particles_end(); ++p ) {
    
    if (  abs( (*p)->pdg_id() ) == 11 )
      {
    h1_mcEle_Energy_->Fill( (*p)->momentum().e() );
    h1_mcElePt_->Fill( (*p)->momentum().perp() );
    h1_mcEleEta_->Fill( (*p)->momentum().eta() );
    h1_mcElePhi_->Fill( (*p)->momentum().phi() );
    
      }//matches if (  abs( (*p)->pdg_id() ) == 11 )

  }//end loop over MC particles
  
}

//========================================================================
void ZeePlots::bookEleHistograms(){

  file_->cd();

  h1_nEleReco_ = new TH1F("h1_nEleReco", "h1_nEleReco", 20, 0, 20 );
  h1_nEleReco_->SetXTitle("Num. of reco electrons");
  
  h1_recoEleEnergy_ = new TH1F("recoEleEnergy","EleEnergy from SC",300,0.,300.);
  h1_recoEleEnergy_->SetXTitle("eleSCEnergy(GeV)");
  h1_recoEleEnergy_->SetYTitle("events");
  
  h1_recoElePt_ = new TH1F("recoElePt","p_{T} of reco electrons",300,0.,300.);
  h1_recoElePt_->SetXTitle("p_{T}(GeV/c)");
  h1_recoElePt_->SetYTitle("events");

  h1_recoEleEta_ = new TH1F("recoEleEta","Eta of reco electrons",100,-4.,4.);
  h1_recoEleEta_->SetXTitle("#eta");
  h1_recoEleEta_->SetYTitle("events");

   
  h1_recoElePhi_ = new TH1F("recoElePhi","Phi of reco electrons",100,-4.,4.);
  h1_recoElePhi_->SetXTitle("#phi");
  h1_recoElePhi_->SetYTitle("events");



}

//========================================================================

void ZeePlots::fillEleInfo(const reco::GsfElectronCollection* electronCollection) {

  file_->cd();

  h1_nEleReco_->Fill(electronCollection->size());
  
  for(reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();   eleIt != electronCollection->end(); eleIt++)
    {
      
  file_->cd();

      h1_recoEleEnergy_->Fill( eleIt->superCluster()->energy() );
      h1_recoElePt_->Fill( eleIt->pt() );
      h1_recoEleEta_->Fill( eleIt->eta() );
      h1_recoElePhi_->Fill( eleIt->phi() );
      
    }//end loop on electrons   

}

//========================================================================

void ZeePlots::writeEleHistograms(){

  file_->cd();

  std::cout << "Start with ZeePlots::writeEleHistograms(), done file_->cd(); " << std::endl;
 
  h1_recoEleEnergy_->Write();
  h1_recoElePt_->Write();
  h1_recoEleEta_->Write();
  h1_recoElePhi_->Write();

  std::cout << "Done with ZeePlots::writeEleHistograms() " << std::endl;

}

//========================================================================

void ZeePlots::writeMCEleHistograms(){

  file_->cd();

  std::cout << "Start with ZeePlots::writeMCEleHistograms(), done file_->cd(); " << std::endl;

  h1_mcEle_Energy_->Write();
  h1_mcElePt_->Write();
  h1_mcEleEta_->Write();
  h1_mcElePhi_->Write();

  std::cout << "Done with ZeePlots::writeMCEleHistograms() " << std::endl;

}

//========================================================================

void ZeePlots::bookHLTHistograms(){

  file_->cd();

  h1_FiredTriggers_= new TH1F("h1_FiredTriggers", "h1_FiredTriggers", 5,0,5);
 
  h1_HLTVisitedEvents_ = new TH1F("h1_HLTVisitedEvents", "h1_HLTVisitedEvents", 5,0,5);

  h1_HLT1Electron_FiredEvents_ =  new TH1F("h1_HLT1Electron_FiredEvents", "h1_HLT1Electron_FiredEvents", 5,0,5);
  h1_HLT2Electron_FiredEvents_ =  new TH1F("h1_HLT2Electron_FiredEvents", "h1_HLT2Electron_FiredEvents", 5,0,5);
  h1_HLT2ElectronRelaxed_FiredEvents_ =  new TH1F("h1_HLT2ElectronRelaxed_FiredEvents", "h1_HLT2ElectronRelaxed_FiredEvents", 5,0,5);

  h1_HLT1Electron_HLT2Electron_FiredEvents_ =  new TH1F("h1_HLT1Electron_HLT2Electron_FiredEvents", "h1_HLT1Electron_HLT2Electron_FiredEvents", 5,0,5);
  h1_HLT1Electron_HLT2ElectronRelaxed_FiredEvents_ =  new TH1F("h1_HLT1Electron_HLT2ElectronRelaxed_FiredEvents", "h1_HLT1Electron_HLT2ElectronRelaxed_FiredEvents", 5,0,5);
  h1_HLT2Electron_HLT2ElectronRelaxed_FiredEvents_ =  new TH1F("h1_HLT2Electron_HLT2ElectronRelaxed_FiredEvents", "h1_HLT2Electron_HLT2ElectronRelaxed_FiredEvents", 5,0,5);

  h1_HLT1Electron_HLT2Electron_HLT2ElectronRelaxed_FiredEvents_ =  new TH1F("h1_HLT1Electron_HLT2Electron_HLT2ElectronRelaxed_FiredEvents", "h1_HLT1Electron_HLT2Electron_HLT2ElectronRelaxed_FiredEvents", 5,0,5);


}


//========================================================================

void ZeePlots::fillHLTInfo( edm::Handle<edm::TriggerResults> hltTriggerResultHandle ){

  file_->cd();
  
  int hltCount = hltTriggerResultHandle->size();

  bool aHLTResults[200] = { false };
    
  for(int i = 0 ; i < hltCount ; i++) {

    aHLTResults[i] = hltTriggerResultHandle->accept(i);
    if(aHLTResults[i])
      h1_FiredTriggers_->Fill(i);

    //HLT bit 32 = HLT1Electron
    //HLT bit 34 = HLT2Electron
    //HLT bit 35 = HLT2ElectronRelaxed

  }

  h1_HLTVisitedEvents_->Fill(1);
    
  if(aHLTResults[32] && !aHLTResults[34] && !aHLTResults[35])
    h1_HLT1Electron_FiredEvents_->Fill(1);

  if(aHLTResults[34] && !aHLTResults[32] && !aHLTResults[35])
    h1_HLT2Electron_FiredEvents_->Fill(1);

  if(aHLTResults[35] && !aHLTResults[32] && !aHLTResults[34])
    h1_HLT2ElectronRelaxed_FiredEvents_->Fill(1);

  if(aHLTResults[32] && aHLTResults[34] && !aHLTResults[35])
    h1_HLT1Electron_HLT2Electron_FiredEvents_->Fill(1);

  if(aHLTResults[32] && aHLTResults[35] && !aHLTResults[34])
    h1_HLT1Electron_HLT2ElectronRelaxed_FiredEvents_->Fill(1);

  if(aHLTResults[34] && aHLTResults[35] && !aHLTResults[32])
    h1_HLT2Electron_HLT2ElectronRelaxed_FiredEvents_->Fill(1);

  if(aHLTResults[32] && aHLTResults[34] && aHLTResults[35])
    h1_HLT1Electron_HLT2Electron_HLT2ElectronRelaxed_FiredEvents_->Fill(1);



}


void ZeePlots::fillEleClassesPlots( calib::CalibElectron* myEle ){

  int myClass = myEle->getRecoElectron()->classification();

  float myEta = myEle->getRecoElectron()->eta();
  
  if(myClass==0 || myClass==100)
    h1_occupancyVsEtaGold_->Fill(myEta);
  
  std::cout<< "[ZeePlots::fillEleClassesPlots]Done gold"<< std::endl;
  
  if(myClass==40 || myClass==140)
    h1_occupancyVsEtaCrack_->Fill(myEta);
  
  std::cout<< "[ZeePlots::fillEleClassesPlots]Done crack"<< std::endl;
  
  if( (myClass>=30 && myClass<=34) || (myClass>=130 && myClass<=134) )
    h1_occupancyVsEtaShower_->Fill(myEta);
  
  std::cout<< "[ZeePlots::fillEleClassesPlots]Done shower"<< std::endl;
  
  if( myClass==10 || myClass==20 || myClass==110 || myClass ==120)
    h1_occupancyVsEtaSilver_->Fill(myEta);
  
  std::cout<< "[ZeePlots::fillEleClassesPlots]Done"<< std::endl;
  
}


void ZeePlots::bookEleClassesPlots(){

  file_->cd();

  h1_occupancyVsEtaGold_ = new TH1F("occupancyVsEtaGold","occupancyVsEtaGold", 200, -4.,4.);
  h1_occupancyVsEtaGold_->SetYTitle("Electron statistics");
  h1_occupancyVsEtaGold_->SetXTitle("Eta channel");

  h1_occupancyVsEtaSilver_ = new TH1F("occupancyVsEtaSilver","occupancyVsEtaSilver", 200, -4.,4.);
  h1_occupancyVsEtaSilver_->SetYTitle("Electron statistics");
  h1_occupancyVsEtaSilver_->SetXTitle("Eta channel");

  h1_occupancyVsEtaShower_ = new TH1F("occupancyVsEtaShower","occupancyVsEtaShower", 200, -4.,4.);
  h1_occupancyVsEtaShower_->SetYTitle("Electron statistics");
  h1_occupancyVsEtaShower_->SetXTitle("Eta channel");

  h1_occupancyVsEtaCrack_ = new TH1F("occupancyVsEtaCrack","occupancyVsEtaCrack", 200, -4.,4.);
  h1_occupancyVsEtaCrack_->SetYTitle("Electron statistics");
  h1_occupancyVsEtaCrack_->SetXTitle("Eta channel");

}

void ZeePlots::writeEleClassesPlots(){

  file_->cd();

  h1_occupancyVsEtaGold_->Write();
  h1_occupancyVsEtaSilver_->Write();
  h1_occupancyVsEtaShower_->Write();
  h1_occupancyVsEtaCrack_->Write();

}