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