[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPIDLQ.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.                  *
 **************************************************************************/


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

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

#include "AliLog.h"
#include "AliPID.h"

#include "AliTRDCalPIDLQ.h"

ClassImp(AliTRDCalPIDLQ)

Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
    
//_________________________________________________________________________
AliTRDCalPIDLQ::AliTRDCalPIDLQ()
  :TNamed()
  ,fNMom(0)
  ,fTrackMomentum(0)
  ,fMeanChargeRatio(0)
  ,fNbins(0)
  ,fBinSize(0)
  ,fHistdEdx(0)
  ,fHistTimeBin(0)
{
  //
  //  The Default constructor
  //
  
}

//_________________________________________________________________________
AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) 
  :TNamed(name,title)
  ,fNMom(0)
  ,fTrackMomentum(0)
  ,fMeanChargeRatio(0)
  ,fNbins(0)
  ,fBinSize(0)
  ,fHistdEdx(0)
  ,fHistTimeBin(0)
{
  //
  //  The main constructor
  //
  
  Init();

}

//_____________________________________________________________________________
AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) 
  :TNamed(c)
  ,fNMom(c.fNMom)
  ,fTrackMomentum(0)
  ,fMeanChargeRatio(c.fMeanChargeRatio)
  ,fNbins(c.fNbins)
  ,fBinSize(c.fBinSize)
  ,fHistdEdx(0)
  ,fHistTimeBin(0)
{
  //
  // Copy constructor
  //

  AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
  
  target.fTrackMomentum = new Double_t[fNMom];
  for (Int_t i=0; i<fNMom; ++i) {
    target.fTrackMomentum[i] = fTrackMomentum[i];
  }
  if (fHistdEdx) {
    target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
  }  

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

}

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

}

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

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

  if (fTrackMomentum) {
    delete[] fTrackMomentum;
    fTrackMomentum = 0;
  }

}

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

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

}

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

  AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
  
  target.CleanUp();
  
  target.fNMom            = fNMom;
  target.fNbins           = fNbins;
  target.fBinSize         = fBinSize;
  target.fMeanChargeRatio = fMeanChargeRatio;
  
  target.fTrackMomentum = new Double_t[fNMom];
  for (Int_t i=0; i<fNMom; ++i) {
    target.fTrackMomentum[i] = fTrackMomentum[i];
  }

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

  TObject::Copy(c);

}

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

  const Int_t kNMom = 11;

  fNMom = kNMom;
  fTrackMomentum = new Double_t[fNMom];
  Double_t trackMomentum[kNMom] = {  0.6,  0.8,  1.0,  1.5,  2.0
                                  ,  3.0,  4.0,  5.0,  6.0,  8.0
                                  , 10.0 };
  for (Int_t imom = 0; imom < kNMom; imom++) {
    fTrackMomentum[imom] = trackMomentum[imom];
  }  

  fHistdEdx    = new TObjArray(AliPID::kSPECIES * fNMom);
  fHistdEdx->SetOwner();
  fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
  fHistTimeBin->SetOwner();  
  
  // ADC Gain normalization
  fMeanChargeRatio = 1.0;
  
  // Number of bins and bin size
  fNbins   = 0;
  fBinSize = 0.0;

}

//_________________________________________________________________________
Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
{
  //
  // Read the TRD dEdx histograms.
  //

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

  // Read histograms
  Char_t text[200];
  for (Int_t particle = 0; particle < AliPID::kSPECIES; ++particle)
  {
    Char_t* particleKey = "";
    switch (particle)
    {
      case AliPID::kElectron: particleKey = "EL"; break;
      case AliPID::kPion: particleKey = "PI"; break;
      case AliPID::kMuon: particleKey = "MU"; break;
      case AliPID::kKaon: particleKey = "KA"; break;
      case AliPID::kProton: particleKey = "PR"; break;
    }
    
    for (Int_t imom = 0; imom < fNMom; imom++) 
    {
      sprintf(text, "h1dEdx%s%01d", particleKey, imom+1);
      TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
      hist->Scale(1.0/hist->Integral());
      fHistdEdx->AddAt(hist, GetHistID(particle, imom));
  
      if (particle == AliPID::kElectron || particle == AliPID::kPion)
      {
        sprintf(text,"h1MaxTimBin%s%01d", particleKey, imom+1);
        TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
        hist->Scale(1.0/hist->Integral());
        fHistTimeBin->AddAt(hist, GetHistID(particle,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  AliTRDCalPIDLQ::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  AliTRDCalPIDLQ::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();

}

//_________________________________________________________________________
TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
{
  //
  // Returns one selected dEdx histogram
  //

  AliInfo(Form("Histogram 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));

}

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

  AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
              ,fpartName[k]
              ,fTrackMomentum[ip]));
  if (k < 0 || k > AliPID::kSPECIES)
    return 0;
  
  return (TH1F*) fHistTimeBin->At(GetHistID(k,ip));

}

//_________________________________________________________________________
Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
{
  //
  // Gets the Probability of having dedx at a given momentum (mom)
  // and particle type k (0 for e) and (2 for pi)
  // from the precalculated de/dx distributions 
  //
  
  Double_t dedx   = dedx1/fMeanChargeRatio;
  Int_t    iEnBin = ((Int_t) (dedx/fBinSize+1));
  if(iEnBin > fNbins) iEnBin = fNbins;

  if (k < 0 || k > AliPID::kSPECIES) {
    return 1;
  }
  
  TH1F* hist1 = 0;
  TH1F* hist2 = 0;
  Double_t mom1 = 0;
  Double_t mom2 = 0;
  
  // Lower limit
  if (mom<=fTrackMomentum[0])  {
    hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,1));
    hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,0));
    mom1 = fTrackMomentum[1];
    mom2 = fTrackMomentum[0];
  }
    
  // Upper Limit
  if(mom>=fTrackMomentum[fNMom-1]) {
    hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-1));
    hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-2));
    mom2 = fTrackMomentum[fNMom-1];
    mom1 = fTrackMomentum[fNMom-2];
  }
    
  // In the range
  for (Int_t ip=1; ip<fNMom; ip++) {
    if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
      hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,ip));
      hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,ip-1));
      mom1 = fTrackMomentum[ip];
      mom2 = fTrackMomentum[ip-1];
    }
  }
  
  Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin)) 
                / (mom1 - mom2);
  return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);

}

//_________________________________________________________________________
Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
{
  //
  // Gets the Probability of having timbin at a given momentum (mom)
  // and particle type k (0 for e) and (2 for pi)
  // from the precalculated timbin distributions 
  //
  
  if (timbin<=0) {
    return 0.0;
  }

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

  if (mom<=fTrackMomentum[0]) {
    return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
  }
  if (mom>=fTrackMomentum[fNMom-1]) { 
    return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
  }

  for (Int_t ip=1; ip<fNMom; ip++) {
    if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
      Double_t slop = (((TH1F*) fHistTimeBin->At(GetHistID(k,ip)))->GetBinContent(iTBin) 
                     - ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin)) 
                    / (fTrackMomentum[ip] - fTrackMomentum[ip-1]);
      // Linear interpolation
      return ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin) 
              + slop * (mom - fTrackMomentum[ip-1]);
    }
  }
  
  return -1;

}
Commit	Line	Data
7754cd1f	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
2745a409	16
	17	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// Container for the distributions of dE/dx and the time bin of the //
	20	// max. cluster for electrons and pions //
	21	// //
	22	// Author: //
	23	// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> //
	24	// //
	25	///////////////////////////////////////////////////////////////////////////////
7754cd1f	26
7754cd1f	27	#include <TH1F.h>
	28	#include <TFile.h>
	29
b92cdef5	30	#include "AliLog.h"
2745a409	31	#include "AliPID.h"
b92cdef5	32
	33	#include "AliTRDCalPIDLQ.h"
	34
7754cd1f	35	ClassImp(AliTRDCalPIDLQ)
	36
	37	Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
	38
	39	//_________________________________________________________________________
2745a409	40	AliTRDCalPIDLQ::AliTRDCalPIDLQ()
	41	:TNamed()
	42	,fNMom(0)
	43	,fTrackMomentum(0)
	44	,fMeanChargeRatio(0)
	45	,fNbins(0)
	46	,fBinSize(0)
	47	,fHistdEdx(0)
	48	,fHistTimeBin(0)
7754cd1f	49	{
	50	//
	51	// The Default constructor
	52	//
	53
	54	}
	55
	56	//_________________________________________________________________________
2745a409	57	AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t name, const Text_t title)
	58	:TNamed(name,title)
	59	,fNMom(0)
	60	,fTrackMomentum(0)
	61	,fMeanChargeRatio(0)
	62	,fNbins(0)
	63	,fBinSize(0)
	64	,fHistdEdx(0)
	65	,fHistTimeBin(0)
7754cd1f	66	{
	67	//
	68	// The main constructor
	69	//
	70
	71	Init();
b92cdef5	72
7754cd1f	73	}
	74
	75	//_____________________________________________________________________________
2745a409	76	AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c)
	77	:TNamed(c)
	78	,fNMom(c.fNMom)
	79	,fTrackMomentum(0)
	80	,fMeanChargeRatio(c.fMeanChargeRatio)
	81	,fNbins(c.fNbins)
	82	,fBinSize(c.fBinSize)
	83	,fHistdEdx(0)
	84	,fHistTimeBin(0)
7754cd1f	85	{
7754cd1f	86	//
b92cdef5	87	// Copy constructor
7754cd1f	88	//
7754cd1f	89
2745a409	90	AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
7754cd1f	91
2745a409	92	target.fTrackMomentum = new Double_t[fNMom];
	93	for (Int_t i=0; i<fNMom; ++i) {
	94	target.fTrackMomentum[i] = fTrackMomentum[i];
	95	}
	96	if (fHistdEdx) {
	97	target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
	98	}
	99
	100	if (fHistTimeBin) {
	101	target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
	102	}
b92cdef5	103
7754cd1f	104	}
	105
	106	//_________________________________________________________________________
	107	AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
	108	{
	109	//
	110	// Destructor
	111	//
	112
	113	CleanUp();
b92cdef5	114
7754cd1f	115	}
	116
	117	//_________________________________________________________________________
	118	void AliTRDCalPIDLQ::CleanUp()
	119	{
b92cdef5	120	//
	121	// Delets all newly created objects
	122	//
	123
2745a409	124	if (fHistdEdx) {
7754cd1f	125	delete fHistdEdx;
	126	fHistdEdx = 0;
	127	}
	128
2745a409	129	if (fHistTimeBin) {
7754cd1f	130	delete fHistTimeBin;
	131	fHistTimeBin = 0;
	132	}
	133
2745a409	134	if (fTrackMomentum) {
7754cd1f	135	delete[] fTrackMomentum;
	136	fTrackMomentum = 0;
	137	}
b92cdef5	138
7754cd1f	139	}
	140
	141	//_____________________________________________________________________________
	142	AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
	143	{
	144	//
	145	// Assignment operator
	146	//
	147
	148	if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
	149	return *this;
b92cdef5	150
7754cd1f	151	}
	152
	153	//_____________________________________________________________________________
	154	void AliTRDCalPIDLQ::Copy(TObject &c) const
	155	{
	156	//
	157	// Copy function
	158	//
	159
	160	AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
	161
	162	target.CleanUp();
	163
2745a409	164	target.fNMom = fNMom;
	165	target.fNbins = fNbins;
	166	target.fBinSize = fBinSize;
	167	target.fMeanChargeRatio = fMeanChargeRatio;
7754cd1f	168
7754cd1f	169	target.fTrackMomentum = new Double_t[fNMom];
2745a409	170	for (Int_t i=0; i<fNMom; ++i) {
7754cd1f	171	target.fTrackMomentum[i] = fTrackMomentum[i];
2745a409	172	}
7754cd1f	173
2745a409	174	if (fHistdEdx) {
7754cd1f	175	target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
2745a409	176	}
2745a409	177	if (fHistTimeBin) {
7754cd1f	178	target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
2745a409	179	}
2745a409	180
7754cd1f	181	TObject::Copy(c);
b92cdef5	182
7754cd1f	183	}
	184
	185	//_________________________________________________________________________
	186	void AliTRDCalPIDLQ::Init()
	187	{
b92cdef5	188	//
	189	// Initialization
	190	//
	191
	192	const Int_t kNMom = 11;
	193
	194	fNMom = kNMom;
7754cd1f	195	fTrackMomentum = new Double_t[fNMom];
b92cdef5	196	Double_t trackMomentum[kNMom] = { 0.6, 0.8, 1.0, 1.5, 2.0
	197	, 3.0, 4.0, 5.0, 6.0, 8.0
	198	, 10.0 };
	199	for (Int_t imom = 0; imom < kNMom; imom++) {
7754cd1f	200	fTrackMomentum[imom] = trackMomentum[imom];
b92cdef5	201	}
	202
	203	fHistdEdx = new TObjArray(AliPID::kSPECIES * fNMom);
7754cd1f	204	fHistdEdx->SetOwner();
	205	fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
	206	fHistTimeBin->SetOwner();
	207
	208	// ADC Gain normalization
b92cdef5	209	fMeanChargeRatio = 1.0;
7754cd1f	210
7754cd1f	211	// Number of bins and bin size
b92cdef5	212	fNbins = 0;
	213	fBinSize = 0.0;
	214
7754cd1f	215	}
	216
	217	//_________________________________________________________________________
	218	Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
	219	{
	220	//
	221	// Read the TRD dEdx histograms.
	222	//
b92cdef5	223
7754cd1f	224	// Read histogram Root file
	225	TFile *histFile = new TFile(responseFile, "READ");
	226	if (!histFile \|\| !histFile->IsOpen()) {
2745a409	227	AliError(Form("Opening TRD histgram file %s failed", responseFile));
7754cd1f	228	return kFALSE;
	229	}
	230	gROOT->cd();
	231
	232	// Read histograms
	233	Char_t text[200];
	234	for (Int_t particle = 0; particle < AliPID::kSPECIES; ++particle)
	235	{
	236	Char_t* particleKey = "";
	237	switch (particle)
	238	{
	239	case AliPID::kElectron: particleKey = "EL"; break;
	240	case AliPID::kPion: particleKey = "PI"; break;
	241	case AliPID::kMuon: particleKey = "MU"; break;
	242	case AliPID::kKaon: particleKey = "KA"; break;
	243	case AliPID::kProton: particleKey = "PR"; break;
	244	}
	245
	246	for (Int_t imom = 0; imom < fNMom; imom++)
	247	{
	248	sprintf(text, "h1dEdx%s%01d", particleKey, imom+1);
	249	TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
	250	hist->Scale(1.0/hist->Integral());
	251	fHistdEdx->AddAt(hist, GetHistID(particle, imom));
	252
	253	if (particle == AliPID::kElectron \|\| particle == AliPID::kPion)
	254	{
	255	sprintf(text,"h1MaxTimBin%s%01d", particleKey, imom+1);
	256	TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
	257	hist->Scale(1.0/hist->Integral());
	258	fHistTimeBin->AddAt(hist, GetHistID(particle,imom));
	259	}
	260	}
	261	}
	262
	263	histFile->Close();
	264	delete histFile;
	265
	266	// Number of bins and bin size
	267	TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
	268	fNbins = hist->GetNbinsX();
	269	fBinSize = hist->GetBinWidth(1);
	270
	271	return kTRUE;
b92cdef5	272
7754cd1f	273	}
	274
	275	//_________________________________________________________________________
	276	Double_t AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
	277	{
	278	//
	279	// Gets mean of de/dx dist. of e
b92cdef5	280	//
b92cdef5	281
2745a409	282	AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
	283	,fpartName[k]
	284	,fTrackMomentum[ip]));
	285	if (k < 0 \|\| k > AliPID::kSPECIES) {
7754cd1f	286	return 0;
2745a409	287	}
2745a409	288
7754cd1f	289	return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
b92cdef5	290
7754cd1f	291	}
	292
	293	//_________________________________________________________________________
	294	Double_t AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
	295	{
	296	//
	297	// Gets Normalization of de/dx dist. of e
b92cdef5	298	//
7754cd1f	299
2745a409	300	AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
	301	,fpartName[k]
	302	,fTrackMomentum[ip]));
	303	if (k < 0 \|\| k > AliPID::kSPECIES) {
7754cd1f	304	return 0;
2745a409	305	}
7754cd1f	306
7754cd1f	307	return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
b92cdef5	308
7754cd1f	309	}
	310
	311	//_________________________________________________________________________
	312	TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
	313	{
	314	//
b92cdef5	315	// Returns one selected dEdx histogram
7754cd1f	316	//
b92cdef5	317
2745a409	318	AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
	319	,fpartName[k]
	320	,fTrackMomentum[ip]));
	321	if (k < 0 \|\| k > AliPID::kSPECIES) {
7754cd1f	322	return 0;
2745a409	323	}
7754cd1f	324
7754cd1f	325	return (TH1F*) fHistdEdx->At(GetHistID(k,ip));
b92cdef5	326
	327	}
	328
	329	//_________________________________________________________________________
	330	TH1F* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
	331	{
	332	//
	333	// Returns one selected time bin max histogram
	334	//
	335
2745a409	336	AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
	337	,fpartName[k]
	338	,fTrackMomentum[ip]));
