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