// A. Bercuci - used for PID calculations
#include "AliTRDcalibDB.h"
-#include "Cal/AliTRDCalPIDLQ.h"
+#include "Cal/AliTRDCalPID.h"
ClassImp(AliTRDtrack)
,fSeedLab(-1)
,fdEdx(0)
,fDE(0)
- ,fClusterOwner(kFALSE) // A.Bercuci
+ ,fClusterOwner(kFALSE)
,fStopped(kFALSE)
,fLhElectron(0)
,fNWrong(0)
fdEdxPlane[i][j] = 0.0;
}
fTimBinPlane[i] = -1;
- // A.Bercuci additions
- fMom[i] = -1.;
- fSnp[i] = 0.;
- fTgl[i] = 0.;
- }
+ fMom[i] = -1.;
+ fSnp[i] = 0.;
+ fTgl[i] = 0.;
+ }
for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
fIndex[i] = 0;
fIndexBackup[i] = 0;
fdQdl[i] = 0;
- //A.Bercuci additions
- fClusters[i] = 0x0;
- }
+ fClusters[i] = 0x0;
+ }
for (Int_t i = 0; i < 3; i++) {
fBudget[i] = 0;
,fSeedLab(-1)
,fdEdx(0)
,fDE(0)
- ,fClusterOwner(kFALSE) // A.Bercuci
+ ,fClusterOwner(kFALSE)
,fStopped(kFALSE)
,fLhElectron(0)
,fNWrong(0)
// The main AliTRDtrack constructor.
//
- Double_t cnv = 1.0/(GetBz() * kB2C);
+ Double_t cnv = 1.0 / (GetBz() * kB2C);
Double_t pp[5] = { p[0]
, p[1]
, p[3]
, p[4]*cnv };
- Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
- Double_t c32 = x*cov[13] - cov[8];
- Double_t c20 = x*cov[10] - cov[3];
- Double_t c21 = x*cov[11] - cov[4];
- Double_t c42 = x*cov[14] - cov[12];
+ Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
+ Double_t c32 = x*cov[13] - cov[ 8];
+ Double_t c20 = x*cov[10] - cov[ 3];
+ Double_t c21 = x*cov[11] - cov[ 4];
+ Double_t c42 = x*cov[14] - cov[12];
- Double_t cc[15] = { cov[0 ]
- , cov[1 ], cov[2 ]
+ Double_t cc[15] = { cov[ 0]
+ , cov[ 1], cov[ 2]
, c20, c21, c22
- , cov[6 ], cov[7 ], c32, cov[9 ]
+ , cov[ 6], cov[ 7], c32, cov[ 9]
, cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
Set(x,alpha,pp,cc);
SetNumberOfClusters(1);
- fIndex[0] = index;
- fClusters[0] = c; // A.Bercuci additions
+ fIndex[0] = index;
+ fClusters[0] = c;
for (Int_t i = 0; i < kNplane; i++) {
for (Int_t j = 0; j < kNslice; j++) {
fdEdxPlane[i][j] = 0.0;
}
fTimBinPlane[i] = -1;
- // A.Bercuci additions
- fMom[i] = -1.;
- fSnp[i] = 0.;
- fTgl[i] = 0.;
+ fMom[i] = -1.;
+ fSnp[i] = 0.;
+ fTgl[i] = 0.;
}
Double_t q = TMath::Abs(c->GetQ());
q *= TMath::Sqrt((1-s*s)/(1+t*t));
}
- fdQdl[0] = q;
+ fdQdl[0] = q;
for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
fdQdl[i] = 0;
fIndex[i] = 0;
fIndexBackup[i] = 0;
- // A.Bercuci additions
- fClusters[i] = 0x0;
+ fClusters[i] = 0x0;
}
for (Int_t i = 0; i < 3;i++) {
,fSeedLab(t.GetSeedLabel())
,fdEdx(t.fdEdx)
,fDE(t.fDE)
- ,fClusterOwner(kTRUE) // A.Bercuci
+ ,fClusterOwner(kTRUE)
,fStopped(t.fStopped)
,fLhElectron(0)
,fNWrong(t.fNWrong)
}
fTimBinPlane[i] = t.fTimBinPlane[i];
fTracklets[i] = t.fTracklets[i];
- // A.Bercuci additions
- fMom[i] = t.fMom[i];
- fSnp[i] = t.fSnp[i];
- fTgl[i] = t.fTgl[i];
+ fMom[i] = t.fMom[i];
+ fSnp[i] = t.fSnp[i];
+ fTgl[i] = t.fTgl[i];
}
Int_t n = t.GetNumberOfClusters();
fIndex[i] = t.fIndex[i];
fIndexBackup[i] = t.fIndex[i];
fdQdl[i] = t.fdQdl[i];
- // A.Bercuci additions
- if(fClusterOwner && t.fClusters[i]) fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
- else fClusters[i] = t.fClusters[i];
- }
+ if (fClusterOwner && t.fClusters[i]) {
+ fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
+ }
+ else {
+ fClusters[i] = t.fClusters[i];
+ }
+ }
for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
fdQdl[i] = 0;
fIndex[i] = 0;
fIndexBackup[i] = 0;
- // A.Bercuci additions
- fClusters[i] = 0x0;
+ fClusters[i] = 0x0;
}
for (Int_t i = 0; i < 3;i++) {
,fSeedLab(-1)
,fdEdx(t.GetPIDsignal())
,fDE(0)
- ,fClusterOwner(kFALSE) // A.Bercuci
+ ,fClusterOwner(kFALSE)
,fStopped(kFALSE)
,fLhElectron(0.0)
,fNWrong(0)
fdEdxPlane[i][j] = 0.0;
}
fTimBinPlane[i] = -1;
- // A.Bercuci additions
- fMom[i] = -1.;
- fSnp[i] = 0.;
- fTgl[i] = 0.;
+ fMom[i] = -1.;
+ fSnp[i] = 0.;
+ fTgl[i] = 0.;
}
for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
fdQdl[i] = 0;
fIndex[i] = 0;
fIndexBackup[i] = 0;
- // A.Bercuci additions
- fClusters[i] = 0x0;
+ fClusters[i] = 0x0;
}
for (Int_t i = 0; i < 3; i++) {
,fSeedLab(-1)
,fdEdx(t.GetTRDsignal())
,fDE(0)
- ,fClusterOwner(kFALSE) // A.Bercuci
+ ,fClusterOwner(kFALSE)
,fStopped(kFALSE)
,fLhElectron(0)
,fNWrong(0)
fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
}
fTimBinPlane[i] = t.GetTRDTimBin(i);
- // A.Bercuci additions
- fMom[i] = -1.;
- fSnp[i] = 0.;
- fTgl[i] = 0.;
+ fMom[i] = -1.;
+ fSnp[i] = 0.;
+ fTgl[i] = 0.;
}
const AliExternalTrackParam *par = &t;
}
Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
-
for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
fdQdl[i] = 0;
- // A.Bercuci additions
- fClusters[i] = 0x0;
- }
+ fClusters[i] = 0x0;
+ }
for (Int_t i = 0; i < 3; i++) {
fBudget[i] = 0;
delete fBackupTrack;
}
fBackupTrack = 0x0;
- if (fClusterOwner){
- UInt_t icluster=0;
- while(icluster<kMAXCLUSTERSPERTRACK && fClusters[icluster]){
- delete fClusters[icluster];
- fClusters[icluster] = 0x0;
- icluster++;
- }
+ if (fClusterOwner) {
+ UInt_t icluster = 0;
+ while ((icluster < kMAXCLUSTERSPERTRACK) && fClusters[icluster]) {
+ delete fClusters[icluster];
+ fClusters[icluster] = 0x0;
+ icluster++;
+ }
}
}
// This part of the function is never reached ????
// What defines these parameters ????
- Int_t status = 0;
- if (GetNumberOfClusters() < 20) return 0;
- if ((fN > 110) &&
- (fChi2/(Float_t(fN)) < 3)) return 3; // Gold
- if ((fNLast > 30) &&
- (fChi2Last/(Float_t(fNLast)) < 3)) return 3; // Gold
- if ((fNLast > 20) &&
- (fChi2Last/(Float_t(fNLast)) < 2)) return 3; // Gold
- if ((fNLast/(fNExpectedLast+3.0) > 0.8) &&
- (fChi2Last/Float_t(fNLast) < 5) &&
- (fNLast > 20)) return 2; // Silber
- if ((fNLast > 5) &&
- (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) &&
- (fChi2Last/(fNLast-5.0) < 6)) return 1;
-
- return status;
+ //Int_t status = 0;
+ //if (GetNumberOfClusters() < 20) return 0;
+ //if ((fN > 110) &&
+ // (fChi2/(Float_t(fN)) < 3)) return 3; // Gold
+ //if ((fNLast > 30) &&
+ // (fChi2Last/(Float_t(fNLast)) < 3)) return 3; // Gold
+ //if ((fNLast > 20) &&
+ // (fChi2Last/(Float_t(fNLast)) < 2)) return 3; // Gold
+ //if ((fNLast/(fNExpectedLast+3.0) > 0.8) &&
+ // (fChi2Last/Float_t(fNLast) < 5) &&
+ // (fNLast > 20)) return 2; // Silver
+ //if ((fNLast > 5) &&
+ // (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) &&
+ // (fChi2Last/(fNLast-5.0) < 6)) return 1;
+ //
+ //return status;
}
//
// Array to sort the dEdx values according to amplitude
-
- Float_t sorted[kMAXCLUSTERSPERTRACK];
- fdEdx = 0.;
+ Float_t sorted[kMAXCLUSTERSPERTRACK];
+ fdEdx = 0.;
// Require at least 10 clusters for a dedx measurement
if (fNdedx < 10) return;
-
// Can fdQdl be negative ????
- for (Int_t i = 0; i < fNdedx; i++) sorted[i] = TMath::Abs(fdQdl[i]);
+ for (Int_t i = 0; i < fNdedx; i++) {
+ sorted[i] = TMath::Abs(fdQdl[i]);
+ }
// Sort the dedx values by amplitude
Int_t *index = new Int_t[fNdedx];
TMath::Sort(fNdedx, sorted, index, kFALSE);
// Sum up the truncated charge between lower and upper bounds
Int_t nl = Int_t(low * fNdedx);
Int_t nu = Int_t( up * fNdedx);
- for (Int_t i = nl; i <= nu; i++) fdEdx += sorted[index[i]];
- fdEdx /= (nu - nl + 1.0);
+ for (Int_t i = nl; i <= nu; i++) {
+ fdEdx += sorted[index[i]];
+ }
+ fdEdx /= (nu - nl + 1.0);
+
+ delete[] index;
- delete[] index;
}
//_____________________________________________________________________________
-void AliTRDtrack::CookdEdxTimBin()
+void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
{
//
// Set fdEdxPlane and fTimBinPlane and also get the
// Time bin for Max. Cluster
- //
- // Authors:
+ //
+ // Authors:
// Prashant Shukla (shukla@physi.uni-heidelberg.de)
// Alexandru Bercuci (A.Bercuci@gsi.de)
+ //
-
- Double_t maxcharge[AliESDtrack::kNPlane]; // max charge in chamber
- // number of clusters attached to track per chamber and slice
- Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
- //number of time bins in chamber
- Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
- 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
-
- // Reset class and local contors/variables
+ Double_t maxcharge[AliESDtrack::kNPlane]; // max charge in chamber
+ // Number of clusters attached to track per chamber and slice
+ Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
+ // Number of time bins in chamber
+ Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+ Int_t plane; // Plane of current cluster
+ Int_t tb; // Time bin of current cluster
+ Int_t slice; // Current slice
+ AliTRDcluster *cluster = 0x0; // Pointer to current cluster
+
+ // Reset class and local counters/variables
for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
- fTimBinPlane[iPlane] = -1;
- maxcharge[iPlane] = 0.;
- for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+ fTimBinPlane[iPlane] = -1;
+ maxcharge[iPlane] = 0.;
+ for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
fdEdxPlane[iPlane][iSlice] = 0.;
nCluster[iPlane][iSlice] = 0;
}
}
- // start looping over clusters attached to this track
- for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
+ // 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());
+ // 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;
- }
+ 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 * AliESDtrack::kNSlice / ntb;
-
- fdEdxPlane[plane][slice] += fdQdl[iClus];
- if(fdQdl[iClus] > maxcharge[plane]) {
- maxcharge[plane] = fdQdl[iClus];
+ slice = tb * AliESDtrack::kNSlice / ntb;
+
+ fdEdxPlane[plane][slice] += fdQdl[iClus];
+ if (fdQdl[iClus] > maxcharge[plane]) {
+ maxcharge[plane] = fdQdl[iClus];
fTimBinPlane[plane] = tb;
}
- nCluster[plane][slice]++;
- } // End of loop over cluster
+ nCluster[plane][slice]++;
+ } // End of loop over cluster
// Normalize fdEdxPlane to number of clusters and set track segments
for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
- if (nCluster[iPlane][iSlice]) fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
+ if (nCluster[iPlane][iSlice]) {
+ fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
+ }
}
+ }
+
+}
+
+//_____________________________________________________________________________
+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)
{
- if(plane<0 || plane>= kNplane){
- AliError(Form("Trying to access out of range plane (%d)", plane));
- return;
- }
+ //
+ // Set the momenta for a track segment in a given plane
+ //
+
+ if ((plane < 0) ||
+ (plane>= kNplane)) {
+ AliError(Form("Trying to access out of range plane (%d)", plane));
+ return;
+ }
- fSnp[plane] = GetSnp();
- fTgl[plane] = GetTgl();
- Double_t p[3]; GetPxPyPz(p);
- fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
+ fSnp[plane] = GetSnp();
+ fTgl[plane] = GetTgl();
+ Double_t p[3];
+ GetPxPyPz(p);
+ fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
+
}
//_____________________________________________________________________________
Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
{
- if(plane < 0 || plane >= kNplane) return 0.;
- return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1.
-- fSnp[plane]*fSnp[plane]) / (1. + fTgl[plane]*fTgl[plane]));
+ //
+ // Calculate the track length of a track segment in a given plane
+ //
-}
+ if ((plane < 0) || (plane >= kNplane)) return 0.;
+ return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
+ / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane])
+ / (1.0 + fTgl[plane]*fTgl[plane]));
+
+}
//_____________________________________________________________________________
-Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
+Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
{
//
// This function calculates the PID probabilities based on TRD signals
// The method produces probabilities based on the charge
// and the position of the maximum time bin in each layer.
// The dE/dx information can be used as global charge or 2 to 3
- // slices. Check AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual
+ // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
// implementation.
//
// Author
// Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
+ //
-
- AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
- if (!calibration) {
- AliError("No access to calibration data");
- return -1;
- }
+ AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+ if (!calibration) {
+ AliError("No access to calibration data");
+ return kFALSE;
+ }
- // Retrieve the CDB container class with the probability distributions
- const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
- if (!pd) {
- AliError("No access to AliTRDCalPIDLQ");
- 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 x0, Double_t rho)
+Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
{
//
- // Propagates a track of particle with mass=pm to a reference plane
- // defined by x=xk through media of density=rho and radiationLength=x0
+ // Propagates this track to a reference plane defined by "xk" [cm]
+ // correcting for the mean crossed material.
//
+ // "xx0" - thickness/rad.length [units of the radiation length]
+ // "xrho" - thickness*density [g/cm^2]
+ //
if (xk == GetX()) {
return kTRUE;
Double_t x = GetX();
Double_t y = GetY();
Double_t z = GetZ();
- Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
- + (y-oldY)*(y-oldY)
- + (z-oldZ)*(z-oldZ));
+
if (oldX < xk) {
+ xrho = -xrho;
if (IsStartedTimeIntegral()) {
- Double_t l2 = d;
+ Double_t l2 = TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ));
Double_t crv = GetC();
if (TMath::Abs(l2*crv) > 0.0001) {
// Make correction for curvature if neccesary
}
}
- Double_t ll = (oldX < xk) ? -d : d;
- if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
+ if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
return kFALSE;
}
{
// Energy losses************************
- Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
+ Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
return kFALSE;
Double_t dE = 0.153e-3 / beta2
* (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
- * d * rho;
- Float_t budget = d * rho;
- fBudget[0] += budget;
+ * xrho;
+ fBudget[0] += xrho;
/*
// Suspicious part - think about it ?
//Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
Double_t sy2 = c->GetSigmaY2() * 4.0;
Double_t sz2 = c->GetSigmaZ2() * 4.0;
- Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
+ Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
if (!AliExternalTrackParam::Update(p,cov)) {
return kFALSE;
}
- Int_t n = GetNumberOfClusters();
+ Int_t n = GetNumberOfClusters();
fIndex[n] = index;
SetNumberOfClusters(n+1);
-
-
+
SetChi2(GetChi2()+chisq);
- return kTRUE;
+ return kTRUE;
}
//_____________________________________________________________________________
-Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/)
+Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
+ , Double_t h01, Int_t /*plane*/, Int_t tid)
{
//
// Assignes found cluster to the track and updates track information
return kFALSE;
}
- // Register cluster to track
- Int_t n = GetNumberOfClusters();
- fIndex[n] = index;
- fClusters[n] = c; // A.Bercuci 25.07.07
+ // Register cluster to track
+ Int_t n = GetNumberOfClusters();
+ fIndex[n] = index;
+ fClusters[n] = c;
SetNumberOfClusters(n+1);
SetChi2(GetChi2() + chisq);
// fCtt-=k30*oldty+k31*oldtz;
// fCct-=k40*oldty+k41*oldtz;
// //
+
// fCcc-=k40*oldcy+k41*oldcz;
// //
{
//
// Propagate track to given x position
- // Works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
+ // Works inside of the 20 degree segmentation
+ // (local cooordinate frame for TRD , TPC, TOF)
//
// Material budget from geo manager
//
const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
// Critical alpha - cross sector indication
- Double_t dir = (GetX()>xr) ? -1.0 : 1.0;
+ Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
// Direction +-
for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
xyz1[2] = z;
Double_t param[7];
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+ AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
if ((param[0] > 0) &&
(param[1] > 0)) {
- PropagateTo(x,param[1],param[0]);
+ PropagateTo(x,param[1],param[0]*param[4]);
}
if (GetY() > GetX()*kTalphac) {
Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
// Direction +-
- Double_t dir = (radius>r) ? -1.0 : 1.0;
+ Double_t dir = (radius > r) ? -1.0 : 1.0;
for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
xyz1[2] = z;
Double_t param[7];
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+ AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
if (param[1] <= 0) {
param[1] = 100000000;
}
- PropagateTo(x,param[1],param[0]);
+ PropagateTo(x,param[1],param[0]*param[4]);
}
xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
xyz1[2] = z;
Double_t param[7];
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+ AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
if (param[1] <= 0) {
param[1] = 100000000;
}
- PropagateTo(r,param[1],param[0]);
+ PropagateTo(r,param[1],param[0]*param[4]);
return 0;
//
return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
- % AliTRDgeometry::kNsect;
+ % AliTRDgeometry::kNsect;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff