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