b92cdef5	339	if (k < 0 \|\| k > AliPID::kSPECIES)
	340	return 0;
	341
	342	return (TH1F*) fHistTimeBin->At(GetHistID(k,ip));
	343
7754cd1f	344	}
	345
	346	//_________________________________________________________________________
	347	Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
	348	{
	349	//
	350	// Gets the Probability of having dedx at a given momentum (mom)
	351	// and particle type k (0 for e) and (2 for pi)
	352	// from the precalculated de/dx distributions
b92cdef5	353	//
7754cd1f	354
2745a409	355	Double_t dedx = dedx1/fMeanChargeRatio;
2745a409	356	Int_t iEnBin = ((Int_t) (dedx/fBinSize+1));
7754cd1f	357	if(iEnBin > fNbins) iEnBin = fNbins;
7754cd1f	358
2745a409	359	if (k < 0 \|\| k > AliPID::kSPECIES) {
7754cd1f	360	return 1;
2745a409	361	}
7754cd1f	362
	363	TH1F* hist1 = 0;
	364	TH1F* hist2 = 0;
	365	Double_t mom1 = 0;
	366	Double_t mom2 = 0;
	367
	368	// Lower limit
2745a409	369	if (mom<=fTrackMomentum[0]) {
7754cd1f	370	hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,1));
	371	hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,0));
	372	mom1 = fTrackMomentum[1];
	373	mom2 = fTrackMomentum[0];
	374	}
	375
	376	// Upper Limit
2745a409	377	if(mom>=fTrackMomentum[fNMom-1]) {
7754cd1f	378	hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-1));
	379	hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-2));
	380	mom2 = fTrackMomentum[fNMom-1];
	381	mom1 = fTrackMomentum[fNMom-2];
	382	}
	383
	384	// In the range
2745a409	385	for (Int_t ip=1; ip<fNMom; ip++) {
2745a409	386	if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
7754cd1f	387	hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,ip));
	388	hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,ip-1));
	389	mom1 = fTrackMomentum[ip];
	390	mom2 = fTrackMomentum[ip-1];
	391	}
	392	}
	393
b92cdef5	394	Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin))
b92cdef5	395	/ (mom1 - mom2);
7754cd1f	396	return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);
b92cdef5	397
7754cd1f	398	}
	399
	400	//_________________________________________________________________________
	401	Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
	402	{
	403	//
	404	// Gets the Probability of having timbin at a given momentum (mom)
	405	// and particle type k (0 for e) and (2 for pi)
	406	// from the precalculated timbin distributions
b92cdef5	407	//
7754cd1f	408
2745a409	409	if (timbin<=0) {
	410	return 0.0;
	411	}
	412
	413	Int_t iTBin = timbin+1;
7754cd1f	414
2745a409	415	// Everything which is not an electron counts as a pion for time bin max
2745a409	416	if (k != AliPID::kElectron) {
7754cd1f	417	k = AliPID::kPion;
2745a409	418	}
7754cd1f	419
2745a409	420	if (mom<=fTrackMomentum[0]) {
7754cd1f	421	return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
2745a409	422	}
2745a409	423	if (mom>=fTrackMomentum[fNMom-1]) {
7754cd1f	424	return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
2745a409	425	}
	426
	427	for (Int_t ip=1; ip<fNMom; ip++) {
	428	if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
b92cdef5	429	Double_t slop = (((TH1F*) fHistTimeBin->At(GetHistID(k,ip)))->GetBinContent(iTBin)
	430	- ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin))
	431	/ (fTrackMomentum[ip] - fTrackMomentum[ip-1]);
2745a409	432	// Linear interpolation
b92cdef5	433	return ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin)
b92cdef5	434	+ slop * (mom - fTrackMomentum[ip-1]);
7754cd1f	435	}
	436	}
	437
	438	return -1;
b92cdef5	439
7754cd1f	440	}