1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Toolkit containing various usefull things
17 // Usable everywhere in the hfe software package
18 // For more information see the cxx file
21 // All authors of the HFE group
25 #include <TParticle.h>
30 #include "TGraphErrors.h"
31 #include "THnSparse.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODpidUtil.h"
38 #include "AliESDpid.h"
40 #include "AliTPCPIDResponse.h"
41 #include "AliTOFPIDResponse.h"
43 #include "AliHFEtools.h"
47 AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
48 Int_t AliHFEtools::fgLogLevel = 0;
50 //__________________________________________
51 AliHFEtools::AliHFEtools():
56 //__________________________________________
57 Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
59 // Helper function for linearly binned array
61 Double_t *bins = new Double_t[nBins + 1];
62 Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
64 for(Int_t ibin = 1; ibin <= nBins; ibin++)
65 bins[ibin] = bins[ibin-1] + stepsize;
69 //__________________________________________
70 void AliHFEtools::FillLinearBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
72 // Helper function for linearly binned array
74 Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
76 for(Int_t ibin = 1; ibin <= nBins; ibin++)
77 bins[ibin] = bins[ibin-1] + stepsize;
80 //__________________________________________
81 Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
83 // Helper function for logartimically binned array
85 Double_t *bins = new Double_t[nBins+1];
87 Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
88 for(Int_t ibin = 1; ibin <= nBins; ibin++){
89 bins[ibin] = factor * bins[ibin-1];
94 //__________________________________________
95 void AliHFEtools::FillLogarithmicBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
97 // Helper function for logartimically binned array
100 Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
101 for(Int_t ibin = 1; ibin <= nBins; ibin++)
102 bins[ibin] = factor * bins[ibin-1];
105 //_________________________________________
106 Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
109 // converts the axis (defined by the dimension) of THx or THnSparse
110 // object to Log scale. Number of bins and bin min and bin max are preserved
115 AliError("Input histogram is null pointer");
120 if(o->InheritsFrom("TH1")){
121 TH1 *h1 = dynamic_cast<TH1F*>(o);
122 if(h1) axis = h1->GetXaxis();
124 else if(o->InheritsFrom("TH2")){
125 TH2 *h2 = dynamic_cast<TH2F*>(o);
127 axis = h2->GetXaxis();
129 else if(h2 && 1 == dim){
130 axis = h2->GetYaxis();
133 AliError("Only dim = 0 or 1 possible for TH2F");
136 else if(o->InheritsFrom("THnSparse")){
137 THnSparse *hs = dynamic_cast<THnSparse*>(o);
138 if(hs) axis = hs->GetAxis(dim);
141 AliError("Type of input object not recognized, please check your code or update this finction");
145 AliError(Form("Axis '%d' could not be identified in the object \n", dim));
149 Int_t bins = axis->GetNbins();
151 Double_t from = axis->GetXmin();
153 AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
155 Double_t to = axis->GetXmax();
156 TArrayD newBins(bins+1);
158 Double_t factor = TMath::Power(to/from, 1./bins);
159 for(Int_t i=1; i<=bins; ++i){
160 newBins[i] = factor * newBins[i-1];
162 axis->Set(bins, newBins.GetArray());
167 //__________________________________________
168 Float_t AliHFEtools::GetRapidity(const TParticle *part){
173 if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
174 else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
178 //__________________________________________
179 Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
183 if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
184 else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
188 //__________________________________________
189 AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
191 // Get the default PID as singleton instance
195 if(isESD) fgDefaultPID = new AliESDpid;
196 else fgDefaultPID = new AliAODpidUtil;
197 Double_t tres = isMC ? 80. : 130.;
198 fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
200 // TPC Bethe Bloch parameters
201 Double_t alephParameters[5];
204 alephParameters[0] = 2.15898e+00/50.;
205 alephParameters[1] = 1.75295e+01;
206 alephParameters[2] = 3.40030e-09;
207 alephParameters[3] = 1.96178e+00;
208 alephParameters[4] = 3.91720e+00;
210 alephParameters[0] = 0.0283086/0.97;
211 //alephParameters[0] = 0.0283086;
212 alephParameters[1] = 2.63394e+01;
213 alephParameters[2] = 5.04114e-11;
214 alephParameters[3] = 2.12543e+00;
215 alephParameters[4] = 4.88663e+00;
217 fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
218 fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
222 printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
223 printf("Error - Arrrrrrrrr...\n");
224 printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
225 printf("Error - TOF PID is adapted to Monte Carlo\n");
231 //__________________________________________
232 void AliHFEtools::DestroyDefaultPID(){
234 // Destroy default PID object if existing
236 if(fgDefaultPID) delete fgDefaultPID;
240 //__________________________________________
241 Int_t AliHFEtools::GetPdg(const AliVParticle *track){
243 // Get MC PDG code (only MC particles supported)
246 if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
247 const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track);
248 pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
249 } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
250 const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
251 pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
256 //__________________________________________
257 Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
259 // Helper function to convert MC PID codes into AliPID codes
262 switch(TMath::Abs(pdg)){
263 case 11: pid = AliPID::kElectron; break;
264 case 13: pid = AliPID::kMuon; break;
265 case 211: pid = AliPID::kPion; break;
266 case 321: pid = AliPID::kKaon; break;
267 case 2212: pid = AliPID::kProton; break;
268 default: pid = -1; break;
273 //__________________________________________
274 TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
276 // Creates a TH1D from TGraph g. The binwidth of the first bin has to
277 // specified. The others bins are calculated automatically. Supports also Graphs
278 // with non constant x steps. The axis of the Graph can be exchanged if
279 // exchange=kTRUE (modifies the Graph).
281 TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
282 if( result == 0) return result;
284 //------------------------------------------
285 // setup the final hist
286 Int_t nBinX = g->GetN();
287 Double_t* err = g->GetEY(); // error y is still ok even if exchanged
288 for(Int_t i = 0; i < nBinX; i ++){
289 result->SetBinError(i + 1, err[i]);
292 AliHFEtools::ExchangeXYGraph(g); // undo what has been done in GraphToHist
293 AliHFEtools::ExchangeXYGraphErrors(g); // now exchange including errors
299 //__________________________________________
300 Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
302 // exchanges x-values and y-values.
303 if(g==0) return kFALSE;
304 Int_t nbin=g->GetN();
306 for(Int_t i = 0; i < nbin; i ++)
315 //__________________________________________
316 Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
318 // exchanges x-values and y-values and
319 // corresponding errors.
320 if(g==0) return kFALSE;
321 Int_t nbin=g->GetN();
324 for(Int_t i = 0; i < nbin; i ++)
327 ex = g->GetErrorX(i);
328 ey = g->GetErrorY(i);
330 g->SetPointError(i,ey,ex);
337 //__________________________________________
338 TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
340 // Creates a TH1D from TGraph g. The binwidth of the first bin has to
341 // specified. The others bins are calculated automatically. Supports also Graphs
342 // with non constant x steps. The axis of the Graph can be exchanged if
343 // exchange=kTRUE (modifies the Graph).
347 if(g == 0) return result;
348 if(firstBinWidth == -1) return result;
350 myname = g->GetName();
352 if(exchange) AliHFEtools::ExchangeXYGraph(g);
354 Int_t nBinX = g->GetN();
355 Double_t* x = g->GetX();
356 Double_t* y = g->GetY();
358 if(nBinX < 1) return result;
360 //------------------------------------------
361 // create the Matrix for the equation system
362 // and init the values
364 Int_t nDim = nBinX - 1;
365 TMatrixD a(nDim,nDim);
368 Double_t* aA = a.GetMatrixArray();
369 Double_t* aB = b.GetMatrixArray();
370 memset(aA,0,nDim * nDim * sizeof(Double_t));
371 memset(aB,0,nDim * sizeof(Double_t));
372 //------------------------------------------
374 //------------------------------------------
375 // setup equation system
376 // width for 1st bin is given therefore
377 // we shift bin parameter (column) by one to the left
378 // to reduce the matrix size
380 Double_t* xAxis = new Double_t [nBinX + 1];
381 Double_t* binW = new Double_t [nBinX ];
382 binW[0] = firstBinWidth;
384 aB[0] = x[1] - x[0] - 0.5 * binW[0];
387 for(Int_t col = 1; col < nDim ; col ++)
390 aB[col] = x[col + 1] - x[ col ];
391 aA[row * nDim + col - 1 ] = 0.5;
392 aA[row * nDim + col ] = 0.5;
394 //------------------------------------------
396 //------------------------------------------
397 // solve the equations
400 //------------------------------------------
402 //------------------------------------------
403 // calculate the bin boundaries
404 xAxis[0] = x[0] - 0.5 * binW[0];
405 memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
406 for(Int_t col = 0; col < nBinX ; col ++) {
407 xAxis[col + 1] = x[col] + 0.5 * binW[col];
409 //------------------------------------------
411 //------------------------------------------
412 // setup the final hist
413 result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
414 for(Int_t i = 0; i < nBinX; i ++){
415 result->SetBinContent(i + 1, y[i]);
417 result->SetMarkerColor(markercolor);
418 result->SetMarkerStyle(markerstyle);
419 result->SetMarkerSize(markersize);
420 //------------------------------------------