]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEtools.cxx
Fix MC efficiencies more
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEtools.cxx
index c2a3d090312662fda433eecaf14c2778d23b9606..5bc40ce5b8af78f56c907154862cdb33c0bd024e 100644 (file)
 // Authors
 //   All authors of the HFE group
 //
+#include <TArrayD.h>
 #include <TMath.h>
 #include <TParticle.h>
+#include <TF1.h>
 #include "TH1.h"
+#include "TH1D.h"
 #include "TH2.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
+#include "TGraphAsymmErrors.h"
 #include "THnSparse.h"
 #include "TAxis.h"
+#include "TMath.h"
+#include "TString.h"
+#include "TFile.h"
+#include "TKey.h"
+#include "TROOT.h"
 
 #include "AliAODMCParticle.h"
 #include "AliAODpidUtil.h"
@@ -60,6 +71,17 @@ Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ym
   return bins;
 }
 
+//__________________________________________
+void AliHFEtools::FillLinearBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
+  //
+  // Helper function for linearly binned array
+  //
+  Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
+  bins[0] = ymin;
+  for(Int_t ibin = 1; ibin <= nBins; ibin++)
+    bins[ibin] = bins[ibin-1] + stepsize;
+}
+
 //__________________________________________
 Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
   //
@@ -74,6 +96,17 @@ Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double
   return bins;
 }
 
+//__________________________________________
+void AliHFEtools::FillLogarithmicBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
+  //
+  // Helper function for logartimically binned array
+  //
+  bins[0] = ymin;
+  Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
+  for(Int_t ibin = 1; ibin <= nBins; ibin++)
+    bins[ibin] = factor * bins[ibin-1];
+}
+
 //_________________________________________
 Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
 
@@ -125,14 +158,13 @@ Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
     AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
   }
   Double_t to = axis->GetXmax();
-  Double_t *newBins = new Double_t[bins+1];
+  TArrayD newBins(bins+1);
   newBins[0] = from;
   Double_t factor = TMath::Power(to/from, 1./bins);
   for(Int_t i=1; i<=bins; ++i){
     newBins[i] = factor * newBins[i-1];
   }
-  axis->Set(bins, newBins);
-  delete[] newBins;
+  axis->Set(bins, newBins.GetArray());
 
   return kTRUE;
 }
@@ -242,3 +274,293 @@ Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
   };
   return pid;
 }
+
+//__________________________________________
+TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
+{
+    // Creates a TH1D from TGraph g. The binwidth of the first bin has to
+    // specified. The others bins are calculated automatically. Supports also Graphs
+    // with non constant x steps. The axis of the Graph can be exchanged if
+    // exchange=kTRUE (modifies the Graph).
+
+    TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
+    if( result == 0) return result;
+
+    //------------------------------------------
+    // setup the final hist
+    Int_t nBinX = g->GetN();
+    Double_t* err = g->GetEY();           // error y is still ok even if exchanged
+    for(Int_t i = 0; i < nBinX; i ++){
+       result->SetBinError(i + 1, err[i]);
+    }
+    if(exchange){
+        AliHFEtools::ExchangeXYGraph(g);        // undo  what has been done in GraphToHist
+       AliHFEtools::ExchangeXYGraphErrors(g);  // now exchange including errors
+    }
+
+    return result;
+}
+
+//__________________________________________
+Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
+{
+    // exchanges x-values and y-values.
+    if(g==0) return kFALSE;
+    Int_t nbin=g->GetN();
+    Double_t x,y;
+    for(Int_t i = 0; i < nbin; i ++)
+    {
+        g->GetPoint(i,x,y);
+        g->SetPoint(i,y,x);
+    }
+
+    return kTRUE;
+}
+
+//__________________________________________
+Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
+{
+    // exchanges x-values and y-values and
+    // corresponding errors.
+    if(g==0) return kFALSE;
+    Int_t nbin=g->GetN();
+    Double_t x,y;
+    Double_t ex,ey;
+    for(Int_t i = 0; i < nbin; i ++)
+    {
+        g->GetPoint(i,x,y);
+        ex = g->GetErrorX(i);
+        ey = g->GetErrorY(i);
+        g->SetPoint(i,y,x);
+        g->SetPointError(i,ey,ex);
+    }
+
+    return kTRUE;
+
+}
+
+//__________________________________________
+TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
+{
+    // Creates a TH1D from TGraph g. The binwidth of the first bin has to
+    // specified. The others bins are calculated automatically. Supports also Graphs
+    // with non constant x steps. The axis of the Graph can be exchanged if
+    // exchange=kTRUE (modifies the Graph).
+
+
+    TH1D* result = 0;
+    if(g == 0)              return result;
+    if(firstBinWidth == -1) return result;
+    TString myname="";
+    myname = g->GetName();
+    myname += "_graph";
+    if(exchange) AliHFEtools::ExchangeXYGraph(g);
+
+    Int_t nBinX = g->GetN();
+    Double_t* x = g->GetX();
+    Double_t* y = g->GetY();
+
+    if(nBinX < 1) return result;
+
+    //------------------------------------------
+    // create the Matrix for the equation system
+    // and init the values
+
+    Int_t nDim = nBinX - 1;
+    TMatrixD a(nDim,nDim);
+    TMatrixD b(nDim,1);
+
+    Double_t* aA = a.GetMatrixArray();
+    Double_t* aB = b.GetMatrixArray();
+    memset(aA,0,nDim * nDim * sizeof(Double_t));
+    memset(aB,0,nDim * sizeof(Double_t));
+    //------------------------------------------
+
+    //------------------------------------------
+    // setup equation system
+    // width for 1st bin is given therefore
+    // we shift bin parameter (column) by one to the left
+    // to reduce the matrix size
+
+    Double_t* xAxis = new Double_t [nBinX + 1];
+    Double_t* binW  = new Double_t [nBinX ];
+    binW[0] = firstBinWidth;
+
+    aB[0] = x[1] - x[0]  - 0.5 * binW[0];
+    aA[0] = 0.5;
+
+    for(Int_t col = 1; col < nDim ; col ++)
+    {
+       Int_t row = col;
+       aB[col] = x[col + 1] - x[ col ];
+       aA[row * nDim + col - 1 ] = 0.5;
+       aA[row * nDim + col     ] = 0.5;
+    }
+    //------------------------------------------
+
+    //------------------------------------------
+    // solve the equations
+    a.Invert();
+    TMatrixD c = a * b;
+    //------------------------------------------
+
+    //------------------------------------------
+    // calculate the bin boundaries
+    xAxis[0] = x[0] - 0.5 * binW[0];
+    memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
+    for(Int_t col = 0; col < nBinX ; col ++) {
+       xAxis[col + 1] = x[col] + 0.5 * binW[col];
+    }
+    //------------------------------------------
+
+    //------------------------------------------
+    // setup the final hist
+    result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
+    for(Int_t i = 0; i < nBinX; i ++){
+       result->SetBinContent(i + 1, y[i]);
+    }
+    result->SetMarkerColor(markercolor);
+    result->SetMarkerStyle(markerstyle);
+    result->SetMarkerSize(markersize);
+    //------------------------------------------
+
+    delete [] xAxis;
+    delete [] binW;
+
+
+    return result;
+}
+
+//__________________________________________
+void AliHFEtools::BinParameterisation(const TF1 &fun, const TArrayD &xbins, TArrayD &bincontent){
+    //
+    // Calculate binned version of a function defined as the integral of x*f(x) in
+    // the integration range xmin,xmax, where xmin and xmax are the bin limits, divided
+    // by the binwidth. The function is important in case of steeply falling functions
+    //
+    // Parameters
+    //   fun:           the function to be binned
+    //   xbins:         the bin limits
+    //   bincontent:    the binned parameterisation
+    //
+    TString expression(Form("x*%s", fun.GetName()));
+    Double_t xmin(0), xmax(0);
+    fun.GetRange(xmin,xmax);
+    // check range
+    xmin = TMath::Min(xmin, xbins[0]);
+    xmax = TMath::Max(xmax, xbins[xbins.GetSize()-1]);
+    TF1 helper("helper",expression.Data(),xmin,xmax);   // make function x*f(x)
+    if(bincontent.GetSize() != xbins.GetSize()-1)
+        bincontent.Set(xbins.GetSize()-1); // Adapt array to number of bins
+    //Caclulate Binned
+    for(Int_t ib = 0; ib < xbins.GetSize()-1; ib++){
+        xmin = xbins[ib];
+        xmax = xbins[ib+1];
+        bincontent[ib] = (helper.Integral(xmin, xmax))/(xmax - xmin);
+    }
+}
+
+
+
+
+//_________________________________________________________________________
+//Function  AliHFEtools::GetHFEResultList() - opens file from argument and returns TList Object containing String "Results"
+//_________________________________________________________________________
+TList *AliHFEtools::GetHFEResultList(const TString str){
+
+    TFile *f = TFile::Open(str.Data());
+    if(!f || f->IsZombie()){
+        printf("Could not read file %s\n",str.Data()); 
+        return NULL ;
+    }
+    gROOT->cd();
+    TKey *k;
+    TIter next(f->GetListOfKeys());
+    while ((k = dynamic_cast<TKey *>(next()))){
+        TString s(k->GetName());
+        if(s.Contains("Results")) break;
+    }
+    if(!k){
+        printf("Output container not found\n");
+        f->Close(); delete f;
+        return NULL;
+    } 
+    TList *returnlist = dynamic_cast<TList *>(k->ReadObj());
+    f->Close(); delete f;
+    return returnlist;
+}
+
+
+//_________________________________________________________________________
+//Function  AliHFEtools::GetHFEQAList() - opens file from argument and returns TList Object containing String "QA"
+//_________________________________________________________________________
+TList *AliHFEtools::GetHFEQAList(const TString str){
+
+    TFile *f = TFile::Open(str.Data());
+    if(!f || f->IsZombie()){
+        printf("Could not read file %s\n",str.Data()); 
+        return NULL ;
+    }
+    gROOT->cd();
+    TKey *k;
+    TIter next(f->GetListOfKeys());
+    while ((k = dynamic_cast<TKey *>(next()))){
+        TString s(k->GetName());
+        if(s.Contains("QA")) break;
+    }
+    if(!k){
+        printf("Output container not found\n");
+        f->Close(); delete f;
+        return NULL;
+    } 
+    TList *returnlist = dynamic_cast<TList *>(k->ReadObj());
+    f->Close(); delete f;
+    return returnlist;
+}
+
+//__________________________________________
+void AliHFEtools::NormaliseBinWidth(TH1 *histo){
+  //
+  // Helper function to correct histograms for the bin width
+  //
+  Double_t binwidth(0.);
+  for(Int_t ipt = 1; ipt <= histo->GetNbinsX(); ipt++){
+    binwidth = histo->GetBinWidth(ipt);
+    histo->SetBinContent(ipt, histo->GetBinContent(ipt)/binwidth);
+    histo->SetBinError(ipt, histo->GetBinError(ipt)/binwidth);
+  }
+}
+
+//__________________________________________
+void AliHFEtools::NormaliseBinWdith(TGraphErrors *graph){
+  //
+  // Helper function to correct graphs with symmetric errors 
+  // for the bin width
+  //
+  Double_t binwidth(0.);
+  Double_t *ypoints = graph->GetY(),
+           *yerrors = graph->GetEY();
+  for(int ipt = 0; ipt < graph->GetN(); ipt++){
+    binwidth = 2*graph->GetEX()[ipt];
+    ypoints[ipt] /= binwidth;
+    yerrors[ipt] /= binwidth;
+  }
+}
+
+//__________________________________________
+void AliHFEtools::NormaliseBinWdithAsymm(TGraphAsymmErrors *graph){
+  //
+  // Helper function to correct graphs with asymmetric errors 
+  // for the bin width
+  //
+  Double_t binwidth(0.);
+  Double_t *ypoints = graph->GetY(),
+           *yerrorslow = graph->GetEYlow(),
+           *yerrorshigh = graph->GetEYhigh();
+  for(int ipt = 0; ipt < graph->GetN(); ipt++){
+    binwidth = graph->GetEXlow()[ipt] + graph->GetEXhigh()[ipt];
+    ypoints[ipt] /= binwidth;
+    yerrorslow[ipt] /= binwidth;
+    yerrorshigh[ipt] /= binwidth;
+  }
+}