* provided "as is" without express or implied warranty. *
**************************************************************************/
-//-----------------------------------------------------------------
-// Class for dE/dx and Time Bin of Max. Cluster for Electrons and
-// pions in TRD.
-// It is instantiated in class AliTRDpidESD for particle identification
-// in TRD
-// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>
-//-----------------------------------------------------------------
+/* $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 "AliTRDCalPIDLQ.h"
#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 "AliTRDCalPIDLQ.h"
+#include "AliTRDcalibDB.h"
ClassImp(AliTRDCalPIDLQ)
Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
+
+Char_t* AliTRDCalPIDLQ::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
+
+Float_t AliTRDCalPIDLQ::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 AliTRDCalPIDLQ::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};
+
//_________________________________________________________________________
-AliTRDCalPIDLQ::AliTRDCalPIDLQ():
- TNamed(),
- fNMom(0),
- fTrackMomentum(0),
- fMeanChargeRatio(0),
- fNbins(0),
- fBinSize(0),
- fHistdEdx(0),
- fHistTimeBin(0)
+AliTRDCalPIDLQ::AliTRDCalPIDLQ()
+ :TNamed("pid", "PID for TRD")
+ ,fMeanChargeRatio(0)
+ ,fHistdEdx(0x0)
+ ,fHistTimeBin(0x0)
{
//
// The Default constructor
//
-
+
+ Init();
+
}
//_________________________________________________________________________
-AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) : TNamed(name, title)
+AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
+ :TNamed(name,title)
+ ,fMeanChargeRatio(0)
+ ,fHistdEdx(0x0)
+ ,fHistTimeBin(0x0)
{
//
// The main constructor
//
Init();
+
}
//_____________________________________________________________________________
-AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) : TNamed(c)
+AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c)
+ :TNamed(c)
+ ,fMeanChargeRatio(c.fMeanChargeRatio)
+ ,fHistdEdx(0x0)
+ ,fHistTimeBin(0x0)
{
//
- // copy constructor
+ // Copy constructor
//
- Init();
+ if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
- ((AliTRDCalPIDLQ &) c).Copy(*this);
}
//_________________________________________________________________________
//
CleanUp();
+
}
//_________________________________________________________________________
void AliTRDCalPIDLQ::CleanUp()
{
- if (fHistdEdx)
- {
+ //
+ // Delets all newly created objects
+ //
+
+ if (fHistdEdx) {
delete fHistdEdx;
- fHistdEdx = 0;
+ fHistdEdx = 0x0;
}
- if (fHistTimeBin)
- {
+ if (fHistTimeBin) {
delete fHistTimeBin;
- fHistTimeBin = 0;
- }
-
- if (fTrackMomentum)
- {
- delete[] fTrackMomentum;
- fTrackMomentum = 0;
+ fHistTimeBin = 0x0;
}
}
if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
return *this;
+
}
//_____________________________________________________________________________
target.CleanUp();
- target.fNMom = fNMom;
-
- target.fTrackMomentum = new Double_t[fNMom];
- for (Int_t i=0; i<fNMom; ++i)
- target.fTrackMomentum[i] = fTrackMomentum[i];
-
target.fMeanChargeRatio = fMeanChargeRatio;
- target.fNbins = fNbins;
- target.fBinSize = fBinSize;
-
- if (fHistdEdx)
+ if (fHistdEdx) {
target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
-
- if (fHistTimeBin)
+ }
+ if (fHistTimeBin) {
target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
-
+ }
+
TObject::Copy(c);
+
}
//_________________________________________________________________________
void AliTRDCalPIDLQ::Init()
{
- fNMom = 11;
-
- fTrackMomentum = new Double_t[fNMom];
- Double_t trackMomentum[11] = {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10};
- for (Int_t imom = 0; imom < 11; imom++)
- fTrackMomentum[imom] = trackMomentum[imom];
-
- fHistdEdx = new TObjArray(AliPID::kSPECIES * fNMom);
+ //
+ // Initialization
+ //
+
+ fHistdEdx = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
fHistdEdx->SetOwner();
- fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
+ 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;
-
- // Number of bins and bin size
- fNbins = 0;
- fBinSize = 0.;
+ fMeanChargeRatio = 1.0;
}
//_________________________________________________________________________
-Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
+Bool_t AliTRDCalPIDLQ::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 = new TFile(responseFile, "READ");
+ TFile *histFile = TFile::Open(refFile, "READ");
if (!histFile || !histFile->IsOpen()) {
- Error("AliTRDCalPIDLQ", "opening TRD histgram file %s failed", responseFile);
+ AliError(Form("Opening TRD histgram file %s failed", refFile));
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));
- }
- }
+ 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);
+ //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
- printf("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::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();
+//
+// }
+
//_________________________________________________________________________
-Double_t AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
+TH1* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
{
//
- // Gets Normalization of de/dx dist. of e
+ // Returns one selected dEdx histogram
+ //
+
+ if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
+ if(ip<0 || ip>= kNMom ) return 0x0;
- printf("Normalization for particle = %s and momentum = %.2f is:\n",fpartName[k], fTrackMomentum[ip]);
- if (k < 0 || k > AliPID::kSPECIES)
- return 0;
+ AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
- return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
+ return (TH1*)fHistdEdx->At(GetHistID(k, ip));
+
}
//_________________________________________________________________________
-TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
+TH1* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
{
//
+ // Returns one selected time bin max histogram
//
- printf("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));
+
+ 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 AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) 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;
+ // Core function of AliTRDCalPIDLQ 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
+ return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
- 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
+Double_t AliTRDCalPIDLQ::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 k (0 for e) and (2 for pi)
+ // and particle type (spec) (0 for e) and (2 for pi)
// from the precalculated timbin distributions
+ //
- if (timbin<=0)
- return 0.;
- Int_t iTBin=timbin+1;
-
- // everything which is not electron counts as pion for time bin
- if (k != AliPID::kElectron)
- k = AliPID::kPion;
+ if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
- if (mom<=fTrackMomentum[0])
- return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
+ Int_t iTBin = timbin+1;
- if (mom>=fTrackMomentum[fNMom-1])
- return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
+ // Everything which is not an electron counts as a pion for time bin max
+ //if(spec != AliPID::kElectron) spec = AliPID::kPion;
- 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;
+ 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);
}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff