PWGHF/hfe/AliHFEtools.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Toolkit containing various usefull things
  17 // Usable everywhere in the hfe software package
  18 // For more information see the cxx file
  19 //
  20 // Authors
  21 //   All authors of the HFE group
  22 //
  23 #include <TArrayD.h>
  24 #include <TMath.h>
  25 #include <TParticle.h>
  26 #include <TF1.h>
  27 #include "TH1.h"
  28 #include "TH1D.h"
  29 #include "TH2.h"
  30 #include "TGraph.h"
  31 #include "TGraphErrors.h"
  32 #include "TGraphAsymmErrors.h"
  33 #include "THnSparse.h"
  34 #include "TAxis.h"
  35 #include "TMath.h"
  36 #include "TString.h"
  37 #include "TFile.h"
  38 #include "TKey.h"
  39 #include "TROOT.h"
  40
  41 #include "AliAODMCParticle.h"
  42 #include "AliAODpidUtil.h"
  43 #include "AliESDpid.h"
  44 #include "AliLog.h"
  45 #include "AliTPCPIDResponse.h"
  46 #include "AliTOFPIDResponse.h"
  47
  48 #include "AliHFEtools.h"
  49
  50 ClassImp(AliHFEtools)
  51
  52 AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
  53 Int_t AliHFEtools::fgLogLevel = 0;
  54
  55 //__________________________________________
  56 AliHFEtools::AliHFEtools():
  57   TObject()
  58 {
  59 }
  60
  61 //__________________________________________
  62 Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  63   //
  64   // Helper function for linearly binned array
  65   //
  66   Double_t *bins = new Double_t[nBins + 1];
  67   Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  68   bins[0] = ymin;
  69   for(Int_t ibin = 1; ibin <= nBins; ibin++)
  70     bins[ibin] = bins[ibin-1] + stepsize;
  71   return bins;
  72 }
  73
  74 //__________________________________________
  75 void AliHFEtools::FillLinearBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
  76   //
  77   // Helper function for linearly binned array
  78   //
  79   Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  80   bins[0] = ymin;
  81   for(Int_t ibin = 1; ibin <= nBins; ibin++)
  82     bins[ibin] = bins[ibin-1] + stepsize;
  83 }
  84
  85 //__________________________________________
  86 Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  87   //
  88   // Helper function for logartimically binned array
  89   //
  90   Double_t *bins = new Double_t[nBins+1];
  91   bins[0] = ymin;
  92   Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
  93   for(Int_t ibin = 1; ibin <= nBins; ibin++){
  94     bins[ibin] = factor * bins[ibin-1];
  95   }
  96   return bins;
  97 }
  98
  99 //__________________________________________
 100 void AliHFEtools::FillLogarithmicBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
 101   //
 102   // Helper function for logartimically binned array
 103   //
 104   bins[0] = ymin;
 105   Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
 106   for(Int_t ibin = 1; ibin <= nBins; ibin++)
 107     bins[ibin] = factor * bins[ibin-1];
 108 }
 109
 110 //_________________________________________
 111 Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
 112
 113   //
 114   // converts the axis (defined by the dimension) of THx or THnSparse
 115   // object to Log scale. Number of bins and bin min and bin max are preserved
 116   //
 117
 118
 119   if(!o){
 120     AliError("Input histogram is null pointer");
 121     return kFALSE;
 122   }
 123
 124   TAxis *axis = 0x0;
 125   if(o->InheritsFrom("TH1")){
 126     TH1 *h1 = dynamic_cast<TH1F*>(o);
 127     if(h1) axis = h1->GetXaxis();
 128   }
 129   else if(o->InheritsFrom("TH2")){
 130     TH2 *h2 = dynamic_cast<TH2F*>(o);
 131     if(h2 && 0 == dim){
 132       axis = h2->GetXaxis();
 133     }
 134     else if(h2 && 1 == dim){
 135       axis = h2->GetYaxis();
 136     }
 137      else{
 138        AliError("Only dim = 0 or 1 possible for TH2F");
 139      }
 140   }
 141   else if(o->InheritsFrom("THnSparse")){
 142     THnSparse *hs = dynamic_cast<THnSparse*>(o);
 143     if(hs) axis = hs->GetAxis(dim);
 144   }
 145   else{
 146     AliError("Type of input object not recognized, please check your code or update this finction");
 147     return kFALSE;
 148   }
 149   if(!axis){
 150     AliError(Form("Axis '%d' could not be identified in the object \n", dim));
 151     return kFALSE;
 152   }
 153
 154   Int_t bins = axis->GetNbins();
 155
 156   Double_t from = axis->GetXmin();
 157   if(from <= 0){
 158     AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
 159   }
 160   Double_t to = axis->GetXmax();
 161   TArrayD newBins(bins+1);
 162   newBins[0] = from;
 163   Double_t factor = TMath::Power(to/from, 1./bins);
 164   for(Int_t i=1; i<=bins; ++i){
 165     newBins[i] = factor * newBins[i-1];
 166   }
 167   axis->Set(bins, newBins.GetArray());
 168
 169   return kTRUE;
 170 }
 171
 172 //__________________________________________
 173 Float_t AliHFEtools::GetRapidity(const TParticle *part){
 174   //
 175   // return rapidity
 176   //
 177   Float_t rapidity;
 178   if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
 179   else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 180   return rapidity;
 181 }
 182
 183 //__________________________________________
 184 Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
 185   // return rapidity
 186
 187   Float_t rapidity;
 188   if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
 189   else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
 190   return rapidity;
 191 }
 192
 193 //__________________________________________
 194 AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
 195   //
 196   // Get the default PID as singleton instance
 197   //
 198   if(!fgDefaultPID){
 199
 200     if(isESD) fgDefaultPID = new AliESDpid;
 201     else fgDefaultPID = new AliAODpidUtil;
 202     Double_t tres = isMC ? 80. : 130.;
 203     fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
 204
 205     // TPC Bethe Bloch parameters
 206     Double_t alephParameters[5];
 207     if(isMC){
 208       // simulation
 209       alephParameters[0] = 2.15898e+00/50.;
 210       alephParameters[1] = 1.75295e+01;
 211       alephParameters[2] = 3.40030e-09;
 212       alephParameters[3] = 1.96178e+00;
 213       alephParameters[4] = 3.91720e+00;
 214     } else {
 215       alephParameters[0] = 0.0283086/0.97;
 216       //alephParameters[0] = 0.0283086;
 217       alephParameters[1] = 2.63394e+01;
 218       alephParameters[2] = 5.04114e-11;
 219       alephParameters[3] = 2.12543e+00;
 220       alephParameters[4] = 4.88663e+00;
 221     }
 222     fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
 223     fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
 224
 225   }
 226   if(fgLogLevel){
 227     printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
 228     printf("Error - Arrrrrrrrr...\n");
 229     printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
 230     printf("Error - TOF PID is adapted to Monte Carlo\n");
 231   }
 232   return fgDefaultPID;
 233 }
 234
 235
 236 //__________________________________________
 237 void AliHFEtools::DestroyDefaultPID(){
 238   //
 239   // Destroy default PID object if existing
 240   //
 241   if(fgDefaultPID) delete fgDefaultPID;
 242   fgDefaultPID = NULL;
 243 }
 244
 245 //__________________________________________
 246 Int_t AliHFEtools::GetPdg(const AliVParticle *track){
 247   //
 248   // Get MC PDG code (only MC particles supported)
 249   //
 250   Int_t pdg = 0;
 251   if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
 252     const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track);
 253     pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
 254   } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
 255     const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
 256     pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
 257   }
 258   return pdg;
 259 }
 260
 261 //__________________________________________
 262 Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
 263   //
 264   // Helper function to convert MC PID codes into AliPID codes
 265   //
 266   Int_t pid = -1;
 267   switch(TMath::Abs(pdg)){
 268     case 11: pid = AliPID::kElectron; break;
 269     case 13: pid = AliPID::kMuon; break;
 270     case 211: pid = AliPID::kPion; break;
 271     case 321: pid = AliPID::kKaon; break;
 272     case 2212: pid = AliPID::kProton; break;
 273     default: pid = -1; break;
 274   };
 275   return pid;
 276 }
 277
 278 //__________________________________________
 279 TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
 280 {
 281     // Creates a TH1D from TGraph g. The binwidth of the first bin has to
 282     // specified. The others bins are calculated automatically. Supports also Graphs
 283     // with non constant x steps. The axis of the Graph can be exchanged if
 284     // exchange=kTRUE (modifies the Graph).
 285
 286     TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
 287     if( result == 0) return result;
 288
 289     //------------------------------------------
 290     // setup the final hist
 291     Int_t nBinX = g->GetN();
 292     Double_t* err = g->GetEY();           // error y is still ok even if exchanged
 293     for(Int_t i = 0; i < nBinX; i ++){
 294         result->SetBinError(i + 1, err[i]);
 295     }
 296     if(exchange){
 297         AliHFEtools::ExchangeXYGraph(g);        // undo  what has been done in GraphToHist
 298         AliHFEtools::ExchangeXYGraphErrors(g);  // now exchange including errors
 299     }
 300
 301     return result;
 302 }
 303
 304 //__________________________________________
 305 Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
 306 {
 307     // exchanges x-values and y-values.
 308     if(g==0) return kFALSE;
 309     Int_t nbin=g->GetN();
 310     Double_t x,y;
 311     for(Int_t i = 0; i < nbin; i ++)
 312     {
 313         g->GetPoint(i,x,y);
 314         g->SetPoint(i,y,x);
 315     }
 316
 317     return kTRUE;
 318 }
 319
 320 //__________________________________________
 321 Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
 322 {
 323     // exchanges x-values and y-values and
 324     // corresponding errors.
 325     if(g==0) return kFALSE;
 326     Int_t nbin=g->GetN();
 327     Double_t x,y;
 328     Double_t ex,ey;
 329     for(Int_t i = 0; i < nbin; i ++)
 330     {
 331         g->GetPoint(i,x,y);
 332         ex = g->GetErrorX(i);
 333         ey = g->GetErrorY(i);
 334         g->SetPoint(i,y,x);
 335         g->SetPointError(i,ey,ex);
 336     }
 337
 338     return kTRUE;
 339
 340 }
 341
 342 //__________________________________________
 343 TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
 344 {
 345     // Creates a TH1D from TGraph g. The binwidth of the first bin has to
 346     // specified. The others bins are calculated automatically. Supports also Graphs
 347     // with non constant x steps. The axis of the Graph can be exchanged if
 348     // exchange=kTRUE (modifies the Graph).
 349
 350
 351     TH1D* result = 0;
 352     if(g == 0)              return result;
 353     if(firstBinWidth == -1) return result;
 354     TString myname="";
 355     myname = g->GetName();
 356     myname += "_graph";
 357     if(exchange) AliHFEtools::ExchangeXYGraph(g);
 358
 359     Int_t nBinX = g->GetN();
 360     Double_t* x = g->GetX();
 361     Double_t* y = g->GetY();
 362
 363     if(nBinX < 1) return result;
 364
 365     //------------------------------------------
 366     // create the Matrix for the equation system
 367     // and init the values
 368
 369     Int_t nDim = nBinX - 1;
 370     TMatrixD a(nDim,nDim);
 371     TMatrixD b(nDim,1);
 372
 373     Double_t* aA = a.GetMatrixArray();
 374     Double_t* aB = b.GetMatrixArray();
 375     memset(aA,0,nDim * nDim * sizeof(Double_t));
 376     memset(aB,0,nDim * sizeof(Double_t));
 377     //------------------------------------------
 378
 379     //------------------------------------------
 380     // setup equation system
 381     // width for 1st bin is given therefore
 382     // we shift bin parameter (column) by one to the left
 383     // to reduce the matrix size
 384
 385     Double_t* xAxis = new Double_t [nBinX + 1];
 386     Double_t* binW  = new Double_t [nBinX ];
 387     binW[0] = firstBinWidth;
 388
 389     aB[0] = x[1] - x[0]  - 0.5 * binW[0];
 390     aA[0] = 0.5;
 391
 392     for(Int_t col = 1; col < nDim ; col ++)
 393     {
 394         Int_t row = col;
 395         aB[col] = x[col + 1] - x[ col ];
 396         aA[row * nDim + col - 1 ] = 0.5;
 397         aA[row * nDim + col     ] = 0.5;
 398     }
 399     //------------------------------------------
 400
 401     //------------------------------------------
 402     // solve the equations
 403     a.Invert();
 404     TMatrixD c = a * b;
 405     //------------------------------------------
 406
 407     //------------------------------------------
 408     // calculate the bin boundaries
 409     xAxis[0] = x[0] - 0.5 * binW[0];
 410     memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
 411     for(Int_t col = 0; col < nBinX ; col ++) {
 412         xAxis[col + 1] = x[col] + 0.5 * binW[col];
 413     }
 414     //------------------------------------------
 415
 416     //------------------------------------------
 417     // setup the final hist
 418     result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
 419     for(Int_t i = 0; i < nBinX; i ++){
 420         result->SetBinContent(i + 1, y[i]);
 421     }
 422     result->SetMarkerColor(markercolor);
 423     result->SetMarkerStyle(markerstyle);
 424     result->SetMarkerSize(markersize);
 425     //------------------------------------------
 426
 427     delete [] xAxis;
 428     delete [] binW;
 429
 430
 431     return result;
 432 }
 433
 434 //__________________________________________
 435 void AliHFEtools::BinParameterisation(const TF1 &fun, const TArrayD &xbins, TArrayD &bincontent){
 436     //
 437     // Calculate binned version of a function defined as the integral of x*f(x) in
 438     // the integration range xmin,xmax, where xmin and xmax are the bin limits, divided
 439     // by the binwidth. The function is important in case of steeply falling functions
 440     //
 441     // Parameters
 442     //   fun:           the function to be binned
 443     //   xbins:         the bin limits
 444     //   bincontent:    the binned parameterisation
 445     //
 446     TString expression(Form("x*%s", fun.GetName()));
 447     Double_t xmin(0), xmax(0);
 448     fun.GetRange(xmin,xmax);
 449     // check range
 450     xmin = TMath::Min(xmin, xbins[0]);
 451     xmax = TMath::Max(xmax, xbins[xbins.GetSize()-1]);
 452     TF1 helper("helper",expression.Data(),xmin,xmax);   // make function x*f(x)
 453     if(bincontent.GetSize() != xbins.GetSize()-1)
 454         bincontent.Set(xbins.GetSize()-1); // Adapt array to number of bins
 455     //Caclulate Binned
 456     for(Int_t ib = 0; ib < xbins.GetSize()-1; ib++){
 457         xmin = xbins[ib];
 458         xmax = xbins[ib+1];
 459         bincontent[ib] = (helper.Integral(xmin, xmax))/(xmax - xmin);
 460     }
 461 }
 462
 463
 464
 465
 466 //_________________________________________________________________________
 467 //Function  AliHFEtools::GetHFEResultList() - opens file from argument and returns TList Object containing String "Results"
 468 //_________________________________________________________________________
 469 TList *AliHFEtools::GetHFEResultList(const TString str){
 470
 471     TFile *f = TFile::Open(str.Data());
 472     if(!f || f->IsZombie()){
 473         printf("Could not read file %s\n",str.Data());
 474         return NULL ;
 475     }
 476     gROOT->cd();
 477     TKey *k;
 478     TIter next(f->GetListOfKeys());
 479     while ((k = dynamic_cast<TKey *>(next()))){
 480         TString s(k->GetName());
 481         if(s.Contains("Results")) break;
 482     }
 483     if(!k){
 484         printf("Output container not found\n");
 485         f->Close(); delete f;
 486         return NULL;
 487     }
 488     TList *returnlist = dynamic_cast<TList *>(k->ReadObj());
 489     f->Close(); delete f;
 490     return returnlist;
 491 }
 492
 493
 494 //_________________________________________________________________________
 495 //Function  AliHFEtools::GetHFEQAList() - opens file from argument and returns TList Object containing String "QA"
 496 //_________________________________________________________________________
 497 TList *AliHFEtools::GetHFEQAList(const TString str){
 498
 499     TFile *f = TFile::Open(str.Data());
 500     if(!f || f->IsZombie()){
 501         printf("Could not read file %s\n",str.Data());
 502         return NULL ;
 503     }
 504     gROOT->cd();
 505     TKey *k;
 506     TIter next(f->GetListOfKeys());
 507     while ((k = dynamic_cast<TKey *>(next()))){
 508         TString s(k->GetName());
 509         if(s.Contains("QA")) break;
 510     }
 511     if(!k){
 512         printf("Output container not found\n");
 513         f->Close(); delete f;
 514         return NULL;
 515     }
 516     TList *returnlist = dynamic_cast<TList *>(k->ReadObj());
 517     f->Close(); delete f;
 518     return returnlist;
 519 }
 520
 521 //__________________________________________
 522 void AliHFEtools::NormaliseBinWidth(TH1 *histo){
 523   //
 524   // Helper function to correct histograms for the bin width
 525   //
 526   Double_t binwidth(0.);
 527   for(Int_t ipt = 1; ipt <= histo->GetNbinsX(); ipt++){
 528     binwidth = histo->GetBinWidth(ipt);
 529     histo->SetBinContent(ipt, histo->GetBinContent(ipt)/binwidth);
 530     histo->SetBinError(ipt, histo->GetBinError(ipt)/binwidth);
 531   }
 532 }
 533
 534 //__________________________________________
 535 void AliHFEtools::NormaliseBinWdith(TGraphErrors *graph){
 536   //
 537   // Helper function to correct graphs with symmetric errors
 538   // for the bin width
 539   //
 540   Double_t binwidth(0.);
 541   Double_t *ypoints = graph->GetY(),
 542            *yerrors = graph->GetEY();
 543   for(int ipt = 0; ipt < graph->GetN(); ipt++){
 544     binwidth = 2*graph->GetEX()[ipt];
 545     ypoints[ipt] /= binwidth;
 546     yerrors[ipt] /= binwidth;
 547   }
 548 }
 549
 550 //__________________________________________
 551 void AliHFEtools::NormaliseBinWdithAsymm(TGraphAsymmErrors *graph){
 552   //
 553   // Helper function to correct graphs with asymmetric errors
 554   // for the bin width
 555   //
 556   Double_t binwidth(0.);
 557   Double_t *ypoints = graph->GetY(),
 558            *yerrorslow = graph->GetEYlow(),
 559            *yerrorshigh = graph->GetEYhigh();
 560   for(int ipt = 0; ipt < graph->GetN(); ipt++){
 561     binwidth = graph->GetEXlow()[ipt] + graph->GetEXhigh()[ipt];
 562     ypoints[ipt] /= binwidth;
 563     yerrorslow[ipt] /= binwidth;
 564     yerrorshigh[ipt] /= binwidth;
 565   }
 566 }