Cleanup the code. Fix memory leak. Now inherit from AliAnalysisTaskSE (Antoine, Phili...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEtools.cxx
CommitLineData
70da6c5a 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**************************************************************************/
27de2dfb 15
16/* $Id$ */
17
70da6c5a 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"
3a72645a 34#include "AliAODpidUtil.h"
faee3b18 35#include "AliESDpid.h"
70da6c5a 36#include "AliLog.h"
3a72645a 37#include "AliTPCPIDResponse.h"
faee3b18 38#include "AliTOFPIDResponse.h"
70da6c5a 39
40#include "AliHFEtools.h"
41
42ClassImp(AliHFEtools)
43
faee3b18 44AliESDpid *AliHFEtools::fgDefaultPID = NULL;
3a72645a 45AliAODpidUtil *AliHFEtools::fgDefaultPIDaod = NULL;
69ac0e6f 46Int_t AliHFEtools::fgLogLevel = 0;
faee3b18 47
70da6c5a 48//__________________________________________
49AliHFEtools::AliHFEtools():
50 TObject()
51{
52}
53
54//__________________________________________
55Double_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//__________________________________________
68Double_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//_________________________________________
82Bool_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")){
bf892a6a 97 TH1 *h1 = dynamic_cast<TH1F*>(o);
98 if(h1) axis = h1->GetXaxis();
70da6c5a 99 }
100 else if(o->InheritsFrom("TH2")){
bf892a6a 101 TH2 *h2 = dynamic_cast<TH2F*>(o);
102 if(h2 && 0 == dim){
103 axis = h2->GetXaxis();
70da6c5a 104 }
bf892a6a 105 else if(h2 && 1 == dim){
106 axis = h2->GetYaxis();
70da6c5a 107 }
108 else{
109 AliError("Only dim = 0 or 1 possible for TH2F");
110 }
111 }
112 else if(o->InheritsFrom("THnSparse")){
bf892a6a 113 THnSparse *hs = dynamic_cast<THnSparse*>(o);
114 if(hs) axis = hs->GetAxis(dim);
70da6c5a 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);
bf892a6a 139 delete[] newBins;
70da6c5a 140
141 return kTRUE;
142}
143
144//__________________________________________
3a72645a 145Float_t AliHFEtools::GetRapidity(const TParticle *part){
70da6c5a 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//__________________________________________
3a72645a 156Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
70da6c5a 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
faee3b18 165//__________________________________________
166AliESDpid* 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]);
67fe7bd0 193 fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
faee3b18 194
195 }
196 if(fgLogLevel){
67fe7bd0 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");
faee3b18 201 }
202 return fgDefaultPID;
203}
204
205//__________________________________________
3a72645a 206AliAODpidUtil* 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//__________________________________________
faee3b18 245void AliHFEtools::DestroyDefaultPID(){
246 //
247 // Destroy default PID object if existing
248 //
249 if(fgDefaultPID) delete fgDefaultPID;
250 fgDefaultPID = NULL;
3a72645a 251 if(fgDefaultPIDaod) delete fgDefaultPIDaod;
252 fgDefaultPIDaod = NULL;
253}
254
255//__________________________________________
256Int_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);
e3ae862b 263 pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
3a72645a 264 } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
265 const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
e3ae862b 266 pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
3a72645a 267 }
268 return pdg;
269}
270
271//__________________________________________
272Int_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;
faee3b18 286}