PlotUtils.C

Go to the documentation of this file.

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <TH1D.h>
#include <TF1.h>
#include <TMath.h>
#include <TGraphErrors.h>

#endif

//-----------------------------------------------------------------------------
static Double_t FunctionCrystalBall(Double_t *x, Double_t *par)
{
  Double_t N = par[0];
  Double_t width = par[1];
  
  Double_t mean = par[2];
  Double_t sigma = par[3];
  Double_t alpha = par[4];
  Double_t n = par[5];
  
  Double_t A = pow(n/fabs(alpha),n)*exp(-pow(alpha,2)/2);
  Double_t B = n/fabs(alpha) - fabs(alpha);
  
  if ((x[0]-mean)/sigma>-alpha)
    return N*width*TMath::Gaus(x[0],mean,sigma,1);
  else
    return N/(sqrt(TMath::TwoPi())*sigma)*width*A*pow(B-(x[0]-mean)/sigma,-n);
}

//-----------------------------------------------------------------------------
static Double_t FunctionGaussPol3(Double_t *x, Double_t *par)
{
  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
                                      (2*par[2]*par[2]) );
  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0] + par[6]*x[0]*x[0]*x[0];
  return gaus+back;
}

//-----------------------------------------------------------------------------
static Double_t FunctionGaussPol2(Double_t *x, Double_t *par)
{
  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
                                      (2*par[2]*par[2]) );
  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
  return gaus+back;
}

//-----------------------------------------------------------------------------
static Double_t FunctionGaussPol1(Double_t *x, Double_t *par)
{
  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
                                      (2*par[2]*par[2]) );
  Double_t back = par[3] + par[4]*x[0];
  return gaus+back;
}

//-----------------------------------------------------------------------------
static Double_t FunctionGaussPol0(Double_t *x, Double_t *par)
{
  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
                                      (2*par[2]*par[2]) );
  Double_t back = 0;//par[3];
  return gaus+back;
}

//-----------------------------------------------------------------------------
static Double_t FunctionTruncatedPolEta(Double_t *x, Double_t *par)
{
  if ( x[0] > 0.425 && x[0] < 0.65 ) 
  {
    TF1::RejectPoint();
    return 0;
  }
  
  return par[0] + par[1]*x[0];
  // + x[0]*x[0]*par[2]  + x[0]*x[0]*x[0]*par[3];//+ x[0]*x[0]*x[0]*x[0]*par[4];
}

//-----------------------------------------------------------------------------
static Double_t FunctionTruncatedPolPi0(Double_t *x, Double_t *par)
{
  //if ( x[0] > 0.07 && x[0] < 0.2 )
  if ( x[0] < 0.25 )
  {
    TF1::RejectPoint();
    return 0;
  }
  
  return par[0] + par[1]*x[0];
  // + x[0]*x[0]*par[2] + x[0]*x[0]*x[0]*par[3] ;//+ x[0]*x[0]*x[0]*x[0]*par[4];
}

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
static Double_t GetFractionError
(Double_t num   , Double_t den, 
 Double_t numErr, Double_t denErr)
{
  if ( num == 0 || den == 0 ) return 0.;
  
  //printf("\t num %e den %e numErr %e denErr %e\n",num,den,numErr,denErr);
  
  return num/den * TMath::Sqrt( ( numErr * numErr ) / ( num * num ) + 
                                ( denErr * denErr ) / ( den * den )   );
}

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
static TGraphErrors * DivideGraphs
(TGraphErrors* gNum, TGraphErrors *gDen)
{
  if ( !gDen || !gNum ) 
  {
    printf("Graph num %p or graph den %p not available\n",gNum,gDen);
    return 0x0;
  }
  const Int_t nBins = gNum->GetN();
  if ( nBins != gDen->GetN() )
  {
    printf("Cannot divide %s with %d bins and %s with %d bins!\n",
           gNum->GetName(),nBins,gDen->GetName(),gDen->GetN());
    return 0x0;
  }
  
  Double_t ratio   [nBins];
  Double_t ratioErr[nBins]; 
  Double_t x       [nBins];
  Double_t xErr    [nBins];
  for (Int_t ibin = 0; ibin < nBins; ibin++) 
  {
    Double_t num    =  gNum->GetY ()[ibin];
    Double_t den    =  gDen->GetY ()[ibin];
    Double_t numErr =  gNum->GetEY()[ibin];
    Double_t denErr =  gDen->GetEY()[ibin];
    
    x   [ibin]      =  gNum->GetX ()[ibin];
    xErr[ibin]      =  gNum->GetEX()[ibin];
    
    if ( num == 0 || den == 0 ) 
    {
      ratio   [ibin] = 0; 
      ratioErr[ibin] = 0; 
      continue;
    }
    
    ratio   [ibin] = num / den ;
    ratioErr[ibin] =  GetFractionError(num,den,numErr,denErr);  
    
    //    printf("bin %d, x %f (%f) num %f (%f), den %f (%f), ratio %f (%f) \n",
    //           ibin,x[ibin],xErr[ibin],num,numErr,den,denErr,ratio[ibin],ratioErr[ibin]);
  } // do the ratio to sum
  
  return new TGraphErrors(nBins,x,ratio,xErr,ratioErr);
}


//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
static void GetRangeIntegralAndError
(TH1D* h, Int_t binMin, Int_t binMax, 
 Double_t & integral, Double_t & integralErr )
{
  integral    = 0;
  integralErr = 0;
  for(Int_t ibin = binMin; ibin <= binMax; ibin++)
  {
    //if ( h->GetBinContent(ibin) == 0 ) continue ; 
    
    integral    += h->GetBinContent(ibin);
    integralErr += h->GetBinError  (ibin) * h->GetBinError(ibin);
  }
//  printf("\t bin range [%d,%d], n %d, sum %2.2e err %2.2e\n",
//         binMin,binMax,n,integral,integralErr);
  if ( integralErr > 0 ) integralErr = TMath::Sqrt(integralErr);
}

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
static void GetRangeIntegralAndError
(TH1D* h, Float_t minX, Float_t maxX, 
 Double_t & integral, Double_t & integralErr )
{
  Int_t binMin = h->FindBin(minX);
  Int_t binMax = h->FindBin(maxX)-1;
  return GetRangeIntegralAndError(h, binMin, binMax, 
                                  integral, integralErr );
}

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
void RemovePointsOutOfRangeOrLargeErrorFromGraph
(TGraphErrors* graph, 
 Double_t min, Double_t max,
 Float_t errFraction = 0.5,
 Float_t value = -1, Float_t valueErr = 0 )
{
  for(Int_t ibin = 0; ibin < graph->GetN(); ibin++)
  {
    Float_t x  = graph->GetX()[ibin];
    Float_t y  = graph->GetY()[ibin];
    Float_t yE = graph->GetEY()[ibin];
    if(x < min || x > max )
    {
      graph->SetPoint     (ibin,x, value   );
      graph->SetPointError(ibin,x, valueErr);
    }
    else 
      if(yE > errFraction*y) // Remove points with more than 50% error
      {
        //printf("ibin %d, y %f, yE %f\n",ibin,y,yE);
        graph->SetPoint     (ibin,x,value   );
        graph->SetPointError(ibin,x,valueErr);
      }
  }
}

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
static void ScaleBinBySize(TH1D* h)
{
  for(Int_t ibin = 1; ibin < h->GetNbinsX();ibin++)
  {
    Double_t width   = h->GetBinWidth(ibin);
    Double_t content = h->GetBinContent(ibin);
    //printf("bin %d, width %f, content %e\n",ibin,width,content);
    h->SetBinContent(ibin,content/width);
  }
}