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