/**************************************************************************
- * 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$ */
#include "AliLog.h"
#include "AliPID.h"
+#include "../../STAT/TKDPDF.h"
#include "AliTRDCalPIDLQ.h"
#include "AliTRDcalibDB.h"
}
+//_________________________________________________________________________
+Bool_t AliTRDCalPIDLQ::LoadPDF(TDirectoryFile *d)
+{
+ // Read histograms
+ TKDPDF *pdf = 0x0;
+ for (Int_t is=0; is < AliPID::kSPECIES; is++){
+ for (Int_t ip = 0; ip < kNMom; ip++){
+ if(!(pdf = (TKDPDF*)d->Get(Form("%s[%d]", AliPID::ParticleShortName(is), ip)))){
+ AliWarning(Form("Reference for %s[%d] missing.", AliPID::ParticleShortName(is), ip));
+ continue;
+ }
+ fModel->AddAt(pdf->Clone(), GetModelID(ip, is, 0));
+ }
+ }
+}
+ //
//_________________________________________________________________________
Bool_t AliTRDCalPIDLQ::LoadReferences(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));
- }
-
-
+ 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()) {
// 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%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
- hist->Scale(1./hist->Integral());
- fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
+ TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
+ hist->Scale(1./hist->Integral());
+ fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
// if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
//
// 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();
//
if (iType < 0 || iType >= AliPID::kSPECIES) return 0x0;
- if(ip<0 || ip>= kNMom ) return 0x0;
+ if(ip<0 || ip>= kNMom ) return 0x0;
AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fgTrackMomentum[ip]));
, Float_t length, Int_t iplane) 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
+ // 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;
+ 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>fgTrackSegLength[ilength]){
- ilength++;
- }
- Int_t imom = 1;
+ ilength++;
+ }
+ Int_t imom = 1;
while(imom<kNMom-1 && mom>fgTrackMomentum[imom]) imom++;
- Int_t nbinsx, nbinsy;
- TAxis *ax = 0x0, *ay = 0x0;
- Double_t lq1, lq2;
- Double_t mom1 = fgTrackMomentum[imom-1], mom2 = fgTrackMomentum[imom];
- TH2 *hist = 0x0;
- 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.", GetModelID(imom-1, spec, iplane)));
- return 0.;
- }
- ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
- ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
- Float_t x = dedx[0]+dedx[1], y = dedx[2];
+ Int_t nbinsx, nbinsy;
+ TAxis *ax = 0x0, *ay = 0x0;
+ Double_t lq1, lq2;
+ Double_t mom1 = fgTrackMomentum[imom-1], mom2 = fgTrackMomentum[imom];
+ TH2 *hist = 0x0;
+ 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.", GetModelID(imom-1, spec, iplane)));
+ 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));
+ 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*)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.", GetModelID(imom, spec, iplane)));
- return lq1;
- }
- if(kX)
- if(kY) lq2 = hist->GetBinContent( hist->FindBin(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*)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.", GetModelID(imom, spec, iplane)));
+ 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);
-
+ 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
if(mom < fgTrackMomentum[0]) return lq1;
else if(mom > fgTrackMomentum[kNMom-1]) return lq2;
{
// returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)
- return spec * AliTRDCalPID::kNMom + mom;
+ return spec * AliTRDCalPID::kNMom + mom;
}
#include "AliTRDCalPID.h"
+class TDirectoryFile;
class AliTRDCalPIDLQ : public AliTRDCalPID
{
AliTRDCalPIDLQ(const Text_t *name, const Text_t *title);
virtual ~AliTRDCalPIDLQ();
+ Bool_t LoadPDF(TDirectoryFile *);
Bool_t LoadReferences(Char_t* refFile);
TObject* GetModel(Int_t ip, Int_t iType, Int_t iPlane) const;
static Double_t GetLength(Int_t il) { return (il<0 || il>=kNLength) ? -1. : fgTrackSegLength[il]; }
// 2. "file" - output file containing reference data saved in directory
// "dir"
- AliCDBManager *man = AliCDBManager::Instance();
- man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
- man->SetRun(0);
+ AliCDBManager *man = AliCDBManager::Instance();
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetRun(0);
- AliTRDCalPIDRefMaker maker;
- maker.BuildLQReferences(file, dir);
+ AliTRDCalPIDRefMaker maker;
+ maker.BuildLQReferences(file, dir);
}
//___________________________________________________________________
-void generatePIDDB(const char *fileNN = "NN.root", const char *fileLQ = "LQ.root")
+void generatePIDDB()
{
// Write TRD PID DB using the reference data from file "file"
- AliCDBManager *man = AliCDBManager::Instance();
- man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
- man->SetRun(0);
+ AliCDBManager *man = AliCDBManager::Instance();
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetSpecificStorage("TRD/Calib/PIDLQ", "local://.");
+ man->SetRun(0);
- AliCDBStorage *gStorLoc = man->GetStorage("local://$ALICE_ROOT/OCDB");
- if (!gStorLoc) return;
+ AliCDBStorage *gStorLoc = man->GetStorage("local://$ALICE_ROOT/OCDB");
+ if (!gStorLoc) return;
-
- AliTRDCalPID *pidLQ = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
- pidLQ->LoadReferences(fileLQ);
- AliCDBMetaData *md= new AliCDBMetaData();
- md->SetObjectClassName("AliTRDCalPIDLQ");
- md->SetResponsible("Alex Bercuci");
- md->SetBeamPeriod(1);
- md->SetAliRootVersion("v4-06-HEAD"); //root version
- md->SetComment("2D PID for TRD");
- gStorLoc->Put(pidLQ, AliCDBId("TRD/Calib/PIDLQ", 0, 999999999, 0, 1), md);
-
- AliTRDCalPID *pidNN = new AliTRDCalPIDNN("pidNN", "NN TRD PID object");
- pidNN->LoadReferences(fileNN);
- md->SetObjectClassName("AliTRDCalPIDNN");
- md->SetResponsible("Alex Wilk");
- md->SetBeamPeriod(1);
- md->SetAliRootVersion("v4-06-HEAD"); //root version
- md->SetComment("NN PID for TRD");
-
- gStorLoc->Put(pidNN, AliCDBId("TRD/Calib/PIDNN", 0, 999999999, 0, 1), md);
+ AliTRDpidRefMakerLQ pidMaker;
+ TObject *o = pidMaker.GetOCDBEntry("20091101");
+ AliCDBMetaData *md= new AliCDBMetaData();
+ md->SetObjectClassName("AliTRDCalPIDLQ");
+ md->SetResponsible("Alexandru Bercuci");
+ md->SetBeamPeriod(1);
+ md->SetAliRootVersion("v4-16-Release"); //root version
+ md->SetComment("2D PID for TRD");
+ gStorLoc->Put(o, AliCDBId("TRD/Calib/PIDLQ", 0, 999999999, 0, 1), md, AliCDBManager::kReference);
}
//___________________________________________________________________
// > AliTRDCalPID *pid = getPIDObject();
// > pid->GetHistogram(0, 3);
- gStyle->SetOptStat(0);
-
+ gStyle->SetOptStat(0);
+
AliCDBManager *CDBManager = AliCDBManager::Instance();
- CDBManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
- CDBManager->SetRun(0);
+ CDBManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ CDBManager->SetRun(0);
- AliCDBEntry *wrap = CDBManager->Get(Form("TRD/Calib/PID%s", method), 0);
- AliTRDCalPID *pid = dynamic_cast<const AliTRDCalPID *>wrap->GetObject();
- AliCDBMetaData *meta = wrap->GetMetaData();
- if(!pid){
- printf("Error while retriving pid object from DB.\n");
- return 0x0;
- }
- meta->PrintMetaData();
- return pid;
+ AliCDBEntry *wrap = CDBManager->Get(Form("TRD/Calib/PID%s", method), 0);
+ AliTRDCalPID *pid = dynamic_cast<const AliTRDCalPID *>wrap->GetObject();
+ AliCDBMetaData *meta = wrap->GetMetaData();
+ if(!pid){
+ printf("Error while retriving pid object from DB.\n");
+ return 0x0;
+ }
+ meta->PrintMetaData();
+ return pid;
}
memset(fdEdx, 0, 10*sizeof(Float_t));
memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
- memset(fTrain, 0, AliTRDCalPID::kNMom*AliTRDgeometry::kNlayer*sizeof(TEventList*));
- memset(fTest, 0, AliTRDCalPID::kNMom*AliTRDgeometry::kNlayer*sizeof(TEventList*));
-
DefineInput(1, TObjArray::Class());
DefineOutput(1, TTree::Class());
}
AliTRDpidRefMaker::~AliTRDpidRefMaker()
{
if(fReconstructor) delete fReconstructor;
- for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
- for(Int_t il=AliTRDgeometry::kNlayer; il--;){
- if(fTrain[ip][il]) delete fTrain[ip][il];
- if(fTest[ip][il]) delete fTest[ip][il];
- }
- }
}
//AliExternalTrackParam *esd = 0x0;
AliTRDseedV1 *TRDtracklet = 0x0;
for(Int_t itrk=0, nTRD=0; itrk<fTracks->GetEntriesFast(); itrk++){
- // reset the pid information
- memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
-
track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
if(!track->HasESDtrack()) continue;
ULong_t status = track->GetStatus();
if(HasMCdata()) labelsacc[nTRD++] = track->GetLabel();
-
-
- switch(fRefPID){
- case kV0:
- SetRefPID(TRDtrack,fPID);
- break;
- case kMC:
- if(!HasMCdata()){
- AliError("Could not retrive reference PID from MC");
- return;
- }
- switch(track -> GetPDG()){
- case kElectron:
- case kPositron:
- fPID[AliPID::kElectron] = 1.;
- break;
- case kMuonPlus:
- case kMuonMinus:
- fPID[AliPID::kMuon] = 1.;
- break;
- case kPiPlus:
- case kPiMinus:
- fPID[AliPID::kPion] = 1.;
- break;
- case kKPlus:
- case kKMinus:
- fPID[AliPID::kKaon] = 1.;
- break;
- case kProton:
- case kProtonBar:
- fPID[AliPID::kProton] = 1.;
- break;
- }
- break;
- default:
- AliWarning("PID reference source not implemented");
- return;
- }
- // set reconstructor
- //Float_t *dedx;
- TRDtrack -> SetReconstructor(fReconstructor);
+ // fill the pid information
+ memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
+ switch(fRefPID){
+ case kV0: SetRefPID(kV0, 0x0/*v0*/, fPID); break;
+ case kMC: SetRefPID(kMC, track, fPID); break;
+ case kRec: SetRefPID(kRec, TRDtrack, fPID); break;
+ }
// fill the momentum and dE/dx information
+ TRDtrack -> SetReconstructor(fReconstructor);
for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
- if(!GetdEdx(TRDtracklet)) continue;
+ if(!CookdEdx(TRDtracklet)) continue;
switch(fRefP){
case kMC:
if(!HasMCdata()){
//________________________________________________________________________
-void AliTRDpidRefMaker::LoadFile(const Char_t *InFile)
+void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid)
{
- //
- // Loads the files and sets the event list
- // for neural network training and
- // building of the 2-dim reference histograms.
- // Usable for training outside of the makeResults.C macro
- //
-
- TFile::Open(InFile, "READ");
- fData = (TTree*)gFile->Get(GetName());
-
- for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- for(Int_t iLy = 0; iLy < AliTRDgeometry::kNlayer; iLy++){
- fTrain[iMom][iLy] = new TEventList(Form("fTrain_P%d_L%d", iMom, iLy), Form("Training list for momentum bin %d layer %d", iMom, iLy));
- fTest[iMom][iLy] = new TEventList(Form("fTest_P%d_L%d", iMom, iLy), Form("Test list for momentum bin %d layer %d", iMom, iLy));
+// !!!! PREMILMINARY FUNCTION !!!!
+//
+// this is the place for the V0 procedure
+// as long as there is no one implemented,
+// just the probabilities
+// of the TRDtrack are used!
+
+ switch(select){
+ case kV0:
+ {
+ AliTRDv0Info *v0 = static_cast<AliTRDv0Info*>(source);
+ v0->Print(); // do something
}
- }
-}
-
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::LoadContainer(const Char_t *InFileCont)
-{
-
- //
- // Loads the container if no container is there.
- // Useable for training outside of the makeResults.C macro
- //
-
- TFile::Open(InFileCont, "READ");
- fContainer = (TObjArray*)gFile->Get(GetName());
-}
-
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::SetRefPID(void *source, Float_t *pid)
-{
- // !!!! PREMILMINARY FUNCTION !!!!
- //
- // this is the place for the V0 procedure
- // as long as there is no one implemented,
- // just the probabilities
- // of the TRDtrack are used!
-
- AliTRDtrackV1 *TRDtrack = static_cast<AliTRDtrackV1*>(source);
- TRDtrack -> SetReconstructor(fReconstructor);
- //fReconstructor -> SetOption("nn");
- TRDtrack -> CookPID();
- for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
- pid[iPart] = TRDtrack -> GetPID(iPart);
- AliDebug(4, Form("PDG is (in V0info) %d %f", iPart, pid[iPart]));
+ break;
+ case kMC:
+ if(!HasMCdata()){
+ AliError("Could not retrive reference PID from MC");
+ return;
+ }
+ {
+ AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
+ switch(track->GetPDG()){
+ case kElectron:
+ case kPositron:
+ fPID[AliPID::kElectron] = 1.;
+ break;
+ case kMuonPlus:
+ case kMuonMinus:
+ fPID[AliPID::kMuon] = 1.;
+ break;
+ case kPiPlus:
+ case kPiMinus:
+ fPID[AliPID::kPion] = 1.;
+ break;
+ case kKPlus:
+ case kKMinus:
+ fPID[AliPID::kKaon] = 1.;
+ break;
+ case kProton:
+ case kProtonBar:
+ fPID[AliPID::kProton] = 1.;
+ break;
+ }
+ }
+ break;
+ case kRec:
+ {
+ AliTRDtrackV1 *TRDtrack = static_cast<AliTRDtrackV1*>(source);
+ TRDtrack -> SetReconstructor(fReconstructor);
+ //fReconstructor -> SetOption("nn");
+ TRDtrack -> CookPID();
+ for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
+ pid[iPart] = TRDtrack -> GetPID(iPart);
+ AliDebug(4, Form("PDG is (in V0info) %d %f", iPart, pid[iPart]));
+ }
+ }
+ break;
+ default:
+ AliWarning("PID reference source not implemented");
+ return;
}
}
class TTree;
class TObjArray;
-class TEventList;
class AliTRDReconstructor;
class AliTRDseedV1;
class AliTRDpidRefMaker : public AliTRDrecoTask
void CreateOutputObjects();
void Exec(Option_t *option);
- void LoadContainer(const Char_t *InFileCont);
- void LoadFile(const Char_t *InFile);
-
void SetAbundance(Float_t train, Float_t test);
- void SetRefPID(void *source, Float_t *pid);
+ void SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid);
void SetSource(ETRDpidRefMakerSource pid, ETRDpidRefMakerSource momentum) {fRefPID = pid; fRefP = momentum;}
void Terminate(Option_t *);
protected:
- virtual Float_t* GetdEdx(AliTRDseedV1*) = 0;
+ virtual Float_t* CookdEdx(AliTRDseedV1*) = 0;
virtual Int_t GetNslices() = 0;
virtual void Fill();
AliTRDReconstructor *fReconstructor; //! reconstructor needed for recalculation the PID
TObjArray *fV0s; //! v0 array
TTree *fData; //! dEdx-P data
- TEventList *fTrain[AliTRDCalPID::kNMom][AliTRDgeometry::kNlayer]; // Training list for each momentum
- TEventList *fTest[AliTRDCalPID::kNMom][AliTRDgeometry::kNlayer]; // Test list for each momentum
ETRDpidRefMakerSource fRefPID; // reference PID source
ETRDpidRefMakerSource fRefP; // reference momentum source
Float_t fTrainFreq; //! training sample relative abundance
//________________________________________________________________________
-Float_t* AliTRDpidRefMakerLQ::GetdEdx(AliTRDseedV1 *trklt)
+Float_t* AliTRDpidRefMakerLQ::CookdEdx(AliTRDseedV1 *trklt)
{
trklt->CookdEdx(AliTRDpidUtil::kLQslices);
const Float_t *dedx = trklt->GetdEdx();
return fdEdx;
}
+#include "../Cal/AliTRDCalPIDLQ.h"
//__________________________________________________________________
-Bool_t AliTRDpidRefMakerLQ::GenerateOCDBEntry(Option_t *)
+TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
{
- return kTRUE;
+ TDirectoryFile *d = 0x0;
+ if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
+ if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
+ AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
+ cal->LoadPDF(d);
+ return cal;
}
//__________________________________________________________________
TList *l=gPad->GetListOfPrimitives();
for(Int_t is=0; is<AliPID::kSPECIES; is++){
((TVirtualPad*)l->At(is))->cd();
- ((TH2*)arr->At(is))->Draw("colz");
+ ((TH2*)arr->At(is))->Draw("cont4z");
}
return kTRUE;
((TH2*)fContainer->At(0))->Fill(fSbin, fPbin);
TH2* h2 = (TH2*)((TObjArray*)fContainer->At(1+fPbin))->At(fSbin);
h2->Fill(fdEdx[0], fdEdx[1]);
- //printf("h[%s] : [%f] [%f]\n", h2->GetName(), fdEdx[0], fdEdx[1]);
}
//________________________________________________________________________
Bool_t AliTRDpidRefMakerLQ::PostProcess()
{
- TFile::Open(Form("TRD.Calib%s.root", GetName()));
+// Analyse merged dedx = f(p) distributions.
+// - select momentum - species bins
+// - rotate to principal components
+// - locally interpolate with TKDPDF
+// - save interpolation to monitoring histograms
+// - write pdf to file for loading to OCDB
+//
+
+
+ TFile *fCalib = TFile::Open(Form("TRD.Calib%s.root", GetName()), "update");
fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
if (!fData) {
AliError("Tree not available");
return kFALSE;
}
+ TDatime d;
+ TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
+ pdfs->Write();
AliDebug(2, Form("Data[%d]", fData->GetEntries()));
- // save PDF representation
+ // save a volatile PDF representation
+ gROOT->cd();
TH2 *h2 = 0x0;
fResults = new TObjArray(AliTRDCalPID::kNMom);
for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
TObjArray *arr = new TObjArray(AliPID::kSPECIES);
arr->SetName(Form("Pbin%02d", ip));
for(Int_t is=AliPID::kSPECIES; is--;) {
- h2 = new TH2I(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, -6., 6., 50, -6., 6.);
+ h2 = new TH2D(Form("i%s%d", AliPID::ParticleShortName(is), ip), Form("Interpolated %s dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, -3., 3., 50, -3., 3.);
h2->GetXaxis()->SetTitle("log(dE/dx^{*}_{am}) [au]");
h2->GetYaxis()->SetTitle("log(dE/dx^{*}_{dr}) [au]");
h2->GetZaxis()->SetTitle("#");
}
fResults->AddAt(arr, ip);
}
-
+ pdfs->cd();
TCanvas *cc = new TCanvas("cc", "", 500, 500);
Float_t *data[] = {0x0, 0x0};
TPrincipal principal(2, "ND");
- for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
+ for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){
for(Int_t is=AliPID::kSPECIES; is--;) {
principal.Clear();
Int_t n = fData->Draw("dEdx[0]:dEdx[1]", Form("p==%d&&s==%d", ip, is), "goff");
- AliDebug(2, Form("pBin[%d] sBin[%d] n[%d]", ip, is, n));
- if(n<1000/*Int_t(kMinStat)*Int_t(kMinBuckets)*/){
- AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
+
+ // estimate bucket statistics
+ Int_t nb(kMinBuckets), // number of buckets
+ ns(Int_t(Float_t(n)/nb)); //statistics/bucket
+
+// if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
+// else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
+
+ AliDebug(2, Form("pBin[%d] sBin[%d] N[%d] n[%d] nb[%d]", ip, is, n, ns, nb));
+ if(ns<Int_t(kMinStat)){
+ AliWarning(Form("Not enough entries [%d] for %s[%d].", ns, AliPID::ParticleShortName(is), ip));
continue;
}
// allocate storage
data[0][n]=rxy[0]; data[1][n]=rxy[1];
}
- // estimate bucket statistics
- Int_t ns(kMinStat), //statistics/bucket
- nb(kMinBuckets); // number of buckets
- if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
- else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
-
// build PDF
TKDPDF pdf(n, 2, ns, data);
pdf.SetStore();
rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
pdf.Eval(rxy, r, e, kTRUE);
}
- pdf.DrawBins(0,1,-6,6,-6,6);cc->Modified(); cc->Update();
+ // visual on-line monitoring
+ pdf.DrawBins(0,1,-5.,5.,-8.,8.);cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
-
- // save a discretization of the PDF for monitoring
+ // save a discretization of the PDF for result monitoring
TH2 *h2s = (TH2D*)((TObjArray*)fResults->At(ip))->At(is);
TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
for(int ix=1; ix<=ax->GetNbins(); ix++){
}
}
-
- //pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
-
+ // write results to output array
+ pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
delete [] data[0]; delete [] data[1];
}
}
+ pdfs->Write();
+ fCalib->Close(); delete fCalib;
return kTRUE; // testing protection
}
class AliTRDpidRefMakerLQ : public AliTRDpidRefMaker {
public:
enum ETRDpidRefMakerLQsteer{
- kMinStat = 100 // minimum statistics/bucket 10%
+ kMinStat = 50 // minimum statistics/bucket 14%
,kMinBuckets = 450 // minimum number of buckets [lambda(6)*alpha(1.5)*regions(50)]
};
AliTRDpidRefMakerLQ();
~AliTRDpidRefMakerLQ();
void CreateOutputObjects();
- Bool_t GenerateOCDBEntry(Option_t *);
+ TObject* GetOCDBEntry(Option_t *);
Bool_t GetRefFigure(Int_t ifig);
Bool_t PostProcess();
protected:
- Float_t* GetdEdx(AliTRDseedV1*);
+ Float_t* CookdEdx(AliTRDseedV1*);
Int_t GetNslices() { return 2;}
void Fill();
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff