}
-//_____________________________________________________________________________
-void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader
- , AliRawReader *rawReader) const
-{
- //
- // Reconstruct clusters
- //
-
- AliInfo("Reconstruct TRD clusters from RAW data [RunLoader, RawReader]");
-
- AliLoader *loader = runLoader->GetLoader("TRDLoader");
- loader->LoadRecPoints("recreate");
-
- runLoader->CdGAFile();
- Int_t nEvents = runLoader->GetNumberOfEvents();
-
- rawReader->Reset();
- rawReader->Select("TRD");
-
- for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
-
- if (!rawReader->NextEvent()) break;
-
- // New (fast) cluster finder
- AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
- clusterer.Open(runLoader->GetFileName(),iEvent);
- clusterer.Raw2ClustersChamber(rawReader);
-
- clusterer.WriteClusters(-1);
-
- }
-
- loader->UnloadRecPoints();
-
-}
-
//_____________________________________________________________________________
void AliTRDReconstructor::Reconstruct(AliRawReader *rawReader
, TTree *clusterTree) const
}
//_____________________________________________________________________________
-void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader) const
-{
- //
- // Reconstruct clusters
- //
-
- AliInfo("Reconstruct TRD clusters [AliRunLoader]");
- AliLoader *loader = runLoader->GetLoader("TRDLoader");
- loader->LoadRecPoints("recreate");
-
- runLoader->CdGAFile();
- Int_t nEvents = runLoader->GetNumberOfEvents();
-
- for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
- AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
- clusterer.Open(runLoader->GetFileName(),iEvent);
- clusterer.ReadDigits();
- clusterer.MakeClusters();
- clusterer.WriteClusters(-1);
- }
-
- loader->UnloadRecPoints();
-
-}
-
-//_____________________________________________________________________________
-AliTracker *AliTRDReconstructor::CreateTracker(AliRunLoader *runLoader) const
+AliTracker *AliTRDReconstructor::CreateTracker() const
{
//
// Create a TRD tracker
//
- runLoader->CdGAFile();
-
- return new AliTRDtracker(gFile);
-
-}
-
-//_____________________________________________________________________________
-void AliTRDReconstructor::FillESD(AliRunLoader* /*runLoader*/
- , AliRawReader* /*rawReader*/
- , AliESDEvent* /*esd*/) const
-{
- //
- // Fill ESD
- //
-
-}
-
-//_____________________________________________________________________________
-void AliTRDReconstructor::FillESD(AliRawReader* /*rawReader*/
- , TTree* /*clusterTree*/
- , AliESDEvent* /*esd*/) const
-{
- //
- // Fill ESD
- //
+ return new AliTRDtracker(NULL);
}
//
}
-
-//_____________________________________________________________________________
-void AliTRDReconstructor::FillESD(AliRunLoader* /*runLoader*/
- , AliESDEvent* /*esd*/) const
-{
- //
- // Fill ESD
- //
-
-}
virtual Bool_t HasDigitConversion() const { return kFALSE; };
virtual void ConvertDigits(AliRawReader *rawReader, TTree *digitsTree) const;
- virtual Bool_t HasLocalReconstruction() const { return kTRUE; };
- virtual void Reconstruct(AliRunLoader *runLoader, AliRawReader *rawReader) const;
virtual void Reconstruct(AliRawReader *rawReader, TTree *clusterTree) const;
virtual void Reconstruct(TTree *digitsTree, TTree *clusterTree) const;
- virtual void Reconstruct(AliRunLoader *runLoader) const;
- virtual AliTracker *CreateTracker(AliRunLoader *runLoader) const;
+ virtual AliTracker *CreateTracker() const;
- virtual void FillESD(AliRunLoader *runLoader, AliRawReader *rawReader, AliESDEvent *esd) const;
- virtual void FillESD(AliRawReader *rawReader, TTree *clusterTree, AliESDEvent *esd) const;
+ virtual void FillESD(AliRawReader */*rawReader*/, TTree *clusterTree, AliESDEvent *esd) const
+ {FillESD((TTree*)NULL,clusterTree,esd);}
virtual void FillESD(TTree *digitsTree, TTree *clusterTree, AliESDEvent *esd) const;
- virtual void FillESD(AliRunLoader *runLoader, AliESDEvent *esd) const;
static void SetSeedingOn(Bool_t seeding) { fgkSeedingOn = seeding; }
static void SetStreamLevel(Int_t level) { fgStreamLevel = level; }
case kIDFEE :
return CacheCDBEntry(kIDFEE ,"TRD/Calib/FEE");
break;
+ case kIDPIDNN :
+ return CacheCDBEntry(kIDPIDNN ,"TRD/Calib/PIDNN");
case kIDPIDLQ :
return CacheCDBEntry(kIDPIDLQ ,"TRD/Calib/PIDLQ");
break;
}
//_____________________________________________________________________________
-const AliTRDCalPID *AliTRDcalibDB::GetPIDLQObject()
+const AliTRDCalPID *AliTRDcalibDB::GetPIDObject(const Int_t method)
{
//
// Returns the object storing the distributions for PID with likelihood
//
-
- return dynamic_cast<const AliTRDCalPID *>
- (GetCachedCDBObject(kIDPIDLQ));
+ switch(method){
+ case 0: return dynamic_cast<const AliTRDCalPID *>
+ (GetCachedCDBObject(kIDPIDNN));
+ case 1: return dynamic_cast<const AliTRDCalPID *>
+ (GetCachedCDBObject(kIDPIDLQ));
+ }
}
Bool_t IsChamberMasked(Int_t det);
const AliTRDCalMonitoring *GetMonitoringObject();
- const AliTRDCalPID *GetPIDLQObject();
+ const AliTRDCalPID *GetPIDObject(const Int_t method);
// Related functions, these depend on calibration data
static Float_t GetOmegaTau(Float_t vdrift, Float_t bz);
protected:
// For caching see also implentation of GetCachedCDBObject in the .cxx file
- enum { kCDBCacheSize = 15 }; // Number of cached objects
+ enum { kCDBCacheSize = 16 }; // Number of cached objects
enum { kIDVdriftPad = 0
, kIDVdriftChamber
, kIDT0Pad
, kIDChamberPos
, kIDStackPos
, kIDSuperModulePos
+ , kIDPIDNN
, kIDPIDLQ
, kIDMonitoringData
, kIDChamberStatus
AliError("No instance of AliTRDcalibDB.");
return;
}
- const AliTRDCalPID *pd = calibration->GetPIDLQObject();
+ const AliTRDCalPID *pd = calibration->GetPIDObject(1);
AliTRDltuTracklet *trk;
Int_t nTracklets = GetNtracklets();
dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
Float_t length = 3.7;
- probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length);
- probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length);
+ probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length, 0);
+ probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length, 0);
}
}
// Retrieve the CDB container class with the probability distributions
- const AliTRDCalPID *pd = calibration->GetPIDLQObject();
+ const AliTRDCalPID *pd = calibration->GetPIDObject(1);
if (!pd) {
AliErrorGeneral("AliTRDpidESD::MakePID()"
,"No access to AliTRDCalPID");
// Get the probabilities for the different particle species
for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
- p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length);
+ p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iPlan);
//p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
}
}
}
//_____________________________________________________________________________
-void AliTRDtrack::CookdEdxTimBin()
+void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
{
//
// Set fdEdxPlane and fTimBinPlane and also get the
}
+//_____________________________________________________________________________
+void AliTRDtrack::CookdEdxNN(Float_t *dedx){
+
+ // This function calcuates the input for the neural networks
+ // which are used for the PID.
+
+ Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();//number of time bins in chamber
+ Int_t plane; // plane of current cluster
+ Int_t tb; // time bin of current cluster
+ Int_t slice; // curent slice
+ AliTRDcluster *cluster = 0x0; // pointer to current cluster
+ const Int_t lMLPscale = 16000; // scaling of the MLP input to be smaller than 1
+
+ // Reset class and local contors/variables
+ for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++){
+ for (Int_t iSlice = 0; iSlice < kNMLPslice; iSlice++) {
+ *(dedx + (kNMLPslice * iPlane) + iSlice) = 0.;
+ }
+ }
+
+ // start looping over clusters attached to this track
+ for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
+ cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
+ if(!cluster) continue;
+
+ // Read info from current cluster
+ plane = AliTRDgeometry::GetPlane(cluster->GetDetector());
+ if (plane < 0 || plane >= AliESDtrack::kNPlane) {
+ AliError(Form("Wrong plane %d", plane));
+ continue;
+ }
+
+ tb = cluster->GetLocalTimeBin();
+ if(tb == 0 || tb >= ntb){
+ AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus));
+ AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
+ continue;
+ }
+
+ slice = tb * kNMLPslice / ntb;
+
+ *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/lMLPscale;
+ } // End of loop over cluster
+
+}
+
+
//_____________________________________________________________________________
void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
{
}
//_____________________________________________________________________________
-Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
+Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
{
//
// This function calculates the PID probabilities based on TRD signals
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
AliError("No access to calibration data");
- return -1;
+ return kFALSE;
}
- // Retrieve the CDB container class with the probability distributions
- const AliTRDCalPID *pd = calibration->GetPIDLQObject();
- if (!pd) {
- AliError("No access to AliTRDCalPID");
- return -1;
- }
+ // Retrieve the CDB container class with the probability distributions
+ const AliTRDCalPID *pd = calibration->GetPIDObject(fPIDmethod == kNN ? 0 : 1);
+ if (!pd) {
+ AliError("No access to AliTRDCalPID");
+ return kFALSE;
+ }
+ // Calculate the input for the NN if fPIDmethod is kNN
+ Float_t ldEdxNN[AliTRDCalPID::kNPlane * kNMLPslice], *dedx = 0x0;
+ if(fPIDmethod == kNN) {
+ CookdEdxNN(&ldEdxNN[0]);
+ }
- Double_t p0 = 1./AliPID::kSPECIES;
- if(AliPID::kSPECIES != 5){
- AliError("Probabilities array defined only for 5 values. Please modify !!");
- return -1;
- }
- Double_t p[] = {p0, p0, p0, p0, p0};
- Float_t length; // track segment length in chamber
- Int_t nPlanePID = 0;
- // Skip tracks which have no TRD signal at all
- if (fdEdx == 0.) return -1;
-
- for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
+ // Sets the a priori probabilities
+ for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
+ fPID[ispec] = 1./AliPID::kSPECIES;
+ }
- length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
- // check data
- if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
- || fTimBinPlane[iPlane] == -1.) continue;
+ if(AliPID::kSPECIES != 5){
+ AliError("Probabilities array defined only for 5 values. Please modify !!");
+ return kFALSE;
+ }
- // this track segment has fulfilled all requierments
- nPlanePID++;
- // Get the probabilities for the different particle species
- for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
- p[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], fdEdxPlane[iPlane], length);
- //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
- }
- }
- if(nPlanePID == 0) return 0;
+ pidQuality = 0;
+ Float_t length; // track segment length in chamber
- // normalize probabilities
- Double_t probTotal = 0.;
- for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
+ // Skip tracks which have no TRD signal at all
+ if (fdEdx == 0.) return kFALSE;
+
+ for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
+
+ length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
+
+ // check data
+ if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
+ || fTimBinPlane[iPlane] == -1.) continue;
+
+ // this track segment has fulfilled all requierments
+ pidQuality++;
+
+ // Get the probabilities for the different particle species
+ for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
+ switch(fPIDmethod){
+ case kLQ:
+ dedx = fdEdxPlane[iPlane];
+ break;
+ case kNN:
+ dedx = &ldEdxNN[iPlane*kNMLPslice];
+ break;
+ }
+ fPID[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane);
+ }
+ }
+ if(pidQuality == 0) return kTRUE;
- if(probTotal <= 0.) {
- AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
- return -1;
- }
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
+ // normalize probabilities
+ Double_t probTotal = 0.;
+ for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += fPID[iSpecies];
+ if(probTotal <= 0.) {
+ AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
+ return kFALSE;
+ }
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
- // book PID to the ESD track
- esd->SetTRDpid(p);
- esd->SetTRDpidQuality(nPlanePID);
-
- return 0;
+ return kTRUE;
}
+
//_____________________________________________________________________________
Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
{
//_____________________________________________________________________________
Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
- , Double_t h01, Int_t /*plane*/)
+ , Double_t h01, Int_t /*plane*/, Int_t tid)
{
//
// Assignes found cluster to the track and updates track information
, kNplane = 6
, kNcham = 5
, kNsect = 18
- , kNslice = 3};
+ , kNslice = 3
+ , kNMLPslice = 8};
+
+ enum AliTRDPIDMethod {
+ kNN = 0
+ , kLQ = 1
+ };
AliTRDtrack();
AliTRDtrack(/*const */AliTRDcluster *c, Int_t index, const Double_t xx[5], const Double_t cc[15], Double_t xr, Double_t alpha);
Int_t GetClusterIndex(Int_t i) const { return fIndex[i]; }
Double_t GetC() const { return AliExternalTrackParam::GetC(GetBz()); }
Double_t GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const;
+
+ inline AliTRDPIDMethod GetPIDMethod() const {return fPIDmethod;}
+ inline void SetPIDMethod(AliTRDPIDMethod method) {fPIDmethod = method;}
+ // end
Float_t GetPIDsignals(Int_t iPlane, Int_t iSlice) const { return fdEdxPlane[iPlane][iSlice]; }
Int_t GetPIDTimBin(Int_t i) const { return fTimBinPlane[i]; }
Double_t GetLikelihoodElectron() const { return fLhElectron; }
// A. Bercuci new functions
void SetTrackSegmentDirMom(const Int_t plane);
void CookdEdx(Double_t low = 0.05, Double_t up = 0.7);
- void CookdEdxTimBin();
- Int_t CookPID(AliESDtrack *p);
+ void CookdEdxTimBin(const Int_t tid);
+ Bool_t CookPID(Int_t &pidQuality);
void SetCluster(AliTRDcluster* cl, Int_t index = -1) {fClusters[(index == -1) ? GetNumberOfClusters()-1 : index] = cl;}
inline AliTRDcluster* GetCluster(Int_t layer){ return (layer >= 0 && layer < GetNumberOfClusters()) ? fClusters[layer] : 0x0;}
inline Float_t GetMomentumPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fMom[plane] : 0.;}
+ inline const Double_t* GetPID() const {return fPID;}
inline Float_t GetSnpPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fSnp[plane] : 0.;}
inline Float_t GetTglPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fTgl[plane] : 0.;}
Float_t GetTrackLengthPlane(Int_t plane) const;
Bool_t Rotate(Double_t angle, Bool_t absolute = kFALSE);
Float_t StatusForTOF();
Bool_t Update(const AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01);
- Int_t UpdateMI(/*const */AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01, Int_t plane);
+ Int_t UpdateMI(/*const */AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01, Int_t plane, Int_t tid=0);
protected:
AliTRDtrack &operator=(const AliTRDtrack &t);
+ void CookdEdxNN(Float_t *dedx);
Double_t GetBz() const;
Bool_t Update(const AliCluster */*c*/, Double_t /*chi2*/, Int_t /*idx*/) { return 0; }
Double_t GetPredictedChi2(const AliCluster* /*c*/) const { return 0.0; }
Float_t fDE; // Integrated delta energy
Float_t fdEdxPlane[kNplane][kNslice]; // dE/dx from all 6 planes in 3 slices each
Int_t fTimBinPlane[kNplane]; // Time bin of Max cluster from all 6 planes
+ Double_t fPID[AliPID::kSPECIES]; // PID probabilities
// A.Bercuci
Float_t fMom[kNplane]; // Track momentum at chamber entrance
AliTRDcluster* fClusters[kMAXCLUSTERSPERTRACK];
Bool_t fClusterOwner; // indicates the track is owner of cluster
// end A.Bercuci
-
+ AliTRDPIDMethod fPIDmethod; // switch between different PID methods
Bool_t fStopped; // Track stop indication
Int_t fIndex[kMAXCLUSTERSPERTRACK]; // Global indexes of clusters
Int_t fIndexBackup[kMAXCLUSTERSPERTRACK]; // Backup indexes of clusters - used in iterations
Float_t fBudget[3]; // Integrated material budget
AliTRDtrack *fBackupTrack; //! Backup track
- ClassDef(AliTRDtrack,8) // TRD reconstructed tracks
+ ClassDef(AliTRDtrack,9) // TRD reconstructed tracks
};
Int_t foundClr = track->GetNumberOfClusters();
if (foundClr >= foundMin) {
track->CookdEdx();
- track->CookdEdxTimBin(); // A.Bercuci 25.07.07
+ track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
CookLabel(track,1 - fgkLabelFraction);
if (track->GetBackupTrack()) {
UseClusters(track->GetBackupTrack());
Float_t foundMin = fgkMinClustersInTrack * timeBins;
Int_t nseed = 0;
Int_t found = 0;
+ Int_t pidQ = 0;
//Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
AliTRDtrack seed2;
//AliTRDtrack *pt = seed2;
//AliTRDtrack &t = *pt;
FollowProlongation(*pt);
- if (pt->GetNumberOfClusters() >= foundMin) {
+ //if (pt->GetNumberOfClusters() >= foundMin) {
//UseClusters(&t);
//CookLabel(pt, 1-fgkLabelFraction);
pt->CookdEdx();
- pt->CookdEdxTimBin();
- pt->CookPID(seed);
+ pt->CookdEdxTimBin(seed->GetID());
+ pt->SetPIDMethod(AliTRDtrack::kLQ); //switch between TRD PID methods
+ pt->CookPID(pidQ);
//pt->Calibrate(); // slot for calibration
- }
+ //}
found++;
Double_t xTPC = 250.0;
delete seed2;
if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
pt2->CookdEdx( );
- pt2->CookdEdxTimBin();
+ pt2->CookdEdxTimBin(seed->GetID());
seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
fHRefit->Fill(6);
// The size of output array has is 2*n
//
- if (n <= 0)
+ if (n <= 0) {
return 0;
+ }
Int_t *sindexS = new Int_t[n]; // Temporary array for sorting
Int_t *sindexF = new Int_t[2*n];
}
else {
track->CookdEdx();
- track->CookdEdxTimBin();
+ track->CookdEdxTimBin(-1);
CookLabel(track,0.9);
}
//////////////////////////////////////////////////////////////////////////
// //
-// 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 <shukla@pi0.physi.uni-heidelberg.de> //
ClassImp(AliTRDCalPID)
-Char_t* AliTRDCalPID::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
+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};
-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};
-
-
//_________________________________________________________________________
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);
-
}
//_________________________________________________________________________
// Destructor
//
- CleanUp();
-
-}
-
-//_________________________________________________________________________
-void AliTRDCalPID::CleanUp()
-{
- //
- // Delets all newly created objects
- //
-
- if (fHistdEdx) {
- delete fHistdEdx;
- fHistdEdx = 0x0;
- }
-
- if (fHistTimeBin) {
- delete fHistTimeBin;
- fHistTimeBin = 0x0;
- }
-}
-
-//_____________________________________________________________________________
-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)
-{
- //
- // 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 = TFile::Open(refFile, "READ");
- if (!histFile || !histFile->IsOpen()) {
- AliError(Form("Opening TRD histgram file %s failed", refFile));
- return kFALSE;
+ if (fModel) {
+ delete fModel;
}
- gROOT->cd();
-
- // 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(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
- 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 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(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);
}
-
///////////////////////////////////////////////////////////////////////////////
// //
-// 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> //
+// //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Alex Wilk <wilka@uni-muenster.de> //
+// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> //
// //
///////////////////////////////////////////////////////////////////////////////
#include <TNamed.h>
#endif
-class TH1;
-class TObjArray;
-class AliTRDCalPID : public TNamed {
-public:
- enum {
- kNMom = 11,
- kNLength = 4
- };
+class AliTRDCalPID : public TNamed
+{
+
+ public:
+
+ enum {
+ kNMom = 11,
+ kNPlane = 6
+ };
AliTRDCalPID();
AliTRDCalPID(const Text_t *name, const Text_t *title);
- AliTRDCalPID(const AliTRDCalPID& pd);
- virtual ~AliTRDCalPID();
- AliTRDCalPID& operator=(const AliTRDCalPID &c);
- virtual void Copy(TObject &c) const;
-
- Bool_t LoadLQReferences(Char_t* refFile);
- Bool_t LoadNNReferences(Char_t* /*refFile*/) {return kTRUE;}
- //void SetMeanChargeRatio(Double_t ratio) { fMeanChargeRatio = ratio; }
-
- //Double_t GetMeanChargeRatio() const { return fMeanChargeRatio; }
- static Double_t GetMomentum(Int_t ip) { return (ip<0 || ip>=kNMom) ? -1. : fTrackMomentum[ip]; }
- static Double_t GetLength(Int_t il) { return (il<0 || il>=kNLength) ? -1. : fTrackSegLength[il]; }
- //Double_t GetMean(Int_t iType, Int_t ip) const;
- //Double_t GetNormalization(Int_t iType, Int_t ip) const;
-
- TH1* GetHistogram(Int_t iType, Int_t ip/*, Int_t il*/) const;
- TH1* GetHistogramT(Int_t iType, Int_t ip) const;
- Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const;
- Double_t GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const;
+ virtual ~AliTRDCalPID();
+
+ virtual Bool_t LoadReferences(Char_t *refFile) = 0;
+ static Double_t GetMomentum(Int_t ip)
+ { return (ip<0 || ip>=kNMom)
+ ? -1.0
+ : fTrackMomentum[ip];
+ }
+ virtual TObject *GetModel(Int_t ip, Int_t iType, Int_t iPlane) const = 0;
+ virtual Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx
+ , Float_t length, Int_t plane) const = 0;
+ static Char_t *GetPartName(Int_t i) { return fPartName[i]; }
+ static Char_t *GetPartSymb(Int_t i) { return fPartSymb[i]; }
+
+ void SetPartName(Int_t i, Char_t *name) { fPartName[i] = name; }
+ void SetPartSymb(Int_t i, Char_t *symb) { fPartSymb[i] = symb; }
protected:
- void Init();
- inline Int_t GetHistID(Int_t part, Int_t mom/*, Int_t length=0*/) const { return part*kNMom + mom; }
- void CleanUp();
-
- public:
- static Char_t *fpartName[5]; //! Names of particle species
- static Char_t *fpartSymb[5]; //! Symbols of particle species
-
+ virtual void Init() = 0;
+ virtual Int_t GetModelID(Int_t mom, Int_t spec, Int_t plane) const = 0;
+
+ private:
+
+ AliTRDCalPID(const AliTRDCalPID& pd);
+ AliTRDCalPID &operator=(const AliTRDCalPID &c);
+
protected:
- static Float_t fTrackMomentum[kNMom]; // Track momenta for which response functions are available
- static Float_t fTrackSegLength[kNLength]; // Track segment lengths for which response functions are available
- Double_t fMeanChargeRatio; // Ratio of mean charge from real Det. to prob. dist.
- TObjArray *fHistdEdx; // Prob. of dEdx for 5 particles and for several momenta
- TObjArray *fHistTimeBin; // Prob. of max time bin for 5 particles and for several momenta
-
- ClassDef(AliTRDCalPID, 1) // The TRD PID response container class
+ static Char_t *fPartName[5]; //! Names of particle species
+ static Char_t *fPartSymb[5]; //! Symbols of particle species
-};
+ static Float_t fTrackMomentum[kNMom]; // Track momenta for which response functions are available
+ TObjArray *fModel; // Model for probability estimate
-#endif
+ ClassDef(AliTRDCalPID, 3) // Base class for TRD PID methods
+};
+#endif
--- /dev/null
+/**************************************************************************
+ * 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 //
+// //
+// Authors: //
+// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> //
+// Alex Bercuci <a.bercuci@gsi.de> //
+// //
+//////////////////////////////////////////////////////////////////////////
+
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TFile.h>
+#include <TROOT.h>
+
+#include "AliLog.h"
+#include "AliPID.h"
+
+#include "AliTRDCalPIDLQ.h"
+#include "AliTRDcalibDB.h"
+
+ClassImp(AliTRDCalPIDLQ)
+
+Float_t AliTRDCalPIDLQ::fTrackSegLength[kNLength] = { 3.7, 3.9, 4.2, 5.0 };
+
+//_________________________________________________________________________
+AliTRDCalPIDLQ::AliTRDCalPIDLQ()
+ :AliTRDCalPID("pid", "LQ PID references for TRD")
+{
+ //
+ // The Default constructor
+ //
+
+ Init();
+
+}
+
+//_________________________________________________________________________
+AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
+ :AliTRDCalPID(name,title)
+{
+ //
+ // The main constructor
+ //
+
+ Init();
+
+}
+
+//_________________________________________________________________________
+AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
+{
+ //
+ // Destructor
+ //
+
+}
+
+//_________________________________________________________________________
+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));
+ }
+
+
+ // 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;
+ }
+ gROOT->cd();
+
+ // 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));
+
+// 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;
+
+
+}
+
+//_________________________________________________________________________
+TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const // iType not needed
+{
+ //
+ // Returns one selected dEdx histogram
+ //
+
+ 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[iType], fTrackMomentum[ip]));
+
+ return fModel->At(GetModelID(ip, iType, iplane));
+}
+
+//_________________________________________________________________________
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, 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
+ //
+
+ 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*)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));
+ //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));
+ //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
+ if(mom < fTrackMomentum[0]) return LQ1;
+ else if(mom > fTrackMomentum[kNMom-1]) return LQ2;
+ else return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
+
+}
+
+//_________________________________________________________________________
+void AliTRDCalPIDLQ::Init()
+{
+ //
+ // Initialization
+ //
+
+ fModel = new TObjArray(AliPID::kSPECIES * kNMom);
+ fModel -> SetOwner();
+
+}
+
+//_________________________________________________________________________
+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)
+
+ return spec * AliTRDCalPID::kNMom + mom;
+}
--- /dev/null
+#ifndef ALITRDCALPIDLQ_H
+#define ALITRDCALPIDLQ_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// PID distributions for the LQ method //
+// //
+// Author: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDCALPID_H
+#include "AliTRDCalPID.h"
+#endif
+
+class AliTRDCalPIDLQ : public AliTRDCalPID
+{
+
+ public:
+
+ enum {
+ kNLength = 4
+ };
+
+ AliTRDCalPIDLQ();
+ AliTRDCalPIDLQ(const Text_t *name, const Text_t *title);
+ virtual ~AliTRDCalPIDLQ();
+
+ 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. : fTrackSegLength[il]; }
+ Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx
+ , Float_t length, Int_t plane) const;
+
+ private:
+
+ AliTRDCalPIDLQ(const AliTRDCalPIDLQ& pd);
+ AliTRDCalPIDLQ& operator=(const AliTRDCalPIDLQ &c);
+
+ void Init();
+ Int_t GetModelID(Int_t mom, Int_t , Int_t) const;
+
+ protected:
+
+ static Float_t fTrackSegLength[kNLength]; // Track segment lengths
+
+ ClassDef(AliTRDCalPIDLQ, 1) // LQ PID reference manager
+
+};
+#endif
+
--- /dev/null
+/**************************************************************************
+ * 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 of the distributions for the neural network //
+// //
+// Author: //
+// Alex Wilk <wilka@uni-muenster.de> //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#include <TFile.h>
+#include <TROOT.h>
+#include <TObject.h>
+#include <TMultiLayerPerceptron.h>
+
+#include "AliLog.h"
+#include "AliPID.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliTRDtrack.h"
+
+#include "AliTRDCalPIDNN.h"
+#include "AliTRDcalibDB.h"
+
+ClassImp(AliTRDCalPIDNN)
+
+//_________________________________________________________________________
+AliTRDCalPIDNN::AliTRDCalPIDNN()
+ :AliTRDCalPID("pid", "NN PID references for TRD")
+{
+ //
+ // The Default constructor
+ //
+
+ Init();
+
+}
+
+//_________________________________________________________________________
+AliTRDCalPIDNN::AliTRDCalPIDNN(const Text_t *name, const Text_t *title)
+ :AliTRDCalPID(name,title)
+{
+ //
+ // The main constructor
+ //
+
+ Init();
+
+}
+
+// //_____________________________________________________________________________
+// AliTRDCalPIDNN::AliTRDCalPIDNN(const AliTRDCalPIDNN &c)
+// :TNamed(c)
+// ,fMeanChargeRatio(c.fMeanChargeRatio)
+// ,fModel(0x0)
+// {
+// //
+// // Copy constructor
+// //
+//
+// if (this != &c) ((AliTRDCalPIDNN &) c).Copy(*this);
+//
+// }
+
+//_________________________________________________________________________
+AliTRDCalPIDNN::~AliTRDCalPIDNN()
+{
+ //
+ // Destructor
+ //
+
+}
+
+//_________________________________________________________________________
+Bool_t AliTRDCalPIDNN::LoadReferences(Char_t *refFile)
+{
+ //
+ // Read the TRD Neural Networks
+ //
+
+ // Read NN Root file
+ TFile *nnFile = TFile::Open(refFile,"READ");
+ if (!nnFile || !nnFile->IsOpen()) {
+ AliError(Form("Opening TRD histgram file %s failed",refFile));
+ return kFALSE;
+ }
+ gROOT->cd();
+
+ // Read Networks
+ for (Int_t imom = 0; imom < kNMom; imom++) {
+ for (Int_t iplane = 0; iplane < kNPlane; iplane++) {
+ TMultiLayerPerceptron *nn = (TMultiLayerPerceptron *)
+ nnFile->Get(Form("NN_Mom%d_Plane%d",imom,iplane));
+ fModel->AddAt(nn,GetModelID(imom,0,iplane));
+ }
+ }
+
+ nnFile->Close();
+ delete nnFile;
+
+ return kTRUE;
+
+}
+
+//_________________________________________________________________________
+TObject *AliTRDCalPIDNN::GetModel(Int_t ip, Int_t, Int_t iplane) const
+{
+ //
+ // Returns one selected TMultiLayerPerceptron. iType not used.
+ //
+
+ if (ip<0 || ip>= kNMom) return 0x0;
+
+ AliInfo(Form("Retrive MultiLayerPerceptron for %5.2f GeV/c for plane %d"
+ ,fTrackMomentum[ip]
+ ,iplane));
+
+ return fModel->At(GetModelID(ip, 0, iplane));
+
+}
+
+//_________________________________________________________________________
+Double_t AliTRDCalPIDNN::GetProbability(Int_t spec, Float_t mom, Float_t *dedx
+ , Float_t, Int_t iplane) const
+{
+ //
+ // Core function of AliTRDCalPID class for calculating the
+ // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
+ // given momentum "mom", a given dE/dx (slice "dedx") yield per
+ // layer in a given layer (iplane)
+ //
+
+ if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
+
+ // find the interval in momentum and track segment length which applies for this data
+
+ Int_t imom = 1;
+ while(imom<AliTRDCalPID::kNMom-1 && mom>fTrackMomentum[imom]) imom++;
+ Double_t LNN1, LNN2;
+ Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
+
+ TMultiLayerPerceptron *nn = 0x0;
+ if(!(nn = (TMultiLayerPerceptron *) fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
+ //if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
+ AliInfo(Form("Looking for mom(%f) plane(%d)", mom-1, iplane));
+ AliError(Form("NN id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
+ return 0.;
+ }
+
+ Double_t ddedx[AliTRDtrack::kNMLPslice];
+ for(int inode=0; inode<AliTRDtrack::kNMLPslice; inode++) ddedx[inode] = (Double_t)dedx[inode];
+ LNN1 = nn->Evaluate(spec, ddedx);
+
+ if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
+ //if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
+ AliInfo(Form("Looking for mom(%f) plane(%d)", mom, iplane));
+ AliError(Form("NN id %d not found in DB.", GetModelID(imom, spec, iplane)));
+ return LNN1;
+ }
+
+ LNN2 = nn->Evaluate(spec, ddedx);
+
+ // return interpolation over momentum binning
+ if(mom < fTrackMomentum[0]) return LNN1;
+ else if(mom > fTrackMomentum[AliTRDCalPID::kNMom-1]) return LNN2;
+ else return LNN1 + (LNN2 - LNN1)*(mom - mom1)/(mom2 - mom1);
+
+}
+
+//_________________________________________________________________________
+void AliTRDCalPIDNN::Init()
+{
+ //
+ // Initialization
+ //
+
+ fModel = new TObjArray(AliTRDCalPID::kNPlane * AliTRDCalPID::kNMom);
+ fModel->SetOwner();
+
+}
+
+//_________________________________________________________________________
+Int_t AliTRDCalPIDNN::GetModelID(Int_t mom, Int_t, Int_t plane) const
+{
+
+ // returns the ID of the NN distribution (66 MLPs, ID from 56 to 121)
+
+ return /*AliPID::kSPECIES * AliTRDCalPID::kNMom + */
+ plane * AliTRDCalPID::kNMom + mom;
+
+}
--- /dev/null
+#ifndef ALITRDCALPIDNN_H
+#define ALITRDCALPIDNN_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// PID distributions for the NN method //
+// //
+// Author: //
+// Alex Wilk <wilka@uni-muenster.de> //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDCALPID_H
+#include "AliTRDCalPID.h"
+#endif
+
+class AliTRDCalPIDNN : public AliTRDCalPID
+{
+
+ public:
+
+ AliTRDCalPIDNN();
+ AliTRDCalPIDNN(const Text_t *name, const Text_t *title);
+ virtual ~AliTRDCalPIDNN();
+ Bool_t LoadReferences(Char_t *refFile);
+ TObject *GetModel(Int_t ip, Int_t iType, Int_t iPlane) const;
+ Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length, Int_t plane) const;
+
+ private:
+
+ AliTRDCalPIDNN(const AliTRDCalPIDNN &pd);
+ AliTRDCalPIDNN &operator=(const AliTRDCalPIDNN &c);
+
+ void Init();
+ Int_t GetModelID(Int_t mom, Int_t , Int_t) const;
+
+ ClassDef(AliTRDCalPIDNN, 1) // NN PID reference manager
+
+};
+#endif
// Print statistics header
Int_t nPart[AliPID::kSPECIES], nTotPart;
printf("P[GeV/c] ");
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::fpartSymb[ispec]);
+ for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::GetPartSymb(ispec));
printf("\n-----------------------------------------------\n");
Float_t trackMomentum[AliTRDCalPID::kNMom];
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
// check PCA data
if(!fPrinc[ispec]){
- AliError(Form("No data defined for %s.", AliTRDCalPID::fpartName[ispec]));
+ AliError(Form("No data defined for %s.", AliTRDCalPID::GetPartName(ispec)));
return;
}
// build reference histograms
xmin = ymin = 0.;
xmax = 8000.; ymax = 6000.;
if(!fH2dEdx[ispec]){
- fH2dEdx[ispec] = new TH2D(Form("h2%s", AliTRDCalPID::fpartSymb[ispec]), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+ fH2dEdx[ispec] = new TH2D(Form("h2%s", AliTRDCalPID::GetPartSymb(ispec)), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
fH2dEdx[ispec]->SetLineColor(color[ispec]);
}
}
// save dE/dx references
TH2 *h2 = 0x0;
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
- h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::fpartSymb[ispec], mom));
- h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::fpartName[ispec], mom));
+ h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::GetPartSymb(ispec), mom));
+ h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::GetPartName(ispec), mom));
h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
h2->GetZaxis()->SetTitle("Entries");
}
//___________________________________________________________________
-void generatePIDDB(const char *file = "Refs.root")
+void generatePIDDB(const char *fileNN = "NN.root", const char *fileLQ = "LQ.root")
{
// Write TRD PID DB using the reference data from file "file"
man->SetRun(0);
AliCDBStorage *gStorLoc = man->GetStorage("local://$ALICE_ROOT");
- if (!gStorLoc) return;
+ if (!gStorLoc) return;
- AliTRDCalPID *pid = new AliTRDCalPID("pid", "TRD PID object");
- pid->LoadLQReferences(file);
-
- 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(pid, AliCDBId("TRD/Calib/PIDLQ", 0, 0), md);
+ 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);
}
//___________________________________________________________________
-AliTRDCalPID* getPIDObject()
+AliTRDCalPID* getPIDObject(const char *method="NN")
{
// Returns PIDLQ object.
// In order to browse histos do:
CDBManager->SetDefaultStorage("local://$ALICE_ROOT");
CDBManager->SetRun(0);
- AliCDBEntry *wrap = CDBManager->Get("TRD/Calib/PIDLQ", 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){
#pragma link C++ class AliTRDCalDet+;
#pragma link C++ class AliTRDCalFEE+;
#pragma link C++ class AliTRDCalPID+;
+#pragma link C++ class AliTRDCalPIDLQ+;
+#pragma link C++ class AliTRDCalPIDNN+;
#pragma link C++ class AliTRDCalPIDRefMaker+;
#pragma link C++ class AliTRDCalMonitoring+;
Cal/AliTRDCalDet.cxx \
Cal/AliTRDCalFEE.cxx \
Cal/AliTRDCalPID.cxx \
+ Cal/AliTRDCalPIDLQ.cxx \
+ Cal/AliTRDCalPIDNN.cxx \
Cal/AliTRDCalPIDRefMaker.cxx \
Cal/AliTRDCalMonitoring.cxx \
Cal/AliTRDCalChamberStatus.cxx \
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff