]>
Commit | Line | Data |
---|---|---|
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" | |
a8ef1999 | 26 | #include "TH1D.h" |
70da6c5a | 27 | #include "TH2.h" |
a8ef1999 | 28 | #include "TGraph.h" |
29 | #include "TGraphErrors.h" | |
70da6c5a | 30 | #include "THnSparse.h" |
31 | #include "TAxis.h" | |
a8ef1999 | 32 | #include "TMath.h" |
33 | #include "TString.h" | |
70da6c5a | 34 | |
35 | #include "AliAODMCParticle.h" | |
3a72645a | 36 | #include "AliAODpidUtil.h" |
faee3b18 | 37 | #include "AliESDpid.h" |
70da6c5a | 38 | #include "AliLog.h" |
3a72645a | 39 | #include "AliTPCPIDResponse.h" |
faee3b18 | 40 | #include "AliTOFPIDResponse.h" |
70da6c5a | 41 | |
42 | #include "AliHFEtools.h" | |
43 | ||
44 | ClassImp(AliHFEtools) | |
45 | ||
8c1c76e9 | 46 | AliPIDResponse *AliHFEtools::fgDefaultPID = NULL; |
69ac0e6f | 47 | Int_t AliHFEtools::fgLogLevel = 0; |
faee3b18 | 48 | |
70da6c5a | 49 | //__________________________________________ |
50 | AliHFEtools::AliHFEtools(): | |
51 | TObject() | |
52 | { | |
53 | } | |
54 | ||
55 | //__________________________________________ | |
56 | Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){ | |
57 | // | |
58 | // Helper function for linearly binned array | |
59 | // | |
60 | Double_t *bins = new Double_t[nBins + 1]; | |
61 | Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins); | |
62 | bins[0] = ymin; | |
63 | for(Int_t ibin = 1; ibin <= nBins; ibin++) | |
64 | bins[ibin] = bins[ibin-1] + stepsize; | |
65 | return bins; | |
66 | } | |
67 | ||
68 | //__________________________________________ | |
69 | Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){ | |
70 | // | |
71 | // Helper function for logartimically binned array | |
72 | // | |
73 | Double_t *bins = new Double_t[nBins+1]; | |
74 | bins[0] = ymin; | |
75 | Double_t factor = TMath::Power(ymax/ymin, 1./nBins); | |
76 | for(Int_t ibin = 1; ibin <= nBins; ibin++){ | |
77 | bins[ibin] = factor * bins[ibin-1]; | |
78 | } | |
79 | return bins; | |
80 | } | |
81 | ||
82 | //_________________________________________ | |
83 | Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){ | |
84 | ||
85 | // | |
86 | // converts the axis (defined by the dimension) of THx or THnSparse | |
87 | // object to Log scale. Number of bins and bin min and bin max are preserved | |
88 | // | |
89 | ||
90 | ||
91 | if(!o){ | |
92 | AliError("Input histogram is null pointer"); | |
93 | return kFALSE; | |
94 | } | |
95 | ||
96 | TAxis *axis = 0x0; | |
97 | if(o->InheritsFrom("TH1")){ | |
bf892a6a | 98 | TH1 *h1 = dynamic_cast<TH1F*>(o); |
99 | if(h1) axis = h1->GetXaxis(); | |
70da6c5a | 100 | } |
101 | else if(o->InheritsFrom("TH2")){ | |
bf892a6a | 102 | TH2 *h2 = dynamic_cast<TH2F*>(o); |
103 | if(h2 && 0 == dim){ | |
104 | axis = h2->GetXaxis(); | |
70da6c5a | 105 | } |
bf892a6a | 106 | else if(h2 && 1 == dim){ |
107 | axis = h2->GetYaxis(); | |
70da6c5a | 108 | } |
109 | else{ | |
110 | AliError("Only dim = 0 or 1 possible for TH2F"); | |
111 | } | |
112 | } | |
113 | else if(o->InheritsFrom("THnSparse")){ | |
bf892a6a | 114 | THnSparse *hs = dynamic_cast<THnSparse*>(o); |
115 | if(hs) axis = hs->GetAxis(dim); | |
70da6c5a | 116 | } |
117 | else{ | |
118 | AliError("Type of input object not recognized, please check your code or update this finction"); | |
119 | return kFALSE; | |
120 | } | |
121 | if(!axis){ | |
122 | AliError(Form("Axis '%d' could not be identified in the object \n", dim)); | |
123 | return kFALSE; | |
124 | } | |
125 | ||
126 | Int_t bins = axis->GetNbins(); | |
127 | ||
128 | Double_t from = axis->GetXmin(); | |
129 | if(from <= 0){ | |
130 | AliError(Form("Log binning not possible for this axis [min = %f]\n", from)); | |
131 | } | |
132 | Double_t to = axis->GetXmax(); | |
133 | Double_t *newBins = new Double_t[bins+1]; | |
134 | newBins[0] = from; | |
135 | Double_t factor = TMath::Power(to/from, 1./bins); | |
136 | for(Int_t i=1; i<=bins; ++i){ | |
137 | newBins[i] = factor * newBins[i-1]; | |
138 | } | |
139 | axis->Set(bins, newBins); | |
bf892a6a | 140 | delete[] newBins; |
70da6c5a | 141 | |
142 | return kTRUE; | |
143 | } | |
144 | ||
145 | //__________________________________________ | |
3a72645a | 146 | Float_t AliHFEtools::GetRapidity(const TParticle *part){ |
70da6c5a | 147 | // |
148 | // return rapidity | |
149 | // | |
150 | Float_t rapidity; | |
151 | if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999; | |
152 | else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); | |
153 | return rapidity; | |
154 | } | |
155 | ||
156 | //__________________________________________ | |
3a72645a | 157 | Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){ |
70da6c5a | 158 | // return rapidity |
159 | ||
160 | Float_t rapidity; | |
161 | if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999; | |
162 | else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz()))); | |
163 | return rapidity; | |
164 | } | |
165 | ||
faee3b18 | 166 | //__________________________________________ |
8c1c76e9 | 167 | AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){ |
faee3b18 | 168 | // |
169 | // Get the default PID as singleton instance | |
170 | // | |
171 | if(!fgDefaultPID){ | |
8c1c76e9 | 172 | |
173 | if(isESD) fgDefaultPID = new AliESDpid; | |
174 | else fgDefaultPID = new AliAODpidUtil; | |
faee3b18 | 175 | Double_t tres = isMC ? 80. : 130.; |
176 | fgDefaultPID->GetTOFResponse().SetTimeResolution(tres); | |
177 | ||
178 | // TPC Bethe Bloch parameters | |
179 | Double_t alephParameters[5]; | |
180 | if(isMC){ | |
181 | // simulation | |
182 | alephParameters[0] = 2.15898e+00/50.; | |
183 | alephParameters[1] = 1.75295e+01; | |
184 | alephParameters[2] = 3.40030e-09; | |
185 | alephParameters[3] = 1.96178e+00; | |
186 | alephParameters[4] = 3.91720e+00; | |
187 | } else { | |
188 | alephParameters[0] = 0.0283086/0.97; | |
189 | //alephParameters[0] = 0.0283086; | |
190 | alephParameters[1] = 2.63394e+01; | |
191 | alephParameters[2] = 5.04114e-11; | |
192 | alephParameters[3] = 2.12543e+00; | |
193 | alephParameters[4] = 4.88663e+00; | |
194 | } | |
195 | fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]); | |
67fe7bd0 | 196 | fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04); |
faee3b18 | 197 | |
198 | } | |
199 | if(fgLogLevel){ | |
67fe7bd0 | 200 | printf("Error - You are using the default PID: You should use the PID coming from the tender\n"); |
201 | printf("Error - Arrrrrrrrr...\n"); | |
202 | printf("Error - Please rethink your program logic. Using default PID is really dangerous\n"); | |
203 | printf("Error - TOF PID is adapted to Monte Carlo\n"); | |
faee3b18 | 204 | } |
205 | return fgDefaultPID; | |
206 | } | |
207 | ||
3a72645a | 208 | |
faee3b18 | 209 | //__________________________________________ |
210 | void AliHFEtools::DestroyDefaultPID(){ | |
211 | // | |
212 | // Destroy default PID object if existing | |
213 | // | |
214 | if(fgDefaultPID) delete fgDefaultPID; | |
215 | fgDefaultPID = NULL; | |
3a72645a | 216 | } |
217 | ||
218 | //__________________________________________ | |
219 | Int_t AliHFEtools::GetPdg(const AliVParticle *track){ | |
220 | // | |
221 | // Get MC PDG code (only MC particles supported) | |
222 | // | |
223 | Int_t pdg = 0; | |
224 | if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){ | |
225 | const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track); | |
e3ae862b | 226 | pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0; |
3a72645a | 227 | } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){ |
228 | const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track); | |
e3ae862b | 229 | pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0; |
3a72645a | 230 | } |
231 | return pdg; | |
232 | } | |
233 | ||
234 | //__________________________________________ | |
235 | Int_t AliHFEtools::PDG2AliPID(Int_t pdg){ | |
236 | // | |
237 | // Helper function to convert MC PID codes into AliPID codes | |
238 | // | |
239 | Int_t pid = -1; | |
240 | switch(TMath::Abs(pdg)){ | |
241 | case 11: pid = AliPID::kElectron; break; | |
242 | case 13: pid = AliPID::kMuon; break; | |
243 | case 211: pid = AliPID::kPion; break; | |
244 | case 321: pid = AliPID::kKaon; break; | |
245 | case 2212: pid = AliPID::kProton; break; | |
246 | default: pid = -1; break; | |
247 | }; | |
248 | return pid; | |
faee3b18 | 249 | } |
a8ef1999 | 250 | |
251 | //__________________________________________ | |
252 | TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize) | |
253 | { | |
254 | // Creates a TH1D from TGraph g. The binwidth of the first bin has to | |
255 | // specified. The others bins are calculated automatically. Supports also Graphs | |
256 | // with non constant x steps. The axis of the Graph can be exchanged if | |
257 | // exchange=kTRUE (modifies the Graph). | |
258 | ||
259 | TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize); | |
260 | if( result == 0) return result; | |
261 | ||
262 | //------------------------------------------ | |
263 | // setup the final hist | |
264 | Int_t nBinX = g->GetN(); | |
265 | Double_t* err = g->GetEY(); // error y is still ok even if exchanged | |
266 | for(Int_t i = 0; i < nBinX; i ++){ | |
267 | result->SetBinError(i + 1, err[i]); | |
268 | } | |
269 | if(exchange){ | |
270 | AliHFEtools::ExchangeXYGraph(g); // undo what has been done in GraphToHist | |
271 | AliHFEtools::ExchangeXYGraphErrors(g); // now exchange including errors | |
272 | } | |
273 | ||
274 | return result; | |
275 | } | |
276 | ||
277 | //__________________________________________ | |
278 | Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g) | |
279 | { | |
280 | // exchanges x-values and y-values. | |
281 | if(g==0) return kFALSE; | |
282 | Int_t nbin=g->GetN(); | |
283 | Double_t x,y; | |
284 | for(Int_t i = 0; i < nbin; i ++) | |
285 | { | |
286 | g->GetPoint(i,x,y); | |
287 | g->SetPoint(i,y,x); | |
288 | } | |
289 | ||
290 | return kTRUE; | |
291 | } | |
292 | ||
293 | //__________________________________________ | |
294 | Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g) | |
295 | { | |
296 | // exchanges x-values and y-values and | |
297 | // corresponding errors. | |
298 | if(g==0) return kFALSE; | |
299 | Int_t nbin=g->GetN(); | |
300 | Double_t x,y; | |
301 | Double_t ex,ey; | |
302 | for(Int_t i = 0; i < nbin; i ++) | |
303 | { | |
304 | g->GetPoint(i,x,y); | |
305 | ex = g->GetErrorX(i); | |
306 | ey = g->GetErrorY(i); | |
307 | g->SetPoint(i,y,x); | |
308 | g->SetPointError(i,ey,ex); | |
309 | } | |
310 | ||
311 | return kTRUE; | |
312 | ||
313 | } | |
314 | ||
315 | //__________________________________________ | |
316 | TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize) | |
317 | { | |
318 | // Creates a TH1D from TGraph g. The binwidth of the first bin has to | |
319 | // specified. The others bins are calculated automatically. Supports also Graphs | |
320 | // with non constant x steps. The axis of the Graph can be exchanged if | |
321 | // exchange=kTRUE (modifies the Graph). | |
322 | ||
323 | ||
324 | TH1D* result = 0; | |
325 | if(g == 0) return result; | |
326 | if(firstBinWidth == -1) return result; | |
327 | TString myname=""; | |
328 | myname = g->GetName(); | |
329 | myname += "_graph"; | |
330 | if(exchange) AliHFEtools::ExchangeXYGraph(g); | |
331 | ||
332 | Int_t nBinX = g->GetN(); | |
333 | Double_t* x = g->GetX(); | |
334 | Double_t* y = g->GetY(); | |
335 | ||
336 | if(nBinX < 1) return result; | |
337 | ||
338 | //------------------------------------------ | |
339 | // create the Matrix for the equation system | |
340 | // and init the values | |
341 | ||
342 | Int_t nDim = nBinX - 1; | |
343 | TMatrixD a(nDim,nDim); | |
344 | TMatrixD b(nDim,1); | |
345 | ||
346 | Double_t* aA = a.GetMatrixArray(); | |
347 | Double_t* aB = b.GetMatrixArray(); | |
348 | memset(aA,0,nDim * nDim * sizeof(Double_t)); | |
349 | memset(aB,0,nDim * sizeof(Double_t)); | |
350 | //------------------------------------------ | |
351 | ||
352 | //------------------------------------------ | |
353 | // setup equation system | |
354 | // width for 1st bin is given therefore | |
355 | // we shift bin parameter (column) by one to the left | |
356 | // to reduce the matrix size | |
357 | ||
358 | Double_t* xAxis = new Double_t [nBinX + 1]; | |
359 | Double_t* binW = new Double_t [nBinX ]; | |
360 | binW[0] = firstBinWidth; | |
361 | ||
362 | aB[0] = x[1] - x[0] - 0.5 * binW[0]; | |
363 | aA[0] = 0.5; | |
364 | ||
365 | for(Int_t col = 1; col < nDim ; col ++) | |
366 | { | |
367 | Int_t row = col; | |
368 | aB[col] = x[col + 1] - x[ col ]; | |
369 | aA[row * nDim + col - 1 ] = 0.5; | |
370 | aA[row * nDim + col ] = 0.5; | |
371 | } | |
372 | //------------------------------------------ | |
373 | ||
374 | //------------------------------------------ | |
375 | // solve the equations | |
376 | a.Invert(); | |
377 | TMatrixD c = a * b; | |
378 | //------------------------------------------ | |
379 | ||
380 | //------------------------------------------ | |
381 | // calculate the bin boundaries | |
382 | xAxis[0] = x[0] - 0.5 * binW[0]; | |
383 | memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t)); | |
384 | for(Int_t col = 0; col < nBinX ; col ++) { | |
385 | xAxis[col + 1] = x[col] + 0.5 * binW[col]; | |
386 | } | |
387 | //------------------------------------------ | |
388 | ||
389 | //------------------------------------------ | |
390 | // setup the final hist | |
391 | result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis); | |
392 | for(Int_t i = 0; i < nBinX; i ++){ | |
393 | result->SetBinContent(i + 1, y[i]); | |
394 | } | |
395 | result->SetMarkerColor(markercolor); | |
396 | result->SetMarkerStyle(markerstyle); | |
397 | result->SetMarkerSize(markersize); | |
398 | //------------------------------------------ | |
399 | ||
400 | delete [] xAxis; | |
401 | delete [] binW; | |
402 | ||
403 | ||
404 | return result; | |
405 | } |