[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPID.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Container for the distributions of dE/dx and the time bin of the     //
// max. cluster for electrons and pions                                 //
//                                                                      //
// Authors:                                                             //
//   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
//   Alex Bercuci <a.bercuci@gsi.de>                                    //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TROOT.h>

#include "AliLog.h"
#include "AliPID.h"
#include "AliESD.h"
#include "AliESDtrack.h"

#include "AliTRDCalPID.h"
#include "AliTRDcalibDB.h"

ClassImp(AliTRDCalPID)

Char_t* AliTRDCalPID::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};

Char_t* AliTRDCalPID::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};

Float_t AliTRDCalPID::fTrackMomentum[kNMom] = {0.6,  0.8,  1.0,  1.5,  2.0,  3.0,  4.0,  5.0,  6.0,  8.0,  10.0};
  
Float_t AliTRDCalPID::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};

    
//_________________________________________________________________________
AliTRDCalPID::AliTRDCalPID()
  :TNamed("pid", "PID for TRD")
  ,fMeanChargeRatio(0)
  ,fHistdEdx(0x0)
  ,fHistTimeBin(0x0)
{
  //
  //  The Default constructor
  //

  Init();

}

//_________________________________________________________________________
AliTRDCalPID::AliTRDCalPID(const Text_t *name, const Text_t *title) 
  :TNamed(name,title)
  ,fMeanChargeRatio(0)
  ,fHistdEdx(0x0)
  ,fHistTimeBin(0x0)
{
  //
  //  The main constructor
  //
  
  Init();

}

//_____________________________________________________________________________
AliTRDCalPID::AliTRDCalPID(const AliTRDCalPID &c) 
  :TNamed(c)
  ,fMeanChargeRatio(c.fMeanChargeRatio)
  ,fHistdEdx(0x0)
  ,fHistTimeBin(0x0)
{
  //
  // Copy constructor
  //

  if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
  
}

//_________________________________________________________________________
AliTRDCalPID::~AliTRDCalPID()
{
  //
  // Destructor
  //
  
  CleanUp();

}

//_________________________________________________________________________
void AliTRDCalPID::CleanUp()
{
  //
  // Delets all newly created objects
  //

  if (fHistdEdx) {
    delete fHistdEdx;
    fHistdEdx = 0x0;
  }
  
  if (fHistTimeBin) {
    delete fHistTimeBin;
    fHistTimeBin = 0x0;
  }
}

//_____________________________________________________________________________
AliTRDCalPID &AliTRDCalPID::operator=(const AliTRDCalPID &c)
{
  //
  // Assignment operator
  //

  if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
  return *this;

}

//_____________________________________________________________________________
void AliTRDCalPID::Copy(TObject &c) const
{
  //
  // Copy function
  //

  AliTRDCalPID& target = (AliTRDCalPID &) c;
  
  target.CleanUp();
  
  target.fMeanChargeRatio = fMeanChargeRatio;

  if (fHistdEdx) {
    target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
  }
  if (fHistTimeBin) {
    target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
  }

  TObject::Copy(c);

}

//_________________________________________________________________________
void AliTRDCalPID::Init()
{
  //
  // Initialization
  //

  fHistdEdx    = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
  fHistdEdx->SetOwner();
  fHistTimeBin = new TObjArray(2 * kNMom);
  fHistTimeBin->SetOwner();  

	// Initialization of estimator at object instantiation because late
	// initialization in function GetProbability() is not working due to
	// constantness of this function. 
	// fEstimator = new AliTRDCalPIDRefMaker();
	
  // ADC Gain normalization
  fMeanChargeRatio = 1.0;
}

//_________________________________________________________________________
Bool_t AliTRDCalPID::LoadLQReferences(Char_t *refFile)
{
  //
  // Read the TRD dEdx histograms.
  //

	Int_t nTimeBins = 22;
	// Get number of time bins from CDB
	AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
	if(!calibration){
		AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
	}else{
		if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
		else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
	}

	
  // Read histogram Root file  
  TFile *histFile = TFile::Open(refFile, "READ");
  if (!histFile || !histFile->IsOpen()) {
    AliError(Form("Opening TRD histgram file %s failed", refFile));
    return kFALSE;
  }
  gROOT->cd();

  // Read histograms
  for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
    for (Int_t imom = 0; imom < kNMom; imom++){
			TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
			hist->Scale(1./hist->Integral());
			fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));

			if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;

			TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
			if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, nTimeBins));
			ht->Scale(1./ht->Integral());
			fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
		}
  }
  
  histFile->Close();
  delete histFile;
  
  // Number of bins and bin size
  //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
  //fNbins   = hist->GetNbinsX();
  //fBinSize = hist->GetBinWidth(1);
  
  return kTRUE;

}

// //_________________________________________________________________________
// Double_t  AliTRDCalPID::GetMean(Int_t k, Int_t ip) const
// {
//   //
//   // Gets mean of de/dx dist. of e
//   //
// 
//   AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
//               ,fpartName[k]
//               ,fTrackMomentum[ip]));
//   if (k < 0 || k > AliPID::kSPECIES) {
//     return 0;
//   }
// 
//   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
// 
// }
// 
// //_________________________________________________________________________
// Double_t  AliTRDCalPID::GetNormalization(Int_t k, Int_t ip) const
// {
//   //
//   // Gets Normalization of de/dx dist. of e
//   //
// 
//   AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
//               ,fpartName[k]
//               ,fTrackMomentum[ip]));
//   if (k < 0 || k > AliPID::kSPECIES) {
//     return 0;
//   }
//   
//   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
// 
// }

//_________________________________________________________________________
TH1* AliTRDCalPID::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
{
  //
  // Returns one selected dEdx histogram
  //

  if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
	if(ip<0 || ip>= kNMom ) return 0x0;

	AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
  
  return (TH1*)fHistdEdx->At(GetHistID(k, ip));

}

//_________________________________________________________________________
TH1* AliTRDCalPID::GetHistogramT(Int_t k, Int_t ip) const
{
  //
  // Returns one selected time bin max histogram
  //

  if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
	if(ip<0 || ip>= kNMom ) return 0x0;
	  
 	AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));

	return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*kNMom+ip);
}


//_________________________________________________________________________
Double_t AliTRDCalPID::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const
{
  //
	// Core function of AliTRDCalPID class for calculating the
	// likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
	// given momentum "mom" and a given dE/dx (slice "dedx") yield per
	// layer
  //

	if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
		
	//Double_t dedx   = dedx1/fMeanChargeRatio;
	
	// find the interval in momentum and track segment length which applies for this data
	Int_t ilength = 1;
  while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
		ilength++;
	}
	Int_t imom = 1;
  while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;

	Int_t nbinsx, nbinsy;
	TAxis *ax = 0x0, *ay = 0x0;
	Double_t LQ1, LQ2;
	Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
	TH2 *hist = 0x0;
	if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom-1/*, ilength*/)))){
		AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
		AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom-1)));
		return 0.;
	}
	ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
	ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
	Float_t x = dedx[0]+dedx[1], y = dedx[2];
        Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
	Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
	if(kX)
		if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y)); 
    //fEstimator->Estimate2D2(hist, x, y);
		else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
	else
		if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
	 	else LQ1 = hist->GetBinContent(nbinsx, nbinsy);


	if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
		AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
		AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
		return LQ1;
	}
	if(kX)
		if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y)); 
    //fEstimator->Estimate2D2(hist, x, y);
		else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
	else
		if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
	 	else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
	
        // return interpolation over momentum binning
        if(mom < fTrackMomentum[0]) return LQ1;
        else if(mom > fTrackMomentum[kNMom-1]) return LQ2;
        else return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);

}

//_________________________________________________________________________
Double_t AliTRDCalPID::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
{
  //
  // Gets the Probability of having timbin at a given momentum (mom)
  // and particle type (spec) (0 for e) and (2 for pi)
  // from the precalculated timbin distributions 
  //
  
	if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;

  Int_t iTBin = timbin+1;
  
  // Everything which is not an electron counts as a pion for time bin max
  //if(spec != AliPID::kElectron) spec = AliPID::kPion;
  

	Int_t imom = 1;
  while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;

	Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
	TH1F *hist = 0x0;
	if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom-1))){
		AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
		AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom-1));
		return 0.;
	}
	Double_t LQ1 = hist->GetBinContent(iTBin);

	if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom))){
		AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
		AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom));
		return LQ1;
	}
	Double_t LQ2 = hist->GetBinContent(iTBin);

	// return interpolation over momentum binning
  return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
}
Commit	Line	Data
720a0a16	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	/* $Id$ */
	17
	18	//////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Container for the distributions of dE/dx and the time bin of the //
	21	// max. cluster for electrons and pions //
	22	// //
	23	// Authors: //
	24	// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> //
	25	// Alex Bercuci <a.bercuci@gsi.de> //
	26	// //
	27	//////////////////////////////////////////////////////////////////////////
	28
	29	#include <TH1F.h>
	30	#include <TH2F.h>
	31	#include <TFile.h>
	32	#include <TROOT.h>
	33
	34	#include "AliLog.h"
	35	#include "AliPID.h"
	36	#include "AliESD.h"
	37	#include "AliESDtrack.h"
	38
	39	#include "AliTRDCalPID.h"
	40	#include "AliTRDcalibDB.h"
	41
	42	ClassImp(AliTRDCalPID)
	43
	44	Char_t* AliTRDCalPID::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
	45
	46	Char_t* AliTRDCalPID::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
	47
	48	Float_t AliTRDCalPID::fTrackMomentum[kNMom] = {0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0};
	49
	50	Float_t AliTRDCalPID::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};
	51
	52
	53	//_________________________________________________________________________
	54	AliTRDCalPID::AliTRDCalPID()
	55	:TNamed("pid", "PID for TRD")
	56	,fMeanChargeRatio(0)
	57	,fHistdEdx(0x0)
	58	,fHistTimeBin(0x0)
	59	{
	60	//
	61	// The Default constructor
	62	//
	63
	64	Init();
65
66	}
67
68	//_________________________________________________________________________
69	AliTRDCalPID::AliTRDCalPID(const Text_t name, const Text_t title)
70	:TNamed(name,title)
71	,fMeanChargeRatio(0)
72	,fHistdEdx(0x0)
73	,fHistTimeBin(0x0)
74	{
75	//
76	// The main constructor
77	//
78
79	Init();
80
81	}
82
83	//_____________________________________________________________________________
84	AliTRDCalPID::AliTRDCalPID(const AliTRDCalPID &c)
85	:TNamed(c)
86	,fMeanChargeRatio(c.fMeanChargeRatio)
87	,fHistdEdx(0x0)
88	,fHistTimeBin(0x0)
89	{
90	//
91	// Copy constructor
92	//
93
94	if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
95
96	}
97
98	//_________________________________________________________________________
99	AliTRDCalPID::~AliTRDCalPID()
100	{
101	//
102	// Destructor
103	//
104
105	CleanUp();
106
107	}
108
109	//_________________________________________________________________________
110	void AliTRDCalPID::CleanUp()
111	{
112	//
113	// Delets all newly created objects
114	//
115
116	if (fHistdEdx) {
117	delete fHistdEdx;
118	fHistdEdx = 0x0;
119	}
120
121	if (fHistTimeBin) {
122	delete fHistTimeBin;
123	fHistTimeBin = 0x0;
124	}
125	}
126
127	//_____________________________________________________________________________
128	AliTRDCalPID &AliTRDCalPID::operator=(const AliTRDCalPID &c)
129	{
130	//
131	// Assignment operator
132	//
133
134	if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
135	return *this;
136
137	}
138
139	//_____________________________________________________________________________
140	void AliTRDCalPID::Copy(TObject &c) const
141	{
142	//
143	// Copy function
144	//
145
146	AliTRDCalPID& target = (AliTRDCalPID &) c;
147
148	target.CleanUp();
149
150	target.fMeanChargeRatio = fMeanChargeRatio;
151
152	if (fHistdEdx) {
153	target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
154	}
155	if (fHistTimeBin) {
156	target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
157	}
158
159	TObject::Copy(c);
160
161	}
162
163	//_________________________________________________________________________
164	void AliTRDCalPID::Init()
165	{
166	//
167	// Initialization
168	//
169
170	fHistdEdx = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
171	fHistdEdx->SetOwner();
172	fHistTimeBin = new TObjArray(2 * kNMom);
173	fHistTimeBin->SetOwner();
174
175	// Initialization of estimator at object instantiation because late
176	// initialization in function GetProbability() is not working due to
177	// constantness of this function.
178	// fEstimator = new AliTRDCalPIDRefMaker();
179
180	// ADC Gain normalization
181	fMeanChargeRatio = 1.0;
182	}
183
184	//_________________________________________________________________________
185	Bool_t AliTRDCalPID::LoadLQReferences(Char_t *refFile)
186	{
187	//
188	// Read the TRD dEdx histograms.
189	//
190
191	Int_t nTimeBins = 22;
192	// Get number of time bins from CDB
193	AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
194	if(!calibration){
195	AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
196	}else{
197	if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
198	else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
199	}
200
201
202	// Read histogram Root file
203	TFile *histFile = TFile::Open(refFile, "READ");
204	if (!histFile \|\| !histFile->IsOpen()) {
205	AliError(Form("Opening TRD histgram file %s failed", refFile));
206	return kFALSE;
207	}
208	gROOT->cd();
209
210	// Read histograms
211	for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
212	for (Int_t imom = 0; imom < kNMom; imom++){
213	TH2D* hist = (TH2D)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/, ilength*/))->Clone();
214	hist->Scale(1./hist->Integral());
215	fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
216
217	if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
218
219	TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
220	if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, nTimeBins));
221	ht->Scale(1./ht->Integral());
222	fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
223	}
224	}
225
226	histFile->Close();
227	delete histFile;
228
229	// Number of bins and bin size
230	//TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
231	//fNbins = hist->GetNbinsX();
232	//fBinSize = hist->GetBinWidth(1);
233
234	return kTRUE;
235
236	}
237
238	// //_________________________________________________________________________
239	// Double_t AliTRDCalPID::GetMean(Int_t k, Int_t ip) const
240	// {
241	// //
242	// // Gets mean of de/dx dist. of e
243	// //
244	//
245	// AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
246	// ,fpartName[k]
247	// ,fTrackMomentum[ip]));
248	// if (k < 0 \|\| k > AliPID::kSPECIES) {
249	// return 0;
250	// }
251	//
252	// return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
253	//
254	// }
255	//
256	// //_________________________________________________________________________
257	// Double_t AliTRDCalPID::GetNormalization(Int_t k, Int_t ip) const
258	// {
259	// //
260	// // Gets Normalization of de/dx dist. of e
261	// //
262	//
263	// AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
264	// ,fpartName[k]
265	// ,fTrackMomentum[ip]));
266	// if (k < 0 \|\| k > AliPID::kSPECIES) {
267	// return 0;
268	// }
269	//
270	// return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
271	//
272	// }
273
274	//_________________________________________________________________________
275	TH1* AliTRDCalPID::GetHistogram(Int_t k, Int_t ip/, Int_t il/) const
276	{
277	//
278	// Returns one selected dEdx histogram
279	//
280
281	if (k < 0 \|\| k >= AliPID::kSPECIES) return 0x0;
282	if(ip<0 \|\| ip>= kNMom ) return 0x0;
283
284	AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
285
286	return (TH1*)fHistdEdx->At(GetHistID(k, ip));
287
288	}
289
290	//_________________________________________________________________________
291	TH1* AliTRDCalPID::GetHistogramT(Int_t k, Int_t ip) const
292	{
293	//
294	// Returns one selected time bin max histogram
295	//
296
297	if (k < 0 \|\| k >= AliPID::kSPECIES) return 0x0;
298	if(ip<0 \|\| ip>= kNMom ) return 0x0;
299
300	AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
301
302	return (TH1)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)kNMom+ip);
303	}
304
305
306
307	//_________________________________________________________________________
308	Double_t AliTRDCalPID::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const
309	{
310	//
311	// Core function of AliTRDCalPID class for calculating the
312	// likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
313	// given momentum "mom" and a given dE/dx (slice "dedx") yield per
314	// layer
315	//
316
317	if (spec < 0 \|\| spec >= AliPID::kSPECIES) return 0.;
318
319	//Double_t dedx = dedx1/fMeanChargeRatio;
320
321	// find the interval in momentum and track segment length which applies for this data
322	Int_t ilength = 1;
323	while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
324	ilength++;
325	}
326	Int_t imom = 1;
327	while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
79143cbf	328
720a0a16	329	Int_t nbinsx, nbinsy;
	330	TAxis ax = 0x0, ay = 0x0;
	331	Double_t LQ1, LQ2;
	332	Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
	333	TH2 *hist = 0x0;
	334	if(!(hist = (TH2D)fHistdEdx->At(GetHistID(spec, imom-1/, ilength*/)))){
	335	AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
	336	AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom-1)));
	337	return 0.;
	338	}
	339	ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
	340	ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
	341	Float_t x = dedx[0]+dedx[1], y = dedx[2];
