1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Toolkit containing various usefull things
17 // Usable everywhere in the hfe software package
18 // For more information see the cxx file
21 // All authors of the HFE group
24 #include <TParticle.h>
27 #include "THnSparse.h"
30 #include "AliAODMCParticle.h"
31 #include "AliAODpidUtil.h"
32 #include "AliESDpid.h"
34 #include "AliTPCPIDResponse.h"
35 #include "AliTOFPIDResponse.h"
37 #include "AliHFEtools.h"
41 AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
42 Int_t AliHFEtools::fgLogLevel = 0;
44 //__________________________________________
45 AliHFEtools::AliHFEtools():
50 //__________________________________________
51 Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
53 // Helper function for linearly binned array
55 Double_t *bins = new Double_t[nBins + 1];
56 Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
58 for(Int_t ibin = 1; ibin <= nBins; ibin++)
59 bins[ibin] = bins[ibin-1] + stepsize;
63 //__________________________________________
64 Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
66 // Helper function for logartimically binned array
68 Double_t *bins = new Double_t[nBins+1];
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];
77 //_________________________________________
78 Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
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
87 AliError("Input histogram is null pointer");
92 if(o->InheritsFrom("TH1")){
93 TH1 *h1 = dynamic_cast<TH1F*>(o);
94 if(h1) axis = h1->GetXaxis();
96 else if(o->InheritsFrom("TH2")){
97 TH2 *h2 = dynamic_cast<TH2F*>(o);
99 axis = h2->GetXaxis();
101 else if(h2 && 1 == dim){
102 axis = h2->GetYaxis();
105 AliError("Only dim = 0 or 1 possible for TH2F");
108 else if(o->InheritsFrom("THnSparse")){
109 THnSparse *hs = dynamic_cast<THnSparse*>(o);
110 if(hs) axis = hs->GetAxis(dim);
113 AliError("Type of input object not recognized, please check your code or update this finction");
117 AliError(Form("Axis '%d' could not be identified in the object \n", dim));
121 Int_t bins = axis->GetNbins();
123 Double_t from = axis->GetXmin();
125 AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
127 Double_t to = axis->GetXmax();
128 Double_t *newBins = new Double_t[bins+1];
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];
134 axis->Set(bins, newBins);
140 //__________________________________________
141 Float_t AliHFEtools::GetRapidity(const TParticle *part){
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())));
151 //__________________________________________
152 Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
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())));
161 //__________________________________________
162 AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
164 // Get the default PID as singleton instance
168 if(isESD) fgDefaultPID = new AliESDpid;
169 else fgDefaultPID = new AliAODpidUtil;
170 Double_t tres = isMC ? 80. : 130.;
171 fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
173 // TPC Bethe Bloch parameters
174 Double_t alephParameters[5];
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;
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;
190 fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
191 fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
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");
204 //__________________________________________
205 void AliHFEtools::DestroyDefaultPID(){
207 // Destroy default PID object if existing
209 if(fgDefaultPID) delete fgDefaultPID;
213 //__________________________________________
214 Int_t AliHFEtools::GetPdg(const AliVParticle *track){
216 // Get MC PDG code (only MC particles supported)
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;
229 //__________________________________________
230 Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
232 // Helper function to convert MC PID codes into AliPID codes
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;