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