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