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