79143cbf	342	Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
720a0a16	343	Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
	344	if(kX)
	345	if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y));
	346	//fEstimator->Estimate2D2(hist, x, y);
	347	else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
	348	else
	349	if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
	350	else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
	351
	352
	353	if(!(hist = (TH2D)fHistdEdx->At(GetHistID(spec, imom/, ilength*/)))){
	354	AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
	355	AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
	356	return LQ1;
	357	}
	358	if(kX)
	359	if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y));
	360	//fEstimator->Estimate2D2(hist, x, y);
	361	else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
	362	else
	363	if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
	364	else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
720a0a16	365
79143cbf	366	// return interpolation over momentum binning
	367	if(mom < fTrackMomentum[0]) return LQ1;
	368	else if(mom > fTrackMomentum[kNMom-1]) return LQ2;
	369	else return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
720a0a16	370
	371	}
	372
	373	//_________________________________________________________________________
	374	Double_t AliTRDCalPID::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
	375	{
	376	//
	377	// Gets the Probability of having timbin at a given momentum (mom)
	378	// and particle type (spec) (0 for e) and (2 for pi)
	379	// from the precalculated timbin distributions
	380	//
	381
	382	if (spec < 0 \|\| spec >= AliPID::kSPECIES) return 0.;
	383
	384	Int_t iTBin = timbin+1;
	385
	386	// Everything which is not an electron counts as a pion for time bin max
	387	//if(spec != AliPID::kElectron) spec = AliPID::kPion;
	388
	389
	390
	391	Int_t imom = 1;
	392	while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
	393
	394	Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
	395	TH1F *hist = 0x0;
	396	if(!(hist = (TH1F) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)kNMom+imom-1))){
	397	AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
	398	AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom-1));
	399	return 0.;
	400	}
	401	Double_t LQ1 = hist->GetBinContent(iTBin);
	402
	403	if(!(hist = (TH1F) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)kNMom+imom))){
	404	AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
	405	AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom));
	406	return LQ1;
	407	}
	408	Double_t LQ2 = hist->GetBinContent(iTBin);
	409
	410	// return interpolation over momentum binning
	411	return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
	412	}
	413