/Users/jmonk/Physics/ForIA/src/AnalysisTools/MinBiasTrackEfficiency.cxx

Go to the documentation of this file.

#include "ForIA/AnalysisTools/MinBiasTrackEfficiency.hh"
#include "ForIA/IDataLoader.hh"
#include "ForIA/Event.hh"
#include "ForIA/Vertex.hh"
#include "ForIA/Units.hh"

#include "TH1D.h"

#include <boost/numeric/conversion/cast.hpp>
#include <boost/lexical_cast.hpp>

#include <limits>

namespace ForIA{
 
  MinBiasTrackEfficiency::MinBiasTrackEfficiency(): m_doInitialise(true), m_doVtx(true), m_doTrigger(true), 
  m_canChangeSettings(true), m_nextEvent(true), m_vtx_z(0.), m_gotTruth(false){

  }
  void MinBiasTrackEfficiency::initialise(){
    
    DataLoaderPtr dl = IDataLoader::create();
    
    dl->open("Root/MinBias/Event7TeV_v20.root");
    
    TH1D *VertexEfficiency_NPart = dl->retrieve<TH1D>("VertexEfficiency_NPart");
    
    double nBins = VertexEfficiency_NPart->GetNbinsX();
    
    double dBin = VertexEfficiency_NPart->GetBinCenter(1) + 0.5;
    m_vtxNTrkOffset = boost::numeric_cast<int>(dBin) - 1;
    
    m_vtxEfficiencyVsNTrk.clear();
    
    for(int bin = 1; bin <= nBins; ++bin){
      double eff = VertexEfficiency_NPart->GetBinContent(bin);
      m_vtxEfficiencyVsNTrk.push_back(eff);
    }
    
    TH1D *VertEff_Z0min_p_BS2_Pt_lt_02 = dl->retrieve<TH1D>("VertEff_Z0min_p_BS2_Pt_lt_02");
    fillMap(VertEff_Z0min_p_BS2_Pt_lt_02, &m_vtxEfficiencyVsZ0_2TrkLt200);
    
    TH1D *VertEff_Z0min_p_BS2_Pt_gt_02 = dl->retrieve<TH1D>("VertEff_Z0min_p_BS2_Pt_gt_02");
    fillMap(VertEff_Z0min_p_BS2_Pt_gt_02, &m_vtxEfficiencyVsZ0_2TrkGt200);
    
    TH1D *VertEff_Z0min_p_BS3 = dl->retrieve<TH1D>("VertEff_Z0min_p_BS3");
    fillMap(VertEff_Z0min_p_BS3, &m_vtxEfficiencyVsZ0_3Trk);
    
    TH1D *VertEff_Z0min_p_BS4 = dl->retrieve<TH1D>("VertEff_Z0min_p_BS4");
    fillMap(VertEff_Z0min_p_BS4, &m_vtxEfficiencyVsZ0_4Trk);
    
    TH1D *TriggerEfficiency_NPart = dl->retrieve<TH1D>("TriggerEfficiency_NPart");
    
    nBins = TriggerEfficiency_NPart->GetNbinsX();
    m_trigEfficiencyVsNTrk.clear();
    for(int bin = 1; bin <= nBins; ++bin){
      double eff = TriggerEfficiency_NPart->GetBinContent(bin);
      m_trigEfficiencyVsNTrk.push_back(eff);
    }
    
    TH1D *VertexWeight_ZPos_Simu = dl->retrieve<TH1D>("VertexWeight_ZPos_Simu");
    fillMap(VertexWeight_ZPos_Simu, &m_vtxWeightVsZ0Truth);
    
    dl->close();
    m_doInitialise = false;
    
    m_trackSelection.setSelection(TrackSelection::MINBIAS_2_BS);
    
    return;
  }
  
  void MinBiasTrackEfficiency::fillMap(const TH1 *histo, map<double, double> *toMap){
    
    toMap->clear();
    
    int nBins = histo->GetNbinsX();
    
    for(int bin=1; bin <= nBins; ++bin){
      double eff = histo->GetBinContent(bin);
      double rightEdge = histo->GetXaxis()->GetBinUpEdge(bin);
      (*toMap)[rightEdge] = eff;
    }
    
    return;
  }
  
  bool MinBiasTrackEfficiency::setEvent(const Event &evt){
        
    if(m_doInitialise) initialise();
    m_canChangeSettings = false;
    m_nextEvent = true;
    m_trackSelection.setEvent(evt);
    
    VertexVector vertices= evt.truthVertices();
    if(vertices.size() != 0){
      m_gotTruth = true;
      m_vtx_z = vertices[0]->z();
    }else{
      m_gotTruth = false;
    }
    
    return true;
  }
  
  double MinBiasTrackEfficiency::eventWeight(){
    
    double eff = eventEfficiency();
    
    // test for zero event efficiency
    // Since the efficiency histos are quantised there is a minimum possible 
    // non-zero efficiency, which is well above this value
    // Of course we should never see a real event with 0 reconstruction efficiency
    // if we do, it indicates the efficiency estimates are incorrect.
    // In MC we may see an event with zero efficiency, depending on other cuts 
    // (e.g. a generated event with only 1 particle)
    if(eff < 0.000001){
            
      return 0.0;
    }
    if(m_gotTruth) return mcZWeight() / eff;
    return 1.0 / eff;
  }
  
  double MinBiasTrackEfficiency::eventEfficiency(){
        
    double eff = 1.0;
    
    if(m_doTrigger) eff *= triggerEfficiency();
    if(m_doVtx)     eff *= vertexEfficiency();
    
    return eff;
  }
  
  double MinBiasTrackEfficiency::triggerEfficiency(){
    
    TrackVector tracks = m_trackSelection.tracks();
    unsigned int nTracks = tracks.size();
  
    double eff = 1.0;
    
    if(nTracks < m_trigEfficiencyVsNTrk.size()){
      eff = m_trigEfficiencyVsNTrk[nTracks];
    }
        
    return eff;
  }
  
  double MinBiasTrackEfficiency::mcZWeight(){
    if(!m_gotTruth) return 1.0;
    
    map<double, double>::const_iterator it = m_vtxWeightVsZ0Truth.upper_bound(m_vtx_z);
    if(it == m_vtxWeightVsZ0Truth.end()) return 0.;
    return it->second;
  }
  
  double MinBiasTrackEfficiency::vertexEfficiency(){
    
    TrackVector tracks = m_trackSelection.tracks();
    unsigned int nTracks = tracks.size();
    
    double eff = 1.;
    
    switch(nTracks){
        
      case 0: 
      case 1:
        eff = 0.;
        break;
        
      case 2: 
      {
        double minPT = (tracks[0]->PT() < tracks[1]->PT())? tracks[0]->PT(): tracks[1]->PT();
        double deltaZ0 = fabs(tracks[0]->z0(Track::BEAM_SPOT) - tracks[1]->z0(Track::BEAM_SPOT));
        
        map<double, double>::const_iterator it;
        
        if(minPT < 200. * MeV){
          it = m_vtxEfficiencyVsZ0_2TrkLt200.upper_bound(deltaZ0);
          if(it == m_vtxEfficiencyVsZ0_2TrkLt200.end()){
            eff = 0.;
          }else{
            eff = it->second;
          }
        }else{
          it = m_vtxEfficiencyVsZ0_2TrkGt200.upper_bound(deltaZ0);
          if(it == m_vtxEfficiencyVsZ0_2TrkGt200.end()){
            eff = 0.;
          }else{
            eff = it->second;
          }
        }
      }
        break;
        
      case 3:
      {
        double deltaZ0 = minDeltaZ();
        map<double, double>::const_iterator bin = m_vtxEfficiencyVsZ0_3Trk.upper_bound(deltaZ0);
        if(bin==m_vtxEfficiencyVsZ0_3Trk.end()){
          eff = 0.;
        }else{
          eff = bin->second;
        }
      } 
        break;
      
      case 4:
      {
        double deltaZ0 = minDeltaZ();
        map<double, double>::const_iterator bin = m_vtxEfficiencyVsZ0_4Trk.upper_bound(deltaZ0);
        if(bin==m_vtxEfficiencyVsZ0_3Trk.end()){
          eff = 0.;
        }else{
          eff = bin->second;
        }        
      } 
        break;

      default:
        
        int bin = nTracks - m_vtxNTrkOffset;
        if(bin < 0){
          string er = "The vertex efficiency Vs. N tracks starts at " + 
          boost::lexical_cast<string>(m_vtxNTrkOffset) + ", but we are trying an event with " + 
          boost::lexical_cast<string>(nTracks) + " tracks!!";
          throw std::runtime_error(er);
        }
        if(bin < (int)m_vtxEfficiencyVsNTrk.size()){
          eff = m_vtxEfficiencyVsNTrk[bin];
        }
        
    }
           
    return eff;
  }
  
  double MinBiasTrackEfficiency::minDeltaZ(){
    double dz = std::numeric_limits<double>::max();
    
    TrackVector tracks = m_trackSelection.tracks();
    
    for(TrackVector::const_iterator track1 = tracks.begin();
        track1 != tracks.end(); ++track1){
      
      TrackVector::const_iterator track2 = track1;
      ++track2;
      while(track2 != tracks.end()){
        double tmp = fabs((*track1)->z0(Track::BEAM_SPOT) - (*track2)->z0(Track::BEAM_SPOT));
        if(tmp < dz) dz = tmp;
        ++track2;
      }
      
    }
    return dz;
  }
  
  
  bool MinBiasTrackEfficiency::setUseVtx(bool vtx){
    if(!m_canChangeSettings) return false;
    m_doVtx = vtx;
    return true;
  }
  
  bool MinBiasTrackEfficiency::setUseTrigger(bool trig){
    if(!m_canChangeSettings) return false;
    m_doTrigger = trig;
    return true;
  }
  
}