MadChi2Calc.cxx

Go to the documentation of this file.

#ifndef madchi2calc_cxx
#define madchi2calc_cxx
#include "TFile.h"
#include "TROOT.h"
#include "TMath.h"
#include "Mad/MadChi2Calc.h"

//******************************************
MadChi2Calc::MadChi2Calc(int ct) {
  
  calcType = ct;

  FluxMapHist = NULL;
  BinErr = 0.02;
  NormErr = 0.02;
  BkdXSecErr = 0.05;
  HeightErr = NULL;

}

//******************************************
MadChi2Calc::~MadChi2Calc() {
  
  TFile *fluxfile;
  if((fluxfile = (TFile*) gROOT->FindObject("FluxFile"))) {
    fluxfile->Close();
  }
  
}

//******************************************
void MadChi2Calc::SetFluxMap(TH2F *hist){
  
  if(hist==0) { //open from default file
    TFile *f = new TFile("CCHistos.root","READ");
    if(f) {
      f->SetName("FluxFile"); //set name so we can close file in destructor
      FluxMapHist = (TH2F*) f->Get("FluxMap");
    }
    else FluxMapHist = NULL;
    return;
  }
  FluxMapHist = hist;
  
}

//******************************************
void MadChi2Calc::SetSysErrors(double err1,double err2,
                               double *err3,double err4){
  this->SetNormalisationError(err1);
  this->SetBinToBinError(err2);
  this->SetBinHeightError(err3);
  this->SetBkdXSecError(err4);
}

//******************************************
double MadChi2Calc::GetChi2(TH1F *sighist,TH1F *mchist,
                            TH1F *bkdhist,double mcfactor) {

  if(!sighist||!mchist) return 0;
  int n = sighist->GetNbinsX();
  double Max = sighist->GetBinLowEdge(n+1);
  double Min = sighist->GetBinLowEdge(1);
  double *sig = new double[n];
  double *mc = new double[n];  
  double *bkd = new double[n];
  for(int i=0;i<n;i++){
    sig[i] = sighist->GetBinContent(i+1);
    mc[i]  = mchist->GetBinContent(i+1);
    if(bkdhist) bkd[i] = bkdhist->GetBinContent(i+1);
    else bkd[i] = 0;
  }

  double chi2 = SumChi2(n,sig,mc,bkd,mcfactor,Min,Max);

  delete [] sig;
  delete [] mc;
  delete [] bkd;

  return chi2; 

}

//*****************************************
double MadChi2Calc::GetChi2(TH2F *sighist,TH2F *mchist,
                            TH2F *bkdhist,double mcfactor) {

  if(!sighist||!mchist) return 0;
  int nx = sighist->GetNbinsX();
  int ny = sighist->GetNbinsY();
  double *sig = new double[nx*ny];
  double *mc = new double[nx*ny];  
  double *bkd = new double[nx*ny];
  for(int i=0;i<nx;i++){
    for(int j=0;j<ny;j++){
      sig[i*ny+j] = sighist->GetBinContent(i+1,j+1);
      mc[i*ny+j]  = mchist->GetBinContent(i+1,j+1);
      if(bkdhist) bkd[i*ny+j] = bkdhist->GetBinContent(i+1,j+1);
      else bkd[i*ny+j] = 0;
    }
  }

  double chi2 = SumChi2(nx*ny,sig,mc,bkd,mcfactor,0,1);

  delete [] sig;
  delete [] mc;
  delete [] bkd;

  return chi2; 

}

//******************************************
double MadChi2Calc::GetChi2(int n, double *sig,double *mc,
                            double *bkd,double mcfactor,double Min,double Max){

  return SumChi2(n,sig,mc,bkd,mcfactor,Min,Max);

}

//******************************************
double MadChi2Calc::SumChi2(int n,double *sig,double *mc,
                            double *bkd,double mcfactor,double Min,double Max){
  
  double chi2 = 0;
  switch(calcType) {
  case 1: 
    for(int i=0;i<n;i++){
      chi2 += BasicChi2(sig[i],mc[i]/mcfactor,bkd[i]/mcfactor);
    }
    break;
  case 2: 
    chi2 = WithSysErrChi2(n,sig,mc,bkd,mcfactor,Min,Max);
    break;
  case 3:
    chi2 = WithSysErrLoopsChi2(n,sig,mc,bkd,mcfactor);
    break;
  default: 
    for(int i=0;i<n;i++){
      chi2 += PoisDistChi2(sig[i],mc[i]/mcfactor,bkd[i]/mcfactor);
    }
    break;
  }

  return chi2;

}

//******************************************
double MadChi2Calc::BasicChi2(double sig,double mc,double bkd){
  
  return (mc+bkd-sig)*(mc+bkd-sig)/sig;

}

//******************************************
double MadChi2Calc::PoisDistChi2(double sig,double mc,double bkd){

  if(sig==0||mc+bkd==0) return 2.*(mc+bkd-sig);
  return 2.*(mc+bkd-sig)+2.*sig*TMath::Log(sig/(mc+bkd));

}

//******************************************
double MadChi2Calc::WithSysErrChi2(int n,double *sig,double *mc,
                                   double *bkd,double mcfactor,
                                   double Min,double Max){

  float fw = (Max-Min)/float(n);
  int nSquare = int(float(n)*float(n));
  float *ErrArray = new float[nSquare];
  float *InvErrArray = new float[nSquare];

  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){   
      
      float DLM = 0.;
      if(i==j) DLM=1.;
      
      float FSum=0.;
      if(FluxMapHist){
        for(int k=0;k<n;k++) {
          int binx_k = FluxMapHist->GetXaxis()->FindBin(float(k+0.5)*fw);
          int binx_j = FluxMapHist->GetXaxis()->FindBin(float(j+0.5)*fw);
          int biny   = FluxMapHist->GetYaxis()->FindBin(float(i+0.5)*fw);
          FSum += FluxMapHist->GetBinContent(binx_k,biny)
            *FluxMapHist->GetBinContent(binx_j,biny);
        }
      }
      
      ErrArray[(j*n)+i] = DLM*(mc[i]+bkd[i])/mcfactor
        + NormErr*NormErr*(mc[i]+bkd[i])*(mc[j]+bkd[j])/(mcfactor*mcfactor)
        + BinErr*BinErr*FSum*(mc[i]+bkd[i])*(mc[j]+bkd[j])/(mcfactor*mcfactor);
      
      if(HeightErr) ErrArray[(j*n)+i] += HeightErr[j]*HeightErr[j]
                      *(mc[i]+bkd[i])*(mc[j]+bkd[j])/(mcfactor*mcfactor);
      else ErrArray[(j*n)+i] += 0.02*0.02
             *(mc[i]+bkd[i])*(mc[j]+bkd[j])/(mcfactor*mcfactor);
    }
  }
  
  TMatrix matrix(n,n);
  matrix.SetMatrixArray(ErrArray);
  matrix.Invert(); //Error Matrix inversion
  matrix.GetMatrix2Array(InvErrArray);
  
  double ChiSq = 0;
  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++) {
      ChiSq+=(sig[i]-mc[i])*InvErrArray[(j*n)+i]*(sig[j]-mc[j]);
    }
  }
  return ChiSq;

}

//******************************************
double MadChi2Calc::WithSysErrLoopsChi2(int n,double *sig,double *mc,
                                        double *bkd,double mcfactor){
  
  int nNormBins = 40;
  double fracNormErr = 0.1; //up to 10% error calculated
  double normLoopStep = 2.*fracNormErr/double(nNormBins);
  
  int nBkdXSecBins = 40;
  double fracBkdXSecErr = 0.2; //up to 20% error calculated
  double bkdXSecLoopStep = 2.*fracBkdXSecErr/double(nBkdXSecBins);
  
  double ChiSqMin = 99999;
  
  //overall normalisation uncertainty:
  for(int i=0;i<nNormBins;i++){
    double normFactor = 1. - fracNormErr + i*normLoopStep;
    //Bkgd cross-section uncertainty
    for(int j=0;j<nBkdXSecBins;j++){
      double bkdXSecFactor = 1. - fracBkdXSecErr + j*bkdXSecLoopStep;
      //calculate sys errors
      double ChiSq = 0;
      for(int k=0;k<n;k++){
        ChiSq += PoisDistChi2(sig[k],normFactor*mc[k]/mcfactor,
                              normFactor*bkdXSecFactor*bkd[k]/mcfactor);
      }
      // Now need to add penalty terms to ChiSq 
      // depending on values of systematic shifts
      if(NormErr>0) ChiSq += TMath::Power((1.-normFactor)/NormErr,2);
      if(BkdXSecErr>0) ChiSq += TMath::Power((1.-bkdXSecFactor)/BkdXSecErr,2);
      if(ChiSq<ChiSqMin) ChiSqMin = ChiSq;
    }
  }
  return ChiSqMin;
}
#endif // #ifdef madchi2calc_cxx

---