PWG3/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 "TH2.h"
  27 #include "THnSparse.h"
  28 #include "TAxis.h"
  29
  30 #include "AliAODMCParticle.h"
  31 #include "AliAODpidUtil.h"
  32 #include "AliESDpid.h"
  33 #include "AliLog.h"
  34 #include "AliTPCPIDResponse.h"
  35 #include "AliTOFPIDResponse.h"
  36
  37 #include "AliHFEtools.h"
  38
  39 ClassImp(AliHFEtools)
  40
  41 AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
  42 Int_t AliHFEtools::fgLogLevel = 0;
  43
  44 //__________________________________________
  45 AliHFEtools::AliHFEtools():
  46   TObject()
  47 {
  48 }
  49
  50 //__________________________________________
  51 Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  52   //
  53   // Helper function for linearly binned array
  54   //
  55   Double_t *bins = new Double_t[nBins + 1];
  56   Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  57   bins[0] = ymin;
  58   for(Int_t ibin = 1; ibin <= nBins; ibin++)
  59     bins[ibin] = bins[ibin-1] + stepsize;
  60   return bins;
  61 }
  62
  63 //__________________________________________
  64 Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  65   //
  66   // Helper function for logartimically binned array
  67   //
  68   Double_t *bins = new Double_t[nBins+1];
  69   bins[0] = ymin;
  70   Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
  71   for(Int_t ibin = 1; ibin <= nBins; ibin++){
  72     bins[ibin] = factor * bins[ibin-1];
  73   }
  74   return bins;
  75 }
  76
  77 //_________________________________________
  78 Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
  79
  80   //
  81   // converts the axis (defined by the dimension) of THx or THnSparse
  82   // object to Log scale. Number of bins and bin min and bin max are preserved
  83   //
  84
  85
  86   if(!o){
  87     AliError("Input histogram is null pointer");
  88     return kFALSE;
  89   }
  90
  91   TAxis *axis = 0x0;
  92   if(o->InheritsFrom("TH1")){
  93     TH1 *h1 = dynamic_cast<TH1F*>(o);
  94     if(h1) axis = h1->GetXaxis();
  95   }
  96   else if(o->InheritsFrom("TH2")){
  97     TH2 *h2 = dynamic_cast<TH2F*>(o);
  98     if(h2 && 0 == dim){
  99       axis = h2->GetXaxis();
 100     }
 101     else if(h2 && 1 == dim){
 102       axis = h2->GetYaxis();
 103     }
 104      else{
 105        AliError("Only dim = 0 or 1 possible for TH2F");
 106      }
 107   }
 108   else if(o->InheritsFrom("THnSparse")){
 109     THnSparse *hs = dynamic_cast<THnSparse*>(o);
 110     if(hs) axis = hs->GetAxis(dim);
 111   }
 112   else{
 113     AliError("Type of input object not recognized, please check your code or update this finction");
 114     return kFALSE;
 115   }
 116   if(!axis){
 117     AliError(Form("Axis '%d' could not be identified in the object \n", dim));
 118     return kFALSE;
 119   }
 120
 121   Int_t bins = axis->GetNbins();
 122
 123   Double_t from = axis->GetXmin();
 124   if(from <= 0){
 125     AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
 126   }
 127   Double_t to = axis->GetXmax();
 128   Double_t *newBins = new Double_t[bins+1];
 129   newBins[0] = from;
 130   Double_t factor = TMath::Power(to/from, 1./bins);
 131   for(Int_t i=1; i<=bins; ++i){
 132     newBins[i] = factor * newBins[i-1];
 133   }
 134   axis->Set(bins, newBins);
 135   delete[] newBins;
 136
 137   return kTRUE;
 138 }
 139
 140 //__________________________________________
 141 Float_t AliHFEtools::GetRapidity(const TParticle *part){
 142   //
 143   // return rapidity
 144   //
 145   Float_t rapidity;
 146   if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
 147   else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 148   return rapidity;
 149 }
 150
 151 //__________________________________________
 152 Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
 153   // return rapidity
 154
 155   Float_t rapidity;
 156   if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
 157   else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
 158   return rapidity;
 159 }
 160
 161 //__________________________________________
 162 AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
 163   //
 164   // Get the default PID as singleton instance
 165   //
 166   if(!fgDefaultPID){
 167
 168     if(isESD) fgDefaultPID = new AliESDpid;
 169     else fgDefaultPID = new AliAODpidUtil;
 170     Double_t tres = isMC ? 80. : 130.;
 171     fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
 172
 173     // TPC Bethe Bloch parameters
 174     Double_t alephParameters[5];
 175     if(isMC){
 176       // simulation
 177       alephParameters[0] = 2.15898e+00/50.;
 178       alephParameters[1] = 1.75295e+01;
 179       alephParameters[2] = 3.40030e-09;
 180       alephParameters[3] = 1.96178e+00;
 181       alephParameters[4] = 3.91720e+00;
 182     } else {
 183       alephParameters[0] = 0.0283086/0.97;
 184       //alephParameters[0] = 0.0283086;
 185       alephParameters[1] = 2.63394e+01;
 186       alephParameters[2] = 5.04114e-11;
 187       alephParameters[3] = 2.12543e+00;
 188       alephParameters[4] = 4.88663e+00;
 189     }
 190     fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
 191     fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
 192
 193   }
 194   if(fgLogLevel){
 195     printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
 196     printf("Error - Arrrrrrrrr...\n");
 197     printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
 198     printf("Error - TOF PID is adapted to Monte Carlo\n");
 199   }
 200   return fgDefaultPID;
 201 }
 202
 203
 204 //__________________________________________
 205 void AliHFEtools::DestroyDefaultPID(){
 206   //
 207   // Destroy default PID object if existing
 208   //
 209   if(fgDefaultPID) delete fgDefaultPID;
 210   fgDefaultPID = NULL;
 211 }
 212
 213 //__________________________________________
 214 Int_t AliHFEtools::GetPdg(const AliVParticle *track){
 215   //
 216   // Get MC PDG code (only MC particles supported)
 217   //
 218   Int_t pdg = 0;
 219   if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
 220     const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track);
 221     pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
 222   } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
 223     const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
 224     pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
 225   }
 226   return pdg;
 227 }
 228
 229 //__________________________________________
 230 Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
 231   //
 232   // Helper function to convert MC PID codes into AliPID codes
 233   //
 234   Int_t pid = -1;
 235   switch(TMath::Abs(pdg)){
 236     case 11: pid = AliPID::kElectron; break;
 237     case 13: pid = AliPID::kMuon; break;
 238     case 211: pid = AliPID::kPion; break;
 239     case 321: pid = AliPID::kKaon; break;
 240     case 2212: pid = AliPID::kProton; break;
 241     default: pid = -1; break;
 242   };
 243   return pid;
 244 }