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

Go to the documentation of this file.

#include "ForIA/AnalysisTools/MinBiasPhiGridCorrections.hh"
#include "ForIA/IDataLoader.hh"

#include "TH2F.h"

#include "boost/lexical_cast.hpp"

namespace ForIA{

  MinBiasPhiGridCorrections::MinBiasPhiGridCorrections(): m_segments(32), m_doInitialise(true),
  m_gotSegment(m_segments, false), m_gotGrid(false), m_randomGenerator(0){

  }

  MinBiasPhiGridCorrections::~MinBiasPhiGridCorrections(){
    if(m_randomGenerator != 0){
      gsl_rng_free(m_randomGenerator);
    }
  }
  
  void MinBiasPhiGridCorrections::initialise(){

    DataLoaderPtr dl = IDataLoader::create();
    dl->open("Root/MinBias/TruthReco.simple.quantiles.root");

    m_quantilePlot = dl->retrieve<TH2F>("TruthRecoCorrelations/quantile/self");
    /*
    for(int seg=0; seg != m_segments; ++seg){
      string path = "TruthRecoCorrelations/seg" + 
        boost::lexical_cast<string>(seg) + "/quantile/";

      TH2F *selfPlot = dl->retrieve<TH2F>(path + "self");
      m_selfPlots.push_back((TH2F*)selfPlot->Clone());

      TH2F *plusPlot = dl->retrieve<TH2F>(path + "diffPlus");
      m_plusPlots.push_back((TH2F*)plusPlot->Clone());

      TH2F *minusPlot = dl->retrieve<TH2F>(path + "diffMinus");
      m_minusPlots.push_back((TH2F*)minusPlot->Clone());

    }
    */
    dl->close();
    
    m_minPT = m_quantilePlot->GetXaxis()->GetBinUpEdge(1);
    
    std::cout<<"correction min pT = "<<m_minPT<<std::endl;
    
    // see http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html
    m_randomGenerator = gsl_rng_alloc (gsl_rng_taus2);
    
    gsl_rng_set(m_randomGenerator, 8912041);
    
    m_doInitialise = false;
    return;
  }

  bool MinBiasPhiGridCorrections::setGrid(const vector<double> &grid){
    
    if(m_doInitialise) initialise();
    
    for(vector<bool>::iterator seg = m_gotSegment.begin();
        seg != m_gotSegment.end(); ++seg){
      (*seg) = false;
    }

    m_gotGrid = false;
    
    m_input = grid;
    
    return true;
  }
  
  const vector<double> &MinBiasPhiGridCorrections::correctedGrid(){
        
    if(m_gotGrid) return m_corrected;
        
    
    vector<double> left;
    vector<double> right;
        
    m_corrected.clear();
    
    for(int seg=0; seg != m_segments; ++seg){
      
      double delta = getDelta(m_quantilePlot, m_input[seg]);
      if(delta < 500. && delta > 0.) delta = 0.;
      double result = m_input[seg] + delta;
      if(result < m_minPT) result = 0.;
      m_corrected.push_back(result);
    }
    /*
    for(int seg=0; seg != m_segments; ++seg){
      int segMinus=(seg==0)? m_segments-1: seg-1;
      double delta = getDelta(m_plusPlots[segMinus], m_input[seg]);
      double diff = m_input[seg] - delta;
      if(diff < 0.){
        diff = 0.;
        delta = m_input[seg];
      }
      left.push_back(delta);
      m_corrected.push_back(diff);
    }
    
    for(int seg=0; seg != m_segments; ++seg){
      int segPlus = seg + 1;
      if(segPlus == m_segments) segPlus = 0;
      double delta = getDelta(m_minusPlots[segPlus], m_input[seg]);
      double diff = m_corrected[seg] - delta;
      if(diff < 0.){
        diff = 0.;
        delta = m_corrected[seg];
      }
      right.push_back(delta);
      m_corrected[seg] = diff;
    }
    
    for(int seg=0; seg != m_segments; ++seg){
      double delta = getDelta(m_selfPlots[seg], m_input[seg]);
      int segPlus = seg + 1;
      if(segPlus == m_segments) segPlus = 0;

      int segMinus=(seg==0)? m_segments-1: seg-1;
     
      double cor = m_input[seg] + delta;// + left[segPlus] + right[segMinus];
      
      if(cor < 0.) cor=0.;
      
      m_corrected[seg] = cor;
      
    }
    */
    m_gotGrid = true;
    return m_corrected;
  }
  /*
  double MinBiasPhiGridCorrections::correctedSegment(int seg){
    
    if(m_gotSegment[seg]) return m_corrected[seg];
    
    int segPlus = seg + 1;
    if(segPlus == m_segments) segPlus = 0;
    int segMinus = seg - 1;
    if(segMinus == -1) segMinus = m_segments - 1;
    
    double deltaSelf  = getDelta(m_selfPlots[seg] , m_input[seg]);
    double deltaPlus  = getDelta(m_plusPlots[seg] , m_input[segPlus]);
    double deltaMinus = getDelta(m_minusPlots[seg], m_input[segMinus]);
    
    m_corrected[seg] = m_input[seg] + deltaSelf + deltaPlus + deltaMinus;
    
    return m_corrected[seg];
  }
*/
  double MinBiasPhiGridCorrections::getDelta(const TH2F *qPlot, double pT){
    
    int xbin = qPlot->GetXaxis()->FindBin(pT);
    
    double ran = gsl_rng_uniform(m_randomGenerator);
    
    int ybin = qPlot->GetYaxis()->FindBin(ran);
    
    if(ybin < 2) ybin=2;
    
    int binUp = qPlot->GetBin(xbin, ybin);
    int binDn = qPlot->GetBin(xbin, ybin - 1);
    
    double deltaUp = qPlot->GetBinContent(binUp);
    double deltaDn = qPlot->GetBinContent(binDn);
    
    double upEdge = qPlot->GetYaxis()->GetBinUpEdge(ybin);
    double dnEdge = qPlot->GetYaxis()->GetBinLowEdge(ybin);
    
    
    return (deltaUp - (upEdge - ran) * (deltaUp - deltaDn) / (upEdge - dnEdge));
  }
  
}