////////////////////////////////////////////////////////////////////////////
#include "TMath.h"
-#include "TLinearFitter.h"
-#include "TClonesArray.h" // tmp
#include <TTreeStream.h>
#include "AliLog.h"
#include "AliMathBase.h"
+#include "AliRieman.h"
#include "AliCDBManager.h"
-#include "AliTracker.h"
+#include "AliTRDReconstructor.h"
#include "AliTRDpadPlane.h"
#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "Cal/AliTRDCalROC.h"
#include "Cal/AliTRDCalDet.h"
-ClassImp(AliTRDseedV1)
+class AliTracker;
-TLinearFitter *AliTRDseedV1::fgFitterY = NULL;
-TLinearFitter *AliTRDseedV1::fgFitterZ = NULL;
+ClassImp(AliTRDseedV1)
//____________________________________________________________________
AliTRDseedV1::AliTRDseedV1(Int_t det)
,fZ(0.)
,fS2Y(0.)
,fS2Z(0.)
- ,fC(0.)
,fChi2(0.)
{
//
//
memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
- memset(fPad, 0, 3*sizeof(Float_t));
+ memset(fPad, 0, 4*sizeof(Float_t));
fYref[0] = 0.; fYref[1] = 0.;
fZref[0] = 0.; fZref[1] = 0.;
fYfit[0] = 0.; fYfit[1] = 0.;
fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
fLabels[2]=0; // number of different labels for tracklet
memset(fRefCov, 0, 7*sizeof(Double_t));
+ // stand alone curvature
+ fC[0] = 0.; fC[1] = 0.;
// covariance matrix [diagonal]
// default sy = 200um and sz = 2.3 cm
fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
,fZ(0.)
,fS2Y(0.)
,fS2Z(0.)
- ,fC(0.)
,fChi2(0.)
{
//
target.fZ = fZ;
target.fS2Y = fS2Y;
target.fS2Z = fS2Z;
- target.fC = fC;
target.fChi2 = fChi2;
memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
- memcpy(target.fPad, fPad, 3*sizeof(Float_t));
+ memcpy(target.fPad, fPad, 4*sizeof(Float_t));
target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t));
+ target.fC[0] = fC[0]; target.fC[1] = fC[1];
memcpy(target.fCov, fCov, 3*sizeof(Double_t));
TObject::Copy(ref);
}
+//____________________________________________________________
+void AliTRDseedV1::Init(const AliRieman *rieman)
+{
+// Initialize this tracklet using the riemann fit information
+
+
+ fZref[0] = rieman->GetZat(fX0);
+ fZref[1] = rieman->GetDZat(fX0);
+ fYref[0] = rieman->GetYat(fX0);
+ fYref[1] = rieman->GetDYat(fX0);
+ if(fkReconstructor && fkReconstructor->IsHLT()){
+ fRefCov[0] = 1;
+ fRefCov[2] = 10;
+ }else{
+ fRefCov[0] = rieman->GetErrY(fX0);
+ fRefCov[2] = rieman->GetErrZ(fX0);
+ }
+ fC[0] = rieman->GetC();
+ fChi2 = rieman->GetChi2();
+}
+
+
//____________________________________________________________
Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
{
fPt=0.;
fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
fS2Y=0.; fS2Z=0.;
- fC=0.; fChi2 = 0.;
+ fC[0]=0.; fC[1]=0.;
+ fChi2 = 0.;
- memset(fPad, 0, 3*sizeof(Float_t));
+ memset(fPad, 0, 4*sizeof(Float_t));
fYref[0] = 0.; fYref[1] = 0.;
fZref[0] = 0.; fZref[1] = 0.;
fYfit[0] = 0.; fYfit[1] = 0.;
Double_t fSnp = trk->GetSnp();
Double_t fTgl = trk->GetTgl();
fPt = trk->Pt();
- Double_t norm =1./TMath::Sqrt(1. - fSnp*fSnp);
+ Double_t norm =1./TMath::Sqrt((1.-fSnp)*(1.+fSnp));
fYref[1] = fSnp*norm;
fZref[1] = fTgl*norm;
SetCovRef(trk->GetCovariance());
// 3. cluster size
//
- Int_t nclusters[kNslices];
- memset(nclusters, 0, kNslices*sizeof(Int_t));
memset(fdEdx, 0, kNslices*sizeof(Float_t));
-
const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
- AliTRDcluster *c = NULL;
+ 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());
-
+
// Filter clusters for dE/dx calculation
-
+
// 1.consider calibration effects for slice determination
Int_t slice;
- if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
+ if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
slice = Int_t(dx * nslices / kDriftLength);
} else slice = c->GetX() < fX0 ? nslices-1 : 0;
// 2. take sharing into account
Float_t w = /*c->IsShared() ? .5 :*/ 1.;
-
+
// 3. take into account large clusters TODO
//w *= c->GetNPads() > 3 ? .8 : 1.;
-
+
//CHECK !!!
fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
- nclusters[slice]++;
} // End of loop over clusters
-
- //if(fkReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
- if(nslices == AliTRDpidUtil::kLQslices){
- // calculate mean charge per slice (only LQ PID)
- for(int is=0; is<nslices; is++){
- if(nclusters[is]) fdEdx[is] /= nclusters[is];
- }
- }
}
//_____________________________________________________________________________
if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
}
+//____________________________________________________________
+Float_t AliTRDseedV1::GetAnodeWireOffset(Float_t zt)
+{
+// Find position inside the amplification cell for reading drift velocity map
+
+ Float_t d = fPad[3] - zt;
+ if(d<0.){
+ AliError(Form("Fail AnodeWireOffset calculation z0[%+7.2f] zt[%+7.2f] d[%+7.2f].", fPad[3], zt, d));
+ return 0.125;
+ }
+ d -= ((Int_t)(2 * d)) / 2.0;
+ if(d > 0.25) d = 0.5 - d;
+ return d;
+}
+
+
+//____________________________________________________________________
+Float_t AliTRDseedV1::GetCharge(Bool_t useOutliers)
+{
+// Computes total charge attached to tracklet. If "useOutliers" is set clusters
+// which are not in chamber are also used (default false)
+
+ AliTRDcluster *c(NULL); Float_t qt(0.);
+ for(int ic=0; ic<kNclusters; ic++){
+ if(!(c=fClusters[ic])) continue;
+ if(c->IsInChamber() && !useOutliers) continue;
+ qt += TMath::Abs(c->GetQ());
+ }
+ return qt;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::GetEstimatedCrossPoint(Float_t &x, Float_t &z) const
+{
+// Algorithm to estimate cross point in the x-z plane for pad row cross tracklets.
+// Returns true in case of success.
+ if(!IsRowCross()) return kFALSE;
+
+ x=0.; z=0.;
+ AliTRDcluster *c(NULL);
+ // Find radial range for first row
+ Float_t x1[] = {0., 1.e3};
+ for(int ic=0; ic<kNtb; ic++){
+ if(!(c=fClusters[ic])) continue;
+ if(!c->IsInChamber()) continue;
+ if(c->GetX() <= x1[1]) x1[1] = c->GetX();
+ if(c->GetX() >= x1[0]) x1[0] = c->GetX();
+ z=c->GetZ();
+ }
+ if((x1[0] - x1[1])<1.e-5) return kFALSE;
+
+ // Find radial range for second row
+ Bool_t kZ(kFALSE);
+ Float_t x2[] = {0., 1.e3};
+ for(int ic=kNtb; ic<kNclusters; ic++){
+ if(!(c=fClusters[ic])) continue;
+ if(!c->IsInChamber()) continue;
+ if(c->GetX() <= x2[1]) x2[1] = c->GetX();
+ if(c->GetX() >= x2[0]) x2[0] = c->GetX();
+ if(!kZ){
+ z+=c->GetZ();
+ z*=0.5;
+ kZ=kTRUE;
+ }
+ }
+ if((x2[0] - x2[1])<1.e-5) return kFALSE;
+
+ // Find intersection of the 2 radial regions
+ x = 0.5*((x1[0]+x1[1] > x2[0]+x2[1]) ? (x1[1]+x2[0]) : (x1[0]+x2[1]));
+ return kTRUE;
+}
//____________________________________________________________________
Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
if(fClusters[ic] && fClusters[ic]->IsInChamber()){
hasClusterInChamber = kTRUE;
dq += TMath::Abs(fClusters[ic]->GetQ());
- }else if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
+ }
+ if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
hasClusterInChamber = kTRUE;
dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
}
}
dx *= TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
if(dl) (*dl) = dx;
- return dq/dx;
+ if(dx>1.e-9) return dq/dx;
+ else return 0.;
}
//____________________________________________________________
return p;
}
+//____________________________________________________________________
+Float_t AliTRDseedV1::GetOccupancyTB() const
+{
+// Returns procentage of TB occupied by clusters
+
+ Int_t n(0);
+ AliTRDcluster *c(NULL);
+ for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
+ if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
+ n++;
+ }
+
+ return Float_t(n)/AliTRDtrackerV1::GetNTimeBins();
+}
//____________________________________________________________________
Float_t* AliTRDseedV1::GetProbability(Bool_t force)
Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/ TMath::Sqrt((1.0 - GetSnp()*GetSnp()) / (1.0 + GetTgl()*GetTgl()));
//calculate dE/dx
- CookdEdx(fkReconstructor->GetNdEdxSlices());
- AliDebug(4, Form("PID p[%f] dEdx[%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f] l[%f]", GetMomentum(), fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7], length));
+ CookdEdx(AliTRDCalPID::kNSlicesNN);
+ AliDebug(4, Form("p=%6.4f[GeV/c] dEdx{%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f} l=%4.2f[cm]", GetMomentum(), fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7], length));
// Sets the a priori probabilities
+ Bool_t kPIDNN(fkReconstructor->GetPIDMethod()==AliTRDpidUtil::kNN);
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
- fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane());
+ fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices());
return kTRUE;
}
sy2 += sys[0];
sz2 += sys[1];
}
- // rotate covariance matrix
- Double_t t2 = GetTilt()*GetTilt();
- Double_t correction = 1./(1. + t2);
- cov[0] = (sy2+t2*sz2)*correction;
- cov[1] = GetTilt()*(sz2 - sy2)*correction;
- cov[2] = (t2*sy2+sz2)*correction;
-
- //printf("C(%6.1f %+6.3f %6.1f) [%s]\n", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?" RC ":"-");
+
+ // rotate covariance matrix if no RC
+ if(!IsRowCross()){
+ Double_t t2 = GetTilt()*GetTilt();
+ Double_t correction = 1./(1. + t2);
+ cov[0] = (sy2+t2*sz2)*correction;
+ cov[1] = GetTilt()*(sz2 - sy2)*correction;
+ cov[2] = (t2*sy2+sz2)*correction;
+ } else {
+ cov[0] = sy2; cov[1] = 0.; cov[2] = sy2;
+ }
+
+ AliDebug(4, Form("C(%6.1f %+6.3f %6.1f) RC[%c]", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?'y':'n'));
}
//____________________________________________________________
-Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
+Int_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
+ const Double_t kZero(1.e-20);
Double_t l[2], // eigenvalues
v[3]; // eigenvectors
// the secular equation and its solution :
// 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.;
+ if(TMath::Abs(det)<kZero) 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.;
-
+ l[0] = .5*(tr + dd*(c[0]>c[2]?-1.:1.));
+ l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.));
+ if(l[0]<kZero || l[1]<kZero) return 2;
// 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 den = (l[0]-c[0])*(l[0]-c[0])+c[1]*c[1];
+ if(den<kZero){ // almost diagonal
+ v[0] = TMath::Sign(0., c[1]);
+ v[1] = TMath::Sign(1., (l[0]-c[0]));
+ v[2] = TMath::Sign(0., c[1]*(l[0]-c[0])*(l[1]-c[2]));
+ } else {
+ Double_t tmp = 1./TMath::Sqrt(den);
+ v[0] = c[1]* tmp;
+ v[1] = (l[0]-c[0])*tmp;
+ if(TMath::Abs(l[1]-c[2])<kZero) v[2] = TMath::Sign(v[0]*(l[0]-c[0])/kZero, (l[1]-c[2]));
+ else v[2] = v[0]*(l[0]-c[0])/(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];
- return 1.;
+ return 0;
}
//____________________________________________________________
//____________________________________________________________________
UShort_t AliTRDseedV1::GetVolumeId() const
{
- Int_t ic=0;
- while(ic<kNclusters && !fClusters[ic]) ic++;
- return fClusters[ic] ? fClusters[ic]->GetVolumeId() : 0;
-}
+// Returns geometry volume id by delegation
-//____________________________________________________________________
-TLinearFitter* AliTRDseedV1::GetFitterY()
-{
- if(!fgFitterY) fgFitterY = new TLinearFitter(1, "pol1");
- fgFitterY->ClearPoints();
- return fgFitterY;
+ for(Int_t ic(0);ic<kNclusters; ic++){
+ if(fClusters[ic]) return fClusters[ic]->GetVolumeId();
+ }
+ return 0;
}
-//____________________________________________________________________
-TLinearFitter* AliTRDseedV1::GetFitterZ()
-{
- if(!fgFitterZ) fgFitterZ = new TLinearFitter(1, "pol1");
- fgFitterZ->ClearPoints();
- return fgFitterZ;
-}
//____________________________________________________________________
void AliTRDseedV1::Calibrate()
fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
fDiffT, fVD);
+ AliDebug(4, Form("Calibration params for Det[%3d] Col[%3d] Row[%2d]\n t0[%f] vd[%f] s2PRF[%f] ExB[%f] Dl[%f] Dt[%f]", fDet, col, row, fT0, fVD, fS2PRF, fExB, fDiffL, fDiffT));
+
+
SetBit(kCalib, kTRUE);
}
}
//____________________________________________________________
-void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
+void AliTRDseedV1::SetPadPlane(AliTRDpadPlane * const p)
{
// Shortcut method to initialize pad geometry.
- if(!p) return;
- SetTilt(TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle()));
- SetPadLength(p->GetLengthIPad());
- SetPadWidth(p->GetWidthIPad());
+ fPad[0] = p->GetLengthIPad();
+ fPad[1] = p->GetWidthIPad();
+ fPad[2] = TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle());
+ fPad[3] = p->GetRow0() + p->GetAnodeWireOffset();
}
// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
// Debug : level >3
- if(!fkReconstructor->GetRecoParam() ){
+ const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
+
+ if(!recoParam){
AliError("Tracklets 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(fkReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
+ Int_t kClmin = Int_t(recoParam->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
- Double_t sysCov[5]; fkReconstructor->GetRecoParam()->GetSysCovMatrix(sysCov);
+ Double_t sysCov[5]; recoParam->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., //fkReconstructor->GetRecoParam() ->GetRoad1y();
- kroadz = GetPadLength() * fkReconstructor->GetRecoParam()->GetRoadzMultiplicator() + 1.;
+ Double_t kroady = 1., //recoParam->GetRoad1y();
+ kroadz = GetPadLength() * recoParam->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);
}
lr = ir; if(nr>=3) break;
}
- if(fkReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
+ if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
UChar_t stat(0);
if(IsKink()) SETBIT(stat, 1);
//
fkReconstructor = rec;
AliTRDgeometry g;
- AliTRDpadPlane *pp = g.GetPadPlane(fDet);
- fPad[0] = pp->GetLengthIPad();
- fPad[1] = pp->GetWidthIPad();
- fPad[3] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
+ SetPadPlane(g.GetPadPlane(fDet));
+
//fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
//fTgl = fZref[1];
Int_t n = 0, nshare = 0, nused = 0;
//____________________________________________________________________
-Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
+Bool_t AliTRDseedV1::Fit(UChar_t opt)
{
//
// Linear fit of the clusters attached to the tracklet
//
// Parameters :
-// - tilt : switch for tilt pad correction of cluster y position based on
-// the z, dzdx info from outside [default false].
-// - zcorr : switch for using z information to correct for anisochronity
-// and a finner error parameterization estimation [default false]
+// - opt : switch for tilt pad correction of cluster y position. Options are
+// 0 no correction [default]
+// 1 full tilt correction [dz/dx and z0]
+// 2 pseudo tilt correction [dz/dx from pad-chamber geometry]
+//
// Output :
// True if successful
//
// Author
// A.Bercuci <A.Bercuci@gsi.de>
+ if(!fkReconstructor){
+ AliError("The tracklet needs the reconstruction setup. Please initialize by SetReconstructor().");
+ return kFALSE;
+ }
if(!IsCalibrated()) Calibrate();
+ if(opt>2){
+ AliWarning(Form("Option [%d] outside range [0, 2]. Using default",opt));
+ opt=0;
+ }
const Int_t kClmin = 8;
-
+ const Float_t kScalePulls = 10.; // factor to scale y pulls - NOT UNDERSTOOD
// get track direction
Double_t y0 = fYref[0];
Double_t dydx = fYref[1];
Double_t z0 = fZref[0];
Double_t dzdx = fZref[1];
- Double_t yt, zt;
AliTRDtrackerV1::AliTRDLeastSquare fitterY;
AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
// book cluster information
Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
- Int_t n = 0;
- AliTRDcluster *c=NULL, **jc = &fClusters[0];
- for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
- xc[ic] = -1.;
- yc[ic] = 999.;
- zc[ic] = 999.;
- sy[ic] = 0.;
+ Bool_t tilt(opt==1) // full tilt correction
+ ,pseudo(opt==2) // pseudo tilt correction
+ ,rc(IsRowCross()) // row cross candidate
+ ,kDZDX(IsPrimary());// switch dzdx calculation for barrel primary tracks
+ Int_t n(0); // clusters used in fit
+ AliTRDcluster *c(NULL), *cc(NULL), **jc = &fClusters[0];
+ const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
+
+ const Char_t *tcName[]={"NONE", "FULL", "HALF"};
+ AliDebug(2, Form("Options : TC[%s] dzdx[%c]", tcName[opt], kDZDX?'Y':'N'));
+
+ for (Int_t ic=0; ic<kNclusters; ic++, ++jc) {
+ xc[ic] = -1.; yc[ic] = 999.; zc[ic] = 999.; sy[ic] = 0.;
if(!(c = (*jc))) continue;
if(!c->IsInChamber()) continue;
+ // compute pseudo tilt correction
+ if(kDZDX){
+ fZfit[0] = c->GetZ();
+ if(rc){
+ for(Int_t kc=AliTRDseedV1::kNtb; kc<AliTRDseedV1::kNclusters; kc++){
+ if(!(cc=fClusters[kc])) continue;
+ if(!cc->IsInChamber()) continue;
+ fZfit[0] += cc->GetZ(); fZfit[0] *= 0.5;
+ break;
+ }
+ }
+ fZfit[1] = fZfit[0]/fX0;
+ if(rc){
+ fZfit[0] += fZfit[1]*0.5*AliTRDgeometry::CdrHght();
+ fZfit[1] = fZfit[0]/fX0;
+ }
+ kDZDX=kFALSE;
+ }
Float_t w = 1.;
if(c->GetNPads()>4) w = .5;
qc[n] = TMath::Abs(c->GetQ());
// pad row of leading
- // Radial cluster position
- //Int_t jc = TMath::Max(fN-3, 0);
- //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
xc[n] = fX0 - c->GetX();
- // extrapolated track to cluster position
- yt = y0 - xc[n]*dydx;
- zt = z0 - xc[n]*dzdx;
-
// Recalculate cluster error based on tracking information
- c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
+ c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], -1./*zcorr?zt:-1.*/, dydx);
+ c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
sy[n] = TMath::Sqrt(c->GetSigmaY2());
- yc[n] = fkReconstructor->GetRecoParam()->UseGAUS() ?
+ yc[n] = recoParam->UseGAUS() ?
c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
zc[n] = c->GetZ();
- //optional tilt correction
- if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt));
- fitterY.AddPoint(&xc[n], yc[n], TMath::Sqrt(sy[n]));
- if(IsRowCross()) fitterZ.AddPoint(&xc[n], qc[n], 1.);
+ //optional r-phi correction
+ //printf(" n[%2d] yc[%7.5f] ", n, yc[n]);
+ Float_t correction(0.);
+ if(tilt) correction = fPad[2]*(xc[n]*dzdx + zc[n] - z0);
+ else if(pseudo) correction = fPad[2]*(xc[n]*fZfit[1] + zc[n]-fZfit[0]);
+ yc[n]-=correction;
+ //printf("corr(%s%s)[%7.5f] yc1[%7.5f]\n", (tilt?"TC":""), (zcorr?"PC":""), correction, yc[n]);
+
+ AliDebug(5, Form(" tb[%2d] dx[%6.3f] y[%6.2f+-%6.3f]", c->GetLocalTimeBin(), xc[n], yc[n], sy[n]));
+ fitterY.AddPoint(&xc[n], yc[n], sy[n]);
+ if(rc) fitterZ.AddPoint(&xc[n], qc[n]*(ic<kNtb?1.:-1.), 1.);
n++;
}
// to few clusters
- if (n < kClmin) return kFALSE;
-
+ if (n < kClmin){
+ AliDebug(1, Form("Not enough clusters to fit. Clusters: Attached[%d] Fit[%d].", GetN(), n));
+ SetErrorMsg(kFitCl);
+ return kFALSE;
+ }
// fit XY
- fitterY.Eval();
+ if(!fitterY.Eval()){
+ AliDebug(1, "Fit Y failed.");
+ SetErrorMsg(kFitFailedY);
+ return kFALSE;
+ }
fYfit[0] = fitterY.GetFunctionParameter(0);
fYfit[1] = -fitterY.GetFunctionParameter(1);
// store covariance
Double_t p[3];
fitterY.GetCovarianceMatrix(p);
- fCov[0] = p[0]; // variance of y0
- fCov[1] = p[2]; // covariance of y0, dydx
- fCov[2] = p[1]; // variance of dydx
+ fCov[0] = kScalePulls*p[1]; // variance of y0
+ fCov[1] = kScalePulls*p[2]; // covariance of y0, dydx
+ fCov[2] = kScalePulls*p[0]; // variance of dydx
// the ref radial position is set at the minimum of
// the y variance of the tracklet
fX = -fCov[1]/fCov[2];
+ fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
- // collect second row clusters
- Int_t m(0);
- if(IsRowCross()){
-/* // THE LEADING CLUSTER METHOD
+ Float_t xs=fX+.5*AliTRDgeometry::CamHght();
+ if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
+ AliDebug(1, Form("Ref radial position ouside chamber x[%5.2f].", fX));
+ SetErrorMsg(kFitFailedY);
+ return kFALSE;
+ }
+
+/* // THE LEADING CLUSTER METHOD for z fit
Float_t xMin = fX0;
Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
fS2Z = fdX*fZref[1];
fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
- // THE FIT X-Q PLANE METHOD
- Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
- for(; ic>kNtb; ic--, --jc){
- if(!(c = (*jc))) continue;
- if(!c->IsInChamber()) continue;
- qc[n] = TMath::Abs(c->GetQ());
- xc[n] = fX0 - c->GetX();
- zc[n] = c->GetZ();
- fitterZ.AddPoint(&xc[n], -qc[n], 1.);
- n--;m++;
+ // fit QZ
+ if(opt!=1 && IsRowCross()){
+ if(!fitterZ.Eval()) SetErrorMsg(kFitFailedZ);
+ if(!HasError(kFitFailedZ) && TMath::Abs(fitterZ.GetFunctionParameter(1))>1.e-10){
+ // TODO - one has to recalculate xy fit based on
+ // better knowledge of z position
+// Double_t x = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
+// Double_t z0 = .5*(zc[0]+zc[n-1]);
+// fZfit[0] = z0 + fZfit[1]*x;
+// fZfit[1] = fZfit[0]/fX0;
+// redo fit on xy plane
}
- }
- // fit XZ
- if(m && IsRowCross()){
- fitterZ.Eval();
- if(fitterZ.GetFunctionParameter(1)!=0.){
- fX = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
- fX=(fX<0.)?0.:fX;
- Float_t dl = .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
- fX=(fX> dl)?dl:fX;
- fX-=.055; // TODO to be understood
- }
-
- fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
// temporary external error parameterization
fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
// TODO correct formula
//fS2Z = sigma_x*TMath::Abs(fZref[1]);
} else {
- if(IsRowCross() && !m){
- AliDebug(1, "Tracklet crossed row but no clusters found in neighbor row.");
- }
- fZfit[0] = zc[0]; fZfit[1] = 0.;
+ //fZfit[0] = zc[0] + dzdx*0.5*AliTRDgeometry::CdrHght();
fS2Z = GetPadLength()*GetPadLength()/12.;
}
- fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
return kTRUE;
}
AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
+ AliInfo(Form("CALIB PARAMS : T0[%5.2f] Vd[%5.2f] s2PRF[%5.2f] ExB[%5.2f] Dl[%5.2f] Dt[%5.2f]", fT0, fVD, fS2PRF, fExB, fDiffL, fDiffT));
Double_t cov[3], x=GetX();
GetCovAt(x, cov);
AliInfo(Form("Fit | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | ----- |", x, GetY(), TMath::Sqrt(cov[0]), GetZ(), TMath::Sqrt(cov[2]), fYfit[1]));
AliInfo(Form("Ref | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | %5.2f |", x, fYref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[0]), fZref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[2]), fYref[1], fZref[1]))
AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
+ if(IsStandAlone()) AliInfo(Form("C Rieman / Vertex [1/cm] = %f / %f", fC[0], fC[1]));
AliInfo(Form("dEdx [a.u.] = %f / %f / %f / %f / %f/ %f / %f / %f", fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7]));
AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
}
- if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
- if ( GetTilt() != inTracklet->GetTilt() ) return kFALSE;
- if ( GetPadLength() != inTracklet->GetPadLength() ) return kFALSE;
+ if ( TMath::Abs(fS2Y - inTracklet->fS2Y)>1.e-10 ) return kFALSE;
+ if ( TMath::Abs(GetTilt() - inTracklet->GetTilt())>1.e-10 ) return kFALSE;
+ if ( TMath::Abs(GetPadLength() - inTracklet->GetPadLength())>1.e-10 ) return kFALSE;
for (Int_t i = 0; i < kNclusters; i++){
// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
//if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
//if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
- if ( fC != inTracklet->fC ) return kFALSE;
+ if ( TMath::Abs(fC[0] - inTracklet->fC[0])>1.e-10 ) return kFALSE;
//if ( fCC != inTracklet->GetCC() ) return kFALSE;
- if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
+ if ( TMath::Abs(fChi2 - inTracklet->fChi2)>1.e-10 ) return kFALSE;
// if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
if ( fDet != inTracklet->fDet ) return kFALSE;
- if ( fPt != inTracklet->fPt ) return kFALSE;
- if ( fdX != inTracklet->fdX ) return kFALSE;
+ if ( TMath::Abs(fPt - inTracklet->fPt)>1.e-10 ) return kFALSE;
+ if ( TMath::Abs(fdX - inTracklet->fdX)>1.e-10 ) return kFALSE;
for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
AliTRDcluster *curCluster = fClusters[iCluster];