#include "AliTRDchamberTimeBin.h"
#include "AliTRDtrackingChamber.h"
#include "AliTRDtrackerV1.h"
-#include "AliTRDReconstructor.h"
#include "AliTRDrecoParam.h"
#include "AliTRDCommonParam.h"
ClassImp(AliTRDseedV1)
+TLinearFitter *AliTRDseedV1::fgFitterY = NULL;
+TLinearFitter *AliTRDseedV1::fgFitterZ = NULL;
+
//____________________________________________________________________
AliTRDseedV1::AliTRDseedV1(Int_t det)
:AliTRDtrackletBase()
- ,fReconstructor(0x0)
- ,fClusterIter(0x0)
+ ,fkReconstructor(NULL)
+ ,fClusterIter(NULL)
,fExB(0.)
,fVD(0.)
,fT0(0.)
//
// Constructor
//
- for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
+ memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
memset(fPad, 0, 3*sizeof(Float_t));
fYref[0] = 0.; fYref[1] = 0.;
//____________________________________________________________________
AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
:AliTRDtrackletBase((AliTRDtrackletBase&)ref)
- ,fReconstructor(0x0)
- ,fClusterIter(0x0)
+ ,fkReconstructor(NULL)
+ ,fClusterIter(NULL)
,fExB(0.)
,fVD(0.)
,fT0(0.)
if(!fClusters[itb]) continue;
//AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
delete fClusters[itb];
- fClusters[itb] = 0x0;
+ fClusters[itb] = NULL;
}
}
}
//AliInfo("");
AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
- target.fReconstructor = fReconstructor;
- target.fClusterIter = 0x0;
+ target.fkReconstructor = fkReconstructor;
+ target.fClusterIter = NULL;
target.fExB = fExB;
target.fVD = fVD;
target.fT0 = fT0;
if((*c)->IsShared() || (*c)->IsUsed()){
if((*c)->IsShared()) SetNShared(GetNShared()-1);
else SetNUsed(GetNUsed()-1);
- (*c) = 0x0;
+ (*c) = NULL;
fIndexes[ic] = -1;
SetN(GetN()-1);
continue;
const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
- AliTRDcluster *c = 0x0;
+ AliTRDcluster *c = NULL;
for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
Float_t dx = TMath::Abs(fX0 - c->GetX());
nclusters[slice]++;
} // End of loop over clusters
- //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
+ //if(fkReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
if(nslices == AliTRDpidUtil::kLQslices){
// calculate mean charge per slice (only LQ PID)
for(int is=0; is<nslices; is++){
Float_t* AliTRDseedV1::GetProbability(Bool_t force)
{
if(!force) return &fProb[0];
- if(!CookPID()) return 0x0;
+ if(!CookPID()) return NULL;
return &fProb[0];
}
// Parameters
//
// Output
-// returns pointer to the probability array and 0x0 if missing DB access
+// returns pointer to the probability array and NULL if missing DB access
//
-// Detailed description
+// Retrieve PID probabilities for e+-, mu+-, K+-, pi+- and p+- from the DB according to tracklet information:
+// - estimated momentum at tracklet reference point
+// - dE/dx measurements
+// - tracklet length
+// - TRD layer
+// According to the steering settings specified in the reconstruction one of the following methods are used
+// - Neural Network [default] - option "nn"
+// - 2D Likelihood - option "!nn"
-
- // retrive calibration db
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
AliError("No access to calibration data");
return kFALSE;
}
- if (!fReconstructor) {
+ if (!fkReconstructor) {
AliError("Reconstructor not set.");
return kFALSE;
}
// Retrieve the CDB container class with the parametric detector response
- const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
+ const AliTRDCalPID *pd = calibration->GetPIDObject(fkReconstructor->GetPIDMethod());
if (!pd) {
AliError("No access to AliTRDCalPID object");
return kFALSE;
}
- //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
+ //AliInfo(Form("Method[%d] : %s", fkReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
// calculate tracklet length TO DO
Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
/// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
//calculate dE/dx
- CookdEdx(fReconstructor->GetNdEdxSlices());
+ CookdEdx(fkReconstructor->GetNdEdxSlices());
// Sets the a priori probabilities
- for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
- fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane());
- }
-
+ for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
+ fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane());
+
return kTRUE;
}
//GetPadLength()*GetPadLength()/12.;
// insert systematic uncertainties
- if(fReconstructor){
+ if(fkReconstructor){
Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
- fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
+ fkReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
sy2 += sys[0];
sz2 += sys[1];
}
}
//____________________________________________________________
-Double_t AliTRDseedV1::GetCovSqrt(Double_t *c, Double_t *d)
+Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
{
// Helper function to calculate the square root of the covariance matrix.
// The input matrix is stored in the vector c and the result in the vector d.
// Author A.Bercuci <A.Bercuci@gsi.de>
// Date Mar 19 2009
- Double_t L[2], // eigenvalues
- V[3]; // eigenvectors
+ Double_t l[2], // eigenvalues
+ v[3]; // eigenvectors
// the secular equation and its solution :
// (c[0]-L)(c[2]-L)-c[1]^2 = 0
// L^2 - L*Tr(c)+DET(c) = 0
// L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
- Double_t Tr = c[0]+c[2], // trace
- DET = c[0]*c[2]-c[1]*c[1]; // determinant
- if(TMath::Abs(DET)<1.e-20) return -1.;
- Double_t DD = TMath::Sqrt(Tr*Tr - 4*DET);
- L[0] = .5*(Tr + DD);
- L[1] = .5*(Tr - DD);
- if(L[0]<0. || L[1]<0.) return -1.;
+ Double_t tr = c[0]+c[2], // trace
+ det = c[0]*c[2]-c[1]*c[1]; // determinant
+ if(TMath::Abs(det)<1.e-20) return -1.;
+ Double_t dd = TMath::Sqrt(tr*tr - 4*det);
+ l[0] = .5*(tr + dd);
+ l[1] = .5*(tr - dd);
+ if(l[0]<0. || l[1]<0.) return -1.;
// the sym V matrix
// | v00 v10|
// | v10 v11|
- Double_t tmp = (L[0]-c[0])/c[1];
- V[0] = TMath::Sqrt(1./(tmp*tmp+1));
- V[1] = tmp*V[0];
- V[2] = V[1]*c[1]/(L[1]-c[2]);
+ Double_t tmp = (l[0]-c[0])/c[1];
+ v[0] = TMath::Sqrt(1./(tmp*tmp+1));
+ v[1] = tmp*v[0];
+ v[2] = v[1]*c[1]/(l[1]-c[2]);
// the VD^{1/2}V is:
- L[0] = TMath::Sqrt(L[0]); L[1] = TMath::Sqrt(L[1]);
- d[0] = V[0]*V[0]*L[0]+V[1]*V[1]*L[1];
- d[1] = V[0]*V[1]*L[0]+V[1]*V[2]*L[1];
- d[2] = V[1]*V[1]*L[0]+V[2]*V[2]*L[1];
+ l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
+ d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
+ d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
+ d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
return 1.;
}
//____________________________________________________________
-Double_t AliTRDseedV1::GetCovInv(Double_t *c, Double_t *d)
+Double_t AliTRDseedV1::GetCovInv(const Double_t * const c, Double_t *d)
{
// Helper function to calculate the inverse of the covariance matrix.
// The input matrix is stored in the vector c and the result in the vector d.
// Author A.Bercuci <A.Bercuci@gsi.de>
// Date Mar 19 2009
- Double_t Det = c[0]*c[2] - c[1]*c[1];
- if(TMath::Abs(Det)<1.e-20) return 0.;
- Double_t InvDet = 1./Det;
- d[0] = c[2]*InvDet;
- d[1] =-c[1]*InvDet;
- d[2] = c[0]*InvDet;
- return Det;
+ Double_t det = c[0]*c[2] - c[1]*c[1];
+ if(TMath::Abs(det)<1.e-20) return 0.;
+ Double_t invDet = 1./det;
+ d[0] = c[2]*invDet;
+ d[1] =-c[1]*invDet;
+ d[2] = c[0]*invDet;
+ return det;
}
//____________________________________________________________________
return fClusters[ic] ? fClusters[ic]->GetVolumeId() : 0;
}
+//____________________________________________________________________
+TLinearFitter* AliTRDseedV1::GetFitterY()
+{
+ if(!fgFitterY) fgFitterY = new TLinearFitter(1, "pol1");
+ fgFitterY->ClearPoints();
+ return fgFitterY;
+}
+
+//____________________________________________________________________
+TLinearFitter* AliTRDseedV1::GetFitterZ()
+{
+ if(!fgFitterZ) fgFitterZ = new TLinearFitter(1, "pol1");
+ fgFitterZ->ClearPoints();
+ return fgFitterZ;
+}
//____________________________________________________________________
void AliTRDseedV1::Calibrate()
//____________________________________________________________________
-Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
+Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt)
{
//
// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
// r_{z} = 1.5*L_{pad}
// END_LATEX
//
-// Author Alexandru Bercuci <A.Bercuci@gsi.de>
+// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
+// Debug : level >3
Bool_t kPRINT = kFALSE;
- if(!fReconstructor->GetRecoParam() ){
+ if(!fkReconstructor->GetRecoParam() ){
AliError("Seed can not be used without a valid RecoParam.");
return kFALSE;
}
// Initialize reco params for this tracklet
// 1. first time bin in the drift region
Int_t t0 = 14;
- Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
+ Int_t kClmin = Int_t(fkReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
+ Double_t sysCov[5]; fkReconstructor->GetRecoParam()->GetSysCovMatrix(sysCov);
Double_t s2yTrk= fRefCov[0],
s2yCl = 0.,
s2zCl = GetPadLength()*GetPadLength()/12.,
syRef = TMath::Sqrt(s2yTrk),
t2 = GetTilt()*GetTilt();
//define roads
- Double_t kroady = 1., //fReconstructor->GetRecoParam() ->GetRoad1y();
- kroadz = GetPadLength() * 1.5 + 1.;
+ Double_t kroady = 1., //fkReconstructor->GetRecoParam() ->GetRoad1y();
+ kroadz = GetPadLength() * fkReconstructor->GetRecoParam()->GetRoadzMultiplicator() + 1.;
// define probing cluster (the perfect cluster) and default calibration
Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
+ if(fkReconstructor->IsHLT())cp.SetRPhiMethod(AliTRDcluster::kCOG);
Calibrate();
if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
// working variables
const Int_t kNrows = 16;
- AliTRDcluster *clst[kNrows][kNclusters];
- Double_t cond[4], dx, dy, yt, zt,
- yres[kNrows][kNclusters];
- Int_t idxs[kNrows][kNclusters], ncl[kNrows], ncls = 0;
+ const Int_t kNcls = 3*kNclusters; // buffer size
+ AliTRDcluster *clst[kNrows][kNcls];
+ Double_t cond[4], dx, dy, yt, zt, yres[kNrows][kNcls];
+ Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
memset(ncl, 0, kNrows*sizeof(Int_t));
- memset(clst, 0, kNrows*kNclusters*sizeof(AliTRDcluster*));
+ memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
+ memset(clst, 0, kNrows*kNcls*sizeof(AliTRDcluster*));
// Do cluster projection
- AliTRDcluster *c = 0x0;
- AliTRDchamberTimeBin *layer = 0x0;
+ AliTRDcluster *c = NULL;
+ AliTRDchamberTimeBin *layer = NULL;
Bool_t kBUFFER = kFALSE;
- for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
+ for (Int_t it = 0; it < kNtb; it++) {
if(!(layer = chamber->GetTB(it))) continue;
if(!Int_t(*layer)) continue;
// get track projection at layers position
// get standard cluster error corrected for tilt
cp.SetLocalTimeBin(it);
cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1./*zt*/, fYref[1]);
- s2yCl = (cp.GetSigmaY2() + t2*s2zCl)/(1.+t2);
+ s2yCl = (cp.GetSigmaY2() + sysCov[0] + t2*s2zCl)/(1.+t2);
// get estimated road
kroady = 3.*TMath::Sqrt(12.*(s2yTrk + s2yCl));
yres[r][ncl[r]] = dy;
ncl[r]++; ncls++;
- if(ncl[r] >= kNclusters) {
- AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kNclusters));
+ if(ncl[r] >= kNcls) {
+ AliWarning(Form("Cluster candidates reached buffer limit %d. Some may be lost.", kNcls));
kBUFFER = kTRUE;
break;
}
continue;
}
+ if(fkReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 3){
+ TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
+ TVectorD vdy(ncl[ir], yres[ir]);
+ UChar_t stat(0);
+ if(IsKink()) SETBIT(stat, 0);
+ if(IsStandAlone()) SETBIT(stat, 1);
+ cstreamer << "AttachClusters"
+ << "stat=" << stat
+ << "det=" << fDet
+ << "pt=" << fPt
+ << "s2y=" << s2yTrk
+ << "dy=" << &vdy
+ << "m=" << mean
+ << "s=" << syDis
+ << "\n";
+ }
+
// TODO check mean and sigma agains cluster resolution !!
if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syDis), syDis);
// select clusters on a 3 sigmaDistr level
Bool_t kFOUND = kFALSE;
for(Int_t ic = ncl[ir]; ic--;){
if(yres[ir][ic] - mean > 3. * syDis){
- clst[ir][ic] = 0x0; continue;
+ clst[ir][ic] = NULL; continue;
}
nrow[nr]++; kFOUND = kTRUE;
}
//
// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
//
- fReconstructor = rec;
+ fkReconstructor = rec;
AliTRDgeometry g;
AliTRDpadPlane *pp = g.GetPadPlane(fDet);
fPad[0] = pp->GetLengthIPad();
Double_t yt, zt;
//AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
- TLinearFitter fitterY(1, "pol1");
- TLinearFitter fitterZ(1, "pol1");
-
+ TLinearFitter& fitterY=*GetFitterY();
+ TLinearFitter& fitterZ=*GetFitterZ();
+
// book cluster information
Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
Int_t n = 0;
- AliTRDcluster *c=0x0, **jc = &fClusters[0];
+ AliTRDcluster *c=NULL, **jc = &fClusters[0];
for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
xc[ic] = -1.;
yc[ic] = 999.;
c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
sy[n] = TMath::Sqrt(c->GetSigmaY2());
- yc[n] = fReconstructor->UseGAUS() ?
+ yc[n] = fkReconstructor->UseGAUS() ?
c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
zc[n] = c->GetZ();
//optional tilt correction
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff