X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FCal%2FAliTRDCalPID.cxx;h=383ea91e989b22f111391959dd3881e2c556e76a;hb=f2391aef7f4e428e69d4be5acb9a09238276cd5b;hp=3541e0a879d45e7fdedd2aa0b26824bb7e6efd85;hpb=720a0a16fa8d7f72fbb392f133f2999631f9dad4;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/Cal/AliTRDCalPID.cxx b/TRD/Cal/AliTRDCalPID.cxx index 3541e0a879d..383ea91e989 100644 --- a/TRD/Cal/AliTRDCalPID.cxx +++ b/TRD/Cal/AliTRDCalPID.cxx @@ -17,8 +17,7 @@ ////////////////////////////////////////////////////////////////////////// // // -// Container for the distributions of dE/dx and the time bin of the // -// max. cluster for electrons and pions // +// Container for PID information // // // // Authors: // // Prashant Shukla // @@ -26,75 +25,42 @@ // // ////////////////////////////////////////////////////////////////////////// -#include -#include -#include -#include - -#include "AliLog.h" -#include "AliPID.h" -#include "AliESD.h" -#include "AliESDtrack.h" +#include #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}; +const Char_t* AliTRDCalPID::fPartName[AliPID::kSPECIES] = { "electron", "muon", "pion", "kaon", "proton"}; +const Char_t* AliTRDCalPID::fPartSymb[AliPID::kSPECIES] = { "EL", "MU", "PI", "KA", "PR"}; +Color_t AliTRDCalPID::fgPartColor[AliPID::kSPECIES] = { kRed, kGreen, kBlue, kYellow, kMagenta}; +Float_t AliTRDCalPID::fgTrackMomentum[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::fgTrackMomentumBinning[kNMom+1] = { 0.5, 0.7, 0.9, 1.25, 1.75, 2.5 + , 3.5, 4.5, 5.5, 7.0, 9.0, 12.0 }; - //_________________________________________________________________________ AliTRDCalPID::AliTRDCalPID() :TNamed("pid", "PID for TRD") - ,fMeanChargeRatio(0) - ,fHistdEdx(0x0) - ,fHistTimeBin(0x0) + ,fModel(0x0) { // // The Default constructor // - Init(); - } -//_________________________________________________________________________ +//_____________________________________________________________________________ AliTRDCalPID::AliTRDCalPID(const Text_t *name, const Text_t *title) :TNamed(name,title) - ,fMeanChargeRatio(0) - ,fHistdEdx(0x0) - ,fHistTimeBin(0x0) + ,fModel(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() { @@ -102,311 +68,21 @@ AliTRDCalPID::~AliTRDCalPID() // Destructor // - CleanUp(); - -} - -//_________________________________________________________________________ -void AliTRDCalPID::CleanUp() -{ - // - // Delets all newly created objects - // - - if (fHistdEdx) { - delete fHistdEdx; - fHistdEdx = 0x0; - } - - if (fHistTimeBin) { - delete fHistTimeBin; - fHistTimeBin = 0x0; + if (fModel) { + delete fModel; } -} -//_____________________________________________________________________________ -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) +Int_t AliTRDCalPID::GetPartIndex(Int_t pdg) { // - // Read the TRD dEdx histograms. + // Return the index to a given particle, defined by its PDG code // - 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; + for(Int_t is=0; iscd(); - - // 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(ilengthfTrackSegLength[ilength]){ - ilength++; - } - Int_t imom = 1; - while(imomfTrackMomentum[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); - -} - -//_________________________________________________________________________ -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(imomfTrackMomentum[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); + return -1; } -