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 "THnSparse.h"
  33 #include "TAxis.h"
  34 #include "TMath.h"
  35 #include "TString.h"
  36
  37 #include "AliAODMCParticle.h"
  38 #include "AliAODpidUtil.h"
  39 #include "AliESDpid.h"
  40 #include "AliLog.h"
  41 #include "AliTPCPIDResponse.h"
  42 #include "AliTOFPIDResponse.h"
  43
  44 #include "AliHFEtools.h"
  45
  46 ClassImp(AliHFEtools)
  47
  48 AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
  49 Int_t AliHFEtools::fgLogLevel = 0;
  50
  51 //__________________________________________
  52 AliHFEtools::AliHFEtools():
  53   TObject()
  54 {
  55 }
  56
  57 //__________________________________________
  58 Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  59   //
  60   // Helper function for linearly binned array
  61   //
  62   Double_t *bins = new Double_t[nBins + 1];
  63   Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  64   bins[0] = ymin;
  65   for(Int_t ibin = 1; ibin <= nBins; ibin++)
  66     bins[ibin] = bins[ibin-1] + stepsize;
  67   return bins;
  68 }
  69
  70 //__________________________________________
  71 void AliHFEtools::FillLinearBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
  72   //
  73   // Helper function for linearly binned array
  74   //
  75   Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  76   bins[0] = ymin;
  77   for(Int_t ibin = 1; ibin <= nBins; ibin++)
  78     bins[ibin] = bins[ibin-1] + stepsize;
  79 }
  80
  81 //__________________________________________
  82 Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  83   //
  84   // Helper function for logartimically binned array
  85   //
  86   Double_t *bins = new Double_t[nBins+1];
  87   bins[0] = ymin;
  88   Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
  89   for(Int_t ibin = 1; ibin <= nBins; ibin++){
  90     bins[ibin] = factor * bins[ibin-1];
  91   }
  92   return bins;
  93 }
  94
  95 //__________________________________________
  96 void AliHFEtools::FillLogarithmicBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
  97   //
  98   // Helper function for logartimically binned array
  99   //
 100   bins[0] = ymin;
 101   Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
 102   for(Int_t ibin = 1; ibin <= nBins; ibin++)
 103     bins[ibin] = factor * bins[ibin-1];
 104 }
 105
 106 //_________________________________________
 107 Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
 108
 109   //
 110   // converts the axis (defined by the dimension) of THx or THnSparse
 111   // object to Log scale. Number of bins and bin min and bin max are preserved
 112   //
 113
 114
 115   if(!o){
 116     AliError("Input histogram is null pointer");
 117     return kFALSE;
 118   }
 119
 120   TAxis *axis = 0x0;
 121   if(o->InheritsFrom("TH1")){
 122     TH1 *h1 = dynamic_cast<TH1F*>(o);
 123     if(h1) axis = h1->GetXaxis();
 124   }
 125   else if(o->InheritsFrom("TH2")){
 126     TH2 *h2 = dynamic_cast<TH2F*>(o);
 127     if(h2 && 0 == dim){
 128       axis = h2->GetXaxis();
 129     }
 130     else if(h2 && 1 == dim){
 131       axis = h2->GetYaxis();
 132     }
 133      else{
 134        AliError("Only dim = 0 or 1 possible for TH2F");
 135      }
 136   }
 137   else if(o->InheritsFrom("THnSparse")){
 138     THnSparse *hs = dynamic_cast<THnSparse*>(o);
 139     if(hs) axis = hs->GetAxis(dim);
 140   }
 141   else{
 142     AliError("Type of input object not recognized, please check your code or update this finction");
 143     return kFALSE;
 144   }
 145   if(!axis){
 146     AliError(Form("Axis '%d' could not be identified in the object \n", dim));
 147     return kFALSE;
 148   }
 149
 150   Int_t bins = axis->GetNbins();
 151
 152   Double_t from = axis->GetXmin();
 153   if(from <= 0){
 154     AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
 155   }
 156   Double_t to = axis->GetXmax();
 157   TArrayD newBins(bins+1);
 158   newBins[0] = from;
 159   Double_t factor = TMath::Power(to/from, 1./bins);
 160   for(Int_t i=1; i<=bins; ++i){
 161     newBins[i] = factor * newBins[i-1];
 162   }
 163   axis->Set(bins, newBins.GetArray());
 164
 165   return kTRUE;
 166 }
 167
 168 //__________________________________________
 169 Float_t AliHFEtools::GetRapidity(const TParticle *part){
 170   //
 171   // return rapidity
 172   //
 173   Float_t rapidity;
 174   if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
 175   else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 176   return rapidity;
 177 }
 178
 179 //__________________________________________
 180 Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
 181   // return rapidity
 182
 183   Float_t rapidity;
 184   if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
 185   else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
 186   return rapidity;
 187 }
 188
 189 //__________________________________________
 190 AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
 191   //
 192   // Get the default PID as singleton instance
 193   //
 194   if(!fgDefaultPID){
 195
 196     if(isESD) fgDefaultPID = new AliESDpid;
 197     else fgDefaultPID = new AliAODpidUtil;
 198     Double_t tres = isMC ? 80. : 130.;
 199     fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
 200
 201     // TPC Bethe Bloch parameters
 202     Double_t alephParameters[5];
 203     if(isMC){
 204       // simulation
 205       alephParameters[0] = 2.15898e+00/50.;
 206       alephParameters[1] = 1.75295e+01;
 207       alephParameters[2] = 3.40030e-09;
 208       alephParameters[3] = 1.96178e+00;
 209       alephParameters[4] = 3.91720e+00;
 210     } else {
 211       alephParameters[0] = 0.0283086/0.97;
 212       //alephParameters[0] = 0.0283086;
 213       alephParameters[1] = 2.63394e+01;
 214       alephParameters[2] = 5.04114e-11;
 215       alephParameters[3] = 2.12543e+00;
 216       alephParameters[4] = 4.88663e+00;
 217     }
 218     fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
 219     fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
 220
 221   }
 222   if(fgLogLevel){
 223     printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
 224     printf("Error - Arrrrrrrrr...\n");
 225     printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
 226     printf("Error - TOF PID is adapted to Monte Carlo\n");
 227   }
 228   return fgDefaultPID;
 229 }
 230
 231
 232 //__________________________________________
 233 void AliHFEtools::DestroyDefaultPID(){
 234   //
 235   // Destroy default PID object if existing
 236   //
 237   if(fgDefaultPID) delete fgDefaultPID;
 238   fgDefaultPID = NULL;
 239 }
 240
 241 //__________________________________________
 242 Int_t AliHFEtools::GetPdg(const AliVParticle *track){
 243   //
 244   // Get MC PDG code (only MC particles supported)
 245   //
 246   Int_t pdg = 0;
 247   if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
 248     const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track);
 249     pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
 250   } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
 251     const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
 252     pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
 253   }
 254   return pdg;
 255 }
 256
 257 //__________________________________________
 258 Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
 259   //
 260   // Helper function to convert MC PID codes into AliPID codes
 261   //
 262   Int_t pid = -1;
 263   switch(TMath::Abs(pdg)){
 264     case 11: pid = AliPID::kElectron; break;
 265     case 13: pid = AliPID::kMuon; break;
 266     case 211: pid = AliPID::kPion; break;
 267     case 321: pid = AliPID::kKaon; break;
 268     case 2212: pid = AliPID::kProton; break;
 269     default: pid = -1; break;
 270   };
 271   return pid;
 272 }
 273
 274 //__________________________________________
 275 TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
 276 {
 277     // Creates a TH1D from TGraph g. The binwidth of the first bin has to
 278     // specified. The others bins are calculated automatically. Supports also Graphs
 279     // with non constant x steps. The axis of the Graph can be exchanged if
 280     // exchange=kTRUE (modifies the Graph).
 281
 282     TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
 283     if( result == 0) return result;
 284
 285     //------------------------------------------
 286     // setup the final hist
 287     Int_t nBinX = g->GetN();
 288     Double_t* err = g->GetEY();           // error y is still ok even if exchanged
 289     for(Int_t i = 0; i < nBinX; i ++){
 290         result->SetBinError(i + 1, err[i]);
 291     }
 292     if(exchange){
 293         AliHFEtools::ExchangeXYGraph(g);        // undo  what has been done in GraphToHist
 294         AliHFEtools::ExchangeXYGraphErrors(g);  // now exchange including errors
 295     }
 296
 297     return result;
 298 }
 299
 300 //__________________________________________
 301 Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
 302 {
 303     // exchanges x-values and y-values.
 304     if(g==0) return kFALSE;
 305     Int_t nbin=g->GetN();
 306     Double_t x,y;
 307     for(Int_t i = 0; i < nbin; i ++)
 308     {
 309         g->GetPoint(i,x,y);
 310         g->SetPoint(i,y,x);
 311     }
 312
 313     return kTRUE;
 314 }
 315
 316 //__________________________________________
 317 Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
 318 {
 319     // exchanges x-values and y-values and
 320     // corresponding errors.
 321     if(g==0) return kFALSE;
 322     Int_t nbin=g->GetN();
 323     Double_t x,y;
 324     Double_t ex,ey;
 325     for(Int_t i = 0; i < nbin; i ++)
 326     {
 327         g->GetPoint(i,x,y);
 328         ex = g->GetErrorX(i);
 329         ey = g->GetErrorY(i);
 330         g->SetPoint(i,y,x);
 331         g->SetPointError(i,ey,ex);
 332     }
 333
 334     return kTRUE;
 335
 336 }
 337
 338 //__________________________________________
 339 TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
 340 {
 341     // Creates a TH1D from TGraph g. The binwidth of the first bin has to
 342     // specified. The others bins are calculated automatically. Supports also Graphs
 343     // with non constant x steps. The axis of the Graph can be exchanged if
 344     // exchange=kTRUE (modifies the Graph).
 345
 346
 347     TH1D* result = 0;
 348     if(g == 0)              return result;
 349     if(firstBinWidth == -1) return result;
 350     TString myname="";
 351     myname = g->GetName();
 352     myname += "_graph";
 353     if(exchange) AliHFEtools::ExchangeXYGraph(g);
 354
 355     Int_t nBinX = g->GetN();
 356     Double_t* x = g->GetX();
 357     Double_t* y = g->GetY();
 358
 359     if(nBinX < 1) return result;
 360
 361     //------------------------------------------
 362     // create the Matrix for the equation system
 363     // and init the values
 364
 365     Int_t nDim = nBinX - 1;
 366     TMatrixD a(nDim,nDim);
 367     TMatrixD b(nDim,1);
 368
 369     Double_t* aA = a.GetMatrixArray();
 370     Double_t* aB = b.GetMatrixArray();
 371     memset(aA,0,nDim * nDim * sizeof(Double_t));
 372     memset(aB,0,nDim * sizeof(Double_t));
 373     //------------------------------------------
 374
 375     //------------------------------------------
 376     // setup equation system
 377     // width for 1st bin is given therefore
 378     // we shift bin parameter (column) by one to the left
 379     // to reduce the matrix size
 380
 381     Double_t* xAxis = new Double_t [nBinX + 1];
 382     Double_t* binW  = new Double_t [nBinX ];
 383     binW[0] = firstBinWidth;
 384
 385     aB[0] = x[1] - x[0]  - 0.5 * binW[0];
 386     aA[0] = 0.5;
 387
 388     for(Int_t col = 1; col < nDim ; col ++)
 389     {
 390         Int_t row = col;
 391         aB[col] = x[col + 1] - x[ col ];
 392         aA[row * nDim + col - 1 ] = 0.5;
 393         aA[row * nDim + col     ] = 0.5;
 394     }
 395     //------------------------------------------
 396
 397     //------------------------------------------
 398     // solve the equations
 399     a.Invert();
 400     TMatrixD c = a * b;
 401     //------------------------------------------
 402
 403     //------------------------------------------
 404     // calculate the bin boundaries
 405     xAxis[0] = x[0] - 0.5 * binW[0];
 406     memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
 407     for(Int_t col = 0; col < nBinX ; col ++) {
 408         xAxis[col + 1] = x[col] + 0.5 * binW[col];
 409     }
 410     //------------------------------------------
 411
 412     //------------------------------------------
 413     // setup the final hist
 414     result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
 415     for(Int_t i = 0; i < nBinX; i ++){
 416         result->SetBinContent(i + 1, y[i]);
 417     }
 418     result->SetMarkerColor(markercolor);
 419     result->SetMarkerStyle(markerstyle);
 420     result->SetMarkerSize(markersize);
 421     //------------------------------------------
 422
 423     delete [] xAxis;
 424     delete [] binW;
 425
 426
 427     return result;
 428 }
 429
 430 //__________________________________________
 431 void AliHFEtools::BinParameterisation(const TF1 &fun, const TArrayD &xbins, TArrayD &bincontent){
 432     //
 433     // Calculate binned version of a function defined as the integral of x*f(x) in
 434     // the integration range xmin,xmax, where xmin and xmax are the bin limits, divided
 435     // by the binwidth. The function is important in case of steeply falling functions
 436     //
 437     // Parameters
 438     //   fun:           the function to be binned
 439     //   xbins:         the bin limits
 440     //   bincontent:    the binned parameterisation
 441     //
 442     TString expression(Form("x*%s", fun.GetName()));
 443     Double_t xmin(0), xmax(0);
 444     fun.GetRange(xmin,xmax);
 445     // check range
 446     xmin = TMath::Min(xmin, xbins[0]);
 447     xmax = TMath::Max(xmax, xbins[xbins.GetSize()-1]);
 448     TF1 helper("helper",expression.Data(),xmin,xmax);   // make function x*f(x)
 449     if(bincontent.GetSize() != xbins.GetSize()-1)
 450         bincontent.Set(xbins.GetSize()-1); // Adapt array to number of bins
 451     //Caclulate Binned
 452     for(Int_t ib = 0; ib < xbins.GetSize()-1; ib++){
 453         xmin = xbins[ib];
 454         xmax = xbins[ib+1];
 455         bincontent[ib] = (helper.Integral(xmin, xmax))/(xmax - xmin);
 456     }
 457 }