#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
-#include <TTree.h>
#include <TROOT.h>
-#include <TPDGCode.h>
-#include <TParticle.h>
-#include <TPrincipal.h>
#include "AliLog.h"
-#include "AliHeader.h"
-#include "AliGenEventHeader.h"
-#include "AliRun.h"
-#include "AliRunLoader.h"
-#include "AliStack.h"
#include "AliPID.h"
-#include "AliESD.h"
-#include "AliESDtrack.h"
#include "AliTRDCalPIDLQ.h"
-#include "AliTRDCalPIDLQRef.h"
-#include "AliTRDpidESD.h"
#include "AliTRDcalibDB.h"
-#include "AliTRDtrack.h"
-#include "AliTRDgeometry.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::fTrackSegLength[kNLength] = { 3.7, 3.9, 4.2, 5.0 };
+
//_________________________________________________________________________
AliTRDCalPIDLQ::AliTRDCalPIDLQ()
- :TNamed("pid", "PID for TRD")
- ,fNMom(0)
- ,fTrackMomentum(0x0)
- ,fNLength(0)
- ,fTrackSegLength(0x0)
- ,fNTimeBins(0)
- ,fMeanChargeRatio(0)
- ,fNbins(0)
- ,fBinSize(0)
- ,fHistdEdx(0x0)
- ,fHistTimeBin(0x0)
- ,fEstimator(0x0)
+ :AliTRDCalPID("pid", "LQ PID references for TRD")
{
//
// The Default constructor
//
- Init();
+ Init();
+
}
//_________________________________________________________________________
-AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
- :TNamed(name,title)
- ,fNMom(0)
- ,fTrackMomentum(0x0)
- ,fNLength(0)
- ,fTrackSegLength(0x0)
- ,fNTimeBins(0)
- ,fMeanChargeRatio(0)
- ,fNbins(0)
- ,fBinSize(0)
- ,fHistdEdx(0x0)
- ,fHistTimeBin(0x0)
- ,fEstimator(0x0)
+AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
+ :AliTRDCalPID(name,title)
{
//
// The main constructor
}
-//_____________________________________________________________________________
-AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c)
- :TNamed(c)
- ,fNMom(c.fNMom)
- ,fTrackMomentum(0x0)
- ,fNLength(c.fNLength)
- ,fTrackSegLength(0x0)
- ,fNTimeBins(c.fNTimeBins)
- ,fMeanChargeRatio(c.fMeanChargeRatio)
- ,fNbins(c.fNbins)
- ,fBinSize(c.fBinSize)
- ,fHistdEdx(0x0)
- ,fHistTimeBin(0x0)
- ,fEstimator(0x0)
-{
- //
- // Copy constructor
- //
-
- if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
-
-}
-
//_________________________________________________________________________
AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
{
// Destructor
//
- CleanUp();
-
-}
-
-//_________________________________________________________________________
-void AliTRDCalPIDLQ::CleanUp()
-{
- //
- // Delets all newly created objects
- //
-
- if (fHistdEdx) {
- delete fHistdEdx;
- fHistdEdx = 0x0;
- }
-
- if (fHistTimeBin) {
- delete fHistTimeBin;
- fHistTimeBin = 0x0;
- }
-
- if (fTrackMomentum) {
- delete[] fTrackMomentum;
- fTrackMomentum = 0x0;
- }
-
- if (fTrackSegLength) {
- delete[] fTrackSegLength;
- fTrackSegLength = 0x0;
- }
-
- if (fEstimator) {
- delete fEstimator;
- fEstimator = 0x0;
- }
-
-}
-
-//_____________________________________________________________________________
-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.fNbins = fNbins;
- target.fBinSize = fBinSize;
- target.fMeanChargeRatio = fMeanChargeRatio;
- target.fNTimeBins = fNTimeBins;
-
- //target.fNMom = fNMom;
- target.fTrackMomentum = new Double_t[fNMom];
- for (Int_t i=0; i<fNMom; ++i) {
- target.fTrackMomentum[i] = fTrackMomentum[i];
- }
-
- //target.fNLength = fNLength;
- target.fTrackSegLength = new Double_t[fNLength];
- for (Int_t i=0; i<fNLength; ++i) {
- target.fTrackSegLength[i] = fTrackSegLength[i];
- }
-
- if (fHistdEdx) {
- target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
- }
- if (fHistTimeBin) {
- target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
- }
-
- target.fEstimator = new AliTRDCalPIDLQRef(*fEstimator);
-
- TObject::Copy(c);
-
}
//_________________________________________________________________________
-void AliTRDCalPIDLQ::Init()
+Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
{
//
- // Initialization
+ // Read the TRD dEdx histograms.
//
- fNMom = 11;
- fNLength = 4;
-
- fTrackMomentum = new Double_t[fNMom];
- fTrackMomentum[0] = 0.6;
- fTrackMomentum[1] = 0.8;
- fTrackMomentum[2] = 1.0;
- fTrackMomentum[3] = 1.5;
- fTrackMomentum[4] = 2.0;
- fTrackMomentum[5] = 3.0;
- fTrackMomentum[6] = 4.0;
- fTrackMomentum[7] = 5.0;
- fTrackMomentum[8] = 6.0;
- fTrackMomentum[9] = 8.0;
- fTrackMomentum[10] = 10.0;
-
- fTrackSegLength = new Double_t[fNLength];
- fTrackSegLength[0] = 3.7;
- fTrackSegLength[1] = 3.9;
- fTrackSegLength[2] = 4.2;
- fTrackSegLength[3] = 5.0;
-
- fHistdEdx = new TObjArray(AliPID::kSPECIES * fNMom/* * fNLength*/);
- fHistdEdx->SetOwner();
- fHistTimeBin = new TObjArray(2 * fNMom);
- 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 AliTRDCalPIDLQRef();
-
- // Number of Time bins
+ Int_t nTimeBins = 22;
+ // Get number of time bins from CDB
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if(!calibration){
- AliWarning("No AliTRDcalibDB available. Using 22 time bins.");
- fNTimeBins = 22;
- } else fNTimeBins = calibration->GetNumberOfTimeBins();
-
- // ADC Gain normalization
- fMeanChargeRatio = 1.0;
-
- // Number of bins and bin size
- fNbins = 0;
- fBinSize = 0.0;
-}
-
-//_________________________________________________________________________
-Bool_t AliTRDCalPIDLQ::ReadReferences(Char_t *responseFile)
-{
- //
- // Read the TRD dEdx histograms.
- //
+ 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()) {
- AliError(Form("Opening TRD histgram file %s failed", responseFile));
+ 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 < fNMom; imom++){
- TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
+ for (Int_t imom = 0; imom < kNMom; imom++){
+ TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
hist->Scale(1./hist->Integral());
- fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
+ fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
- if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
-
- TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
- if(ht->GetNbinsX() != fNTimeBins) 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, fNTimeBins));
- ht->Scale(1./ht->Integral());
- fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*fNMom + 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);
}
}
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();
-//
-// }
+}
//_________________________________________________________________________
-TH1* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
+TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const // iType not needed
{
//
// Returns one selected dEdx histogram
//
- if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
- if(ip<0 || ip>= fNMom ) return 0x0;
+ if (iType < 0 || iType >= 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]));
+ AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fTrackMomentum[ip]));
- return (TH1*)fHistdEdx->At(GetHistID(k, ip));
-
+ return fModel->At(GetModelID(ip, iType, iplane));
}
//_________________________________________________________________________
-TH1* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length, Int_t iplane) const
{
//
- // Returns one selected time bin max histogram
- //
-
- if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
- if(ip<0 || ip>= fNMom ) 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)*fNMom+ip);
-}
-
-//_________________________________________________________________________
-Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx, Double_t length) const
-{
- //
- // Core function of AliTRDCalPIDLQ class for calculating the
+ // 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
// find the interval in momentum and track segment length which applies for this data
Int_t ilength = 1;
- while(ilength<fNLength-1 && length>fTrackSegLength[ilength]){
+ while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
ilength++;
}
Int_t imom = 1;
- while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
-
+ 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*/)))){
+ if(!(hist = (TH2D*)fModel->At(GetModelID(imom-1, spec, iplane)))){
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)));
+ AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
return 0.;
}
ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
- Bool_t kX = (dedx[0] < ax->GetBinUpEdge(nbinsx));
- Bool_t kY = (dedx[1] < ay->GetBinUpEdge(nbinsy));
+ 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(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
- else LQ1 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
+ 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(dedx[1]));
+ if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
- if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
+ if(!(hist = (TH2D*)fModel->At(GetModelID(imom, spec, iplane)))){
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)));
+ AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom, spec, iplane)));
return LQ1;
}
if(kX)
- if(kY) LQ2 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
- else LQ2 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
+ 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(dedx[1]));
+ 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);
+ // 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 AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
+void AliTRDCalPIDLQ::Init()
{
//
- // 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
+ // Initialization
//
-
- if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
- if (timbin<0 || timbin >= fNTimeBins) 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<fNMom-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)*fNMom+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)*fNMom+imom-1));
- return 0.;
- }
- Double_t LQ1 = hist->GetBinContent(iTBin);
-
- if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+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)*fNMom+imom));
- return LQ1;
- }
- Double_t LQ2 = hist->GetBinContent(iTBin);
-
- // return interpolation over momentum binning
- return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
-}
-
-//__________________________________________________________________
-Bool_t AliTRDCalPIDLQ::WriteReferences(Char_t *File, Char_t *dir)
-{
- // Build, Fill and write to file the histograms used for PID.
- // The simulations are looked in the
- // directories with the general form Form("p%3.1f", momentum)
- // starting from dir (default .). Here momentum belongs to the list
- // of known momentum to PID (see fTrackMomentum).
- // The output histograms are
- // written to the file "File" in cwd (default
- // TRDPIDHistograms.root). In order to build a DB entry
- // consider running $ALICE_ROOT/Cal/AliTRDCreateDummyCDB.C
- //
- // Author:
- // Alex Bercuci (A.Bercuci@gsi.de)
-
- const Int_t nDirs = 1;
- Int_t partCode[AliPID::kSPECIES] =
- {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
-
- // minimal test of simulation location
- TFile *f = new TFile(Form("%s/p%3.1f/galice.root", dir, 2.));
- if(!f || f->IsZombie()){
- AliError(Form("Could not access file galice in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
- return kFALSE;
- }
- f->Close(); delete f;
- f = new TFile(Form("%s/p%3.1f/AliESDs.root", dir, 2.));
- if(!f || f->IsZombie()){
- AliError(Form("Could not access file AliESDs in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
- return kFALSE;
- }
- f->Close(); delete f;
-
- // Init statistics
- Int_t nPart[AliPID::kSPECIES], nTotPart;
- printf("P[GeV/c] ");
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", fpartSymb[ispec]);
- printf("\n-----------------------------------------------\n");
-
-
+ fModel = new TObjArray(AliPID::kSPECIES * kNMom);
+ fModel -> SetOwner();
- // Build PID reference histograms and reference object
- const Int_t color[] = {4, 3, 2, 7, 6};
- for (Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
- if (ispec != AliPID::kElectron && ispec != AliPID::kPion) continue;
-
- h1MaxTB[(ispec>0)?1:0] = new TH1F(Form("h1%s", fpartSymb[ispec]), "", fNTimeBins, -.5, fNTimeBins-.5);
- h1MaxTB[(ispec>0)?1:0]->SetLineColor(color[ispec]);
- }
- AliTRDCalPIDLQRef ref;
-
- // momentum loop
- AliRunLoader *fRunLoader = 0x0;
- TFile *esdFile = 0x0;
- TTree *esdTree = 0x0;
- AliESD *esd = 0x0;
- AliESDtrack *esdTrack = 0x0;
- for (Int_t imom = 0; imom < fNMom; imom++) {
- //for (Int_t imom = 4; imom < 5; imom++) {
- ref.Reset();
-
- for(Int_t idir = 0; idir<nDirs; idir++){
- // open run loader and load gAlice, kinematics and header
- fRunLoader = AliRunLoader::Open(Form("%s/p%3.1f/galice.root", dir, fTrackMomentum[imom]));
- if (!fRunLoader) {
- AliError(Form("Getting run loader for momentum %3.1f failed.", fTrackMomentum[imom]));
- return kFALSE;
- }
- TString s; s.Form("%s/p%3.1f/", dir, fTrackMomentum[imom]);
- fRunLoader->SetDirName(s);
- fRunLoader->LoadgAlice();
- gAlice = fRunLoader->GetAliRun();
- if (!gAlice) {
- AliError(Form("galice object not found for momentum %3.1f.", fTrackMomentum[imom]));
- return kFALSE;
- }
- fRunLoader->LoadKinematics();
- fRunLoader->LoadHeader();
-
- // open the ESD file
- esdFile = TFile::Open(Form("%s/p%3.1f/AliESDs.root", dir, fTrackMomentum[imom]));
- if (!esdFile || esdFile->IsZombie()) {
- AliError(Form("Opening ESD file failed for momentum", fTrackMomentum[imom]));
- return kFALSE;
- }
- esd = new AliESD;
- esdTree = (TTree*)esdFile->Get("esdTree");
- if (!esdTree) {
- AliError(Form("ESD tree not found for momentum %3.1f.", fTrackMomentum[imom]));
- return kFALSE;
- }
- esdTree->SetBranchAddress("ESD", &esd);
- nTotPart = 0;
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
-
-
- // Event loop
- for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
- fRunLoader->GetEvent(iEvent);
-
- // read stack info
- AliStack* stack = gAlice->Stack();
- TArrayF vertex(3);
- fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
-
- // Load event summary data
- esdTree->GetEvent(iEvent);
- if (!esd) {
- AliWarning(Form("ESD object not found for event %d. [@ momentum %3.1f]", iEvent, imom));
- continue;
- }
-
- for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
- esdTrack = esd->GetTrack(iTrack);
-
- if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
- //if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
- //if(esdTrack->GetConstrainedChi2() > 1E9) continue;
- //if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
- if (esdTrack->GetTRDsignal() == 0.) continue;
-
- // read MC info
- Int_t label = esdTrack->GetLabel();
- if(label<0) continue;
- if (label > stack->GetNtrack()) continue; // background
- TParticle* particle = stack->Particle(label);
- if(!particle){
- AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ event %d track %d]", label, iEvent, iTrack));
- continue;
- }
- if(particle->Pt() < 1.E-3) continue;
- // if (TMath::Abs(particle->Eta()) > 0.3) continue;
- TVector3 dVertex(particle->Vx() - vertex[0],
- particle->Vy() - vertex[1],
- particle->Vz() - vertex[2]);
- if (dVertex.Mag() > 1.E-4){
- //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
- //particle->Print();
- continue;
- }
- Int_t iGen = -1;
- for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
- if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
- iGen = ispec;
- break;
- }
- if(iGen<0) continue;
-
- nPart[iGen]++; nTotPart++;
-
- Float_t mom, length;
- Double_t dedx[AliTRDtrack::kNslice], dEdx;
- Int_t timebin;
- for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){
- // read data for track segment
- for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
- dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice);
- dEdx = esdTrack->GetTRDsignals(iPlane, -1);
- timebin = esdTrack->GetTRDTimBin(iPlane);
-
- // check data
- if ((dEdx <= 0.) || (timebin <= -1.)) continue;
-
- // retrive kinematic info for this track segment
- if(!AliTRDpidESD::GetTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
-
- // find segment length and momentum bin
- Int_t jmom = 1, refMom = -1;
- while(jmom<fNMom-1 && mom>fTrackMomentum[jmom]) jmom++;
- if(TMath::Abs(fTrackMomentum[jmom-1] - mom) < fTrackMomentum[jmom-1] * .2) refMom = jmom-1;
- else if(TMath::Abs(fTrackMomentum[jmom] - mom) < fTrackMomentum[jmom] * .2) refMom = jmom;
- if(refMom<0){
- AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ event %d track %d]", iPlane, iEvent, iTrack));
- continue;
- }
- /*while(jleng<fNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
-
- // this track segment has fulfilled all requierments
- //nPlanePID++;
-
- if(dedx[0] > 0. && dedx[1] > 0.){
- dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
- ref.GetPrincipal(iGen)->AddRow(dedx);
- }
- h1MaxTB[(iGen>0)?1:0]->Fill(timebin);
- } // end plane loop
- } // end track loop
- } // end events loop
-
- delete esd; esd = 0x0;
- esdFile->Close();
- delete esdFile; esdFile = 0x0;
-
- fRunLoader->UnloadHeader();
- fRunLoader->UnloadKinematics();
- delete fRunLoader; fRunLoader = 0x0;
- } // end directory loop
-
- // use data to prepare references
- ref.Prepare2DReferences();
- // save this dEdx references
- ref.SaveReferences(imom, File);
- // save MaxTB references
- SaveMaxTimeBin(imom, File);
-
-
- // print momentum statistics
- printf(" %3.1f ", fTrackMomentum[imom]);
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
- printf("\n");
- } // end momentum loop
-
- TFile *fSave = 0x0;
- TListIter it((TList*)gROOT->GetListOfFiles());
- while((fSave=(TFile*)it.Next()))
- if(strcmp(File, fSave->GetName())==0) break;
-
- fSave->cd();
- fSave->Close();
- delete fSave;
-
- return kTRUE;
}
-//__________________________________________________________________
-void AliTRDCalPIDLQ::SaveMaxTimeBin(const Int_t mom, const char *fn)
-{
- //
- // Save the histograms
- //
+//_________________________________________________________________________
+Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t) const
+{
+ // returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)
- TFile *fSave = 0x0;
- TListIter it((TList*)gROOT->GetListOfFiles());
- TDirectory *pwd = gDirectory;
- Bool_t kFOUND = kFALSE;
- while((fSave=(TFile*)it.Next()))
- if(strcmp(fn, fSave->GetName())==0){
- kFOUND = kTRUE;
- break;
- }
- if(!kFOUND) fSave = new TFile(fn, "RECREATE");
- fSave->cd();
-
- TH1 *h;
- h = (TH1F*)h1MaxTB[0]->Clone(Form("h1MaxTBEL%02d", mom));
- h->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fTrackMomentum[mom]));
- h->GetXaxis()->SetTitle("time [100 ns]");
- h->GetYaxis()->SetTitle("Probability");
- h->Write();
-
- h = (TH1F*)h1MaxTB[1]->Clone(Form("h1MaxTBPI%02d", mom));
- h->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fTrackMomentum[mom]));
- h->GetXaxis()->SetTitle("time [100 ns]");
- h->GetYaxis()->SetTitle("Probability");
- h->Write();
-
- pwd->cd();
+ return spec * AliTRDCalPID::kNMom + mom;
}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff