/**************************************************************************
- * 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. *
- **************************************************************************/
+* 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 //
-// //
-// Author: //
-// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> //
-// //
-///////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+// //
+// 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 <TFile.h>
#include <TROOT.h>
+#include <TSystem.h>
#include "AliLog.h"
#include "AliPID.h"
+#include "TKDPDF.h"
+
#include "AliTRDCalPIDLQ.h"
+#include "AliTRDcalibDB.h"
ClassImp(AliTRDCalPIDLQ)
-Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
-
+Float_t AliTRDCalPIDLQ::fgTrackSegLength[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)
+ :AliTRDCalPID("pid", "LQ PID references for TRD")
{
//
// 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();
- }
+ //Init();
}
//_________________________________________________________________________
-AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
+AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
+ :AliTRDCalPID(name,title)
{
//
- // 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
+ // The main constructor
//
-
- 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;
+ Init();
}
//_________________________________________________________________________
-Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
+Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
{
//
// 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));
+ if(gSystem->AccessPathName(refFile, kReadPermission)){
+ AliError(Form("File %s.root not readable", 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));
- }
- }
- }
-
- 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;
+ if(!TFile::Open(refFile)){
+ AliError(Form("File %s corrupted", refFile));
+ return kFALSE;
}
-
- 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;
+ TObjArray *pdf(NULL);
+ if (!( pdf = dynamic_cast<TObjArray*>(gFile->Get("PDF_LQ")))) {
+ AliError("PID data not available");
+ return kFALSE;
}
-
- return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
-
+ fModel=(TObjArray*)pdf->Clone("PDF");
+ gFile->Close();
+ return kTRUE;
}
//_________________________________________________________________________
-TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
+TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t is, Int_t slices) 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));
+ if (is < 0 || is >= AliPID::kSPECIES) return NULL;
+ if(ip<0 || ip>= kNMom ) return NULL;
+ AliDebug(2, Form("Retrive dEdx distribution for %s @ p=%5.2f [GeV/c].", AliPID::ParticleShortName(is), fgTrackMomentum[ip]));
+ return fModel->At(GetModelID(ip, is, slices));
}
//_________________________________________________________________________
-TH1F* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
+Int_t AliTRDCalPIDLQ::GetNrefs()
{
- //
- // 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));
-
+// Returns the number of PDF distribution
+ return AliTRDCalPID::kNMom*AliPID::kSPECIES*2;
}
//_________________________________________________________________________
-Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
+ , const Float_t * const dedx
+ , Float_t length
+ , Int_t slices) 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];
- }
+//
+// 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.;
+
+ Bool_t k2D(slices==2);
+ Double_t x[]={0., 0.};
+ if(!CookdEdx(dedx, x, k2D)) return 0.;
- // 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];
- }
+ // find the interval in momentum and track segment length which applies for this data
+/* Int_t ilength = 1;
+ while(ilength<kNLength-1 && length>fgTrackSegLength[ilength]){
+ ilength++;
+ }*/
+ Int_t imom = 1;
+ while(imom<kNMom-1 && mom>fgTrackMomentum[imom]) imom++;
+
+ Double_t p[2], e[2], r;
+ TKDPDF *pdf(NULL);
+
+ AliDebug(2, Form("Looking %cD for %s@p=%6.4f[GeV/c] log(dEdx)={%7.2f %7.2f}[a.u.] l=%4.2f[cm].", k2D?'2':'1', AliPID::ParticleShortName(spec), mom, x[0], x[1], length));
+ if(!(pdf = dynamic_cast<TKDPDF*>(fModel->At(GetModelID(imom-1, spec, slices))))) {
+ AliError(Form("%cD Ref data @ idx[%d] not found in DB.", k2D?'2':'1', GetModelID(imom-1, spec, slices)));
+ return 0.;
+ }
+ if(!pdf->GetSize()) pdf->Bootstrap();
+ pdf->Eval(x, r, e[0], kFALSE);
+ p[0]=TMath::Abs(r); // conversion from interpolation to PDF
+ AliDebug(2, Form("LQ=%6.3f+-%5.2f%% @ %4.1f[GeV/c]", p[0], 1.E2*e[0]/p[0], fgTrackMomentum[imom-1]));
+
+ if(!(pdf = dynamic_cast<TKDPDF*>(fModel->At(GetModelID(imom, spec, slices))))){
+ AliError(Form("%cD Ref data @ idx[%d] not found in DB.", k2D?'2':'1', GetModelID(imom, spec, slices)));
+ return p[0];
+ }
+ if(!pdf->GetSize()) pdf->Bootstrap();
+ pdf->Eval(x, r, e[1], kFALSE);
+ p[1]=TMath::Abs(r); // conversion from interpolation to PDF
+ AliDebug(2, Form("LQ=%6.3f+-%5.2f%% @ %4.1f[GeV/c]", p[1], 1.E2*e[1]/p[1], fgTrackMomentum[imom]));
+
+ // return interpolation over momentum binning
+ if(mom < fgTrackMomentum[0]) return p[0];
+ else if(mom > fgTrackMomentum[kNMom-1]) return p[1];
+ else{
+ Double_t lmom[2] = {fgTrackMomentum[imom-1], fgTrackMomentum[imom]};
+ return p[0] + (p[1] - p[0])*(mom - lmom[0])/(lmom[1] - lmom[0]);
}
-
- 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
+void AliTRDCalPIDLQ::Init()
{
- //
- // 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;
-
+//
+// PID LQ list initialization
+//
+ fModel = new TObjArray(AliPID::kSPECIES * kNMom * 2);
+ fModel->SetOwner();
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff