#include "AliLog.h"
#include "AliMathBase.h"
+#include "AliCDBManager.h"
+#include "AliTracker.h"
+#include "AliTRDpadPlane.h"
#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliTRDtrackerV1.h"
#include "AliTRDReconstructor.h"
#include "AliTRDrecoParam.h"
-#include "AliTRDgeometry.h"
+
#include "Cal/AliTRDCalPID.h"
+#include "Cal/AliTRDCalROC.h"
+#include "Cal/AliTRDCalDet.h"
ClassImp(AliTRDseedV1)
//____________________________________________________________________
-AliTRDseedV1::AliTRDseedV1(Int_t plane)
+AliTRDseedV1::AliTRDseedV1(Int_t det)
:AliTRDseed()
,fReconstructor(0x0)
- ,fPlane(plane)
+ ,fClusterIter(0x0)
+ ,fClusterIdx(0)
+ ,fDet(det)
,fMom(0.)
,fSnp(0.)
,fTgl(0.)
,fdX(0.)
+ ,fXref(0.)
+ ,fExB(0.)
{
//
// Constructor
//
- //printf("AliTRDseedV1::AliTRDseedV1()\n");
-
for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
+ fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
+ // covariance matrix [diagonal]
+ // default sy = 200um and sz = 2.3 cm
+ fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
}
//____________________________________________________________________
AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
:AliTRDseed((AliTRDseed&)ref)
,fReconstructor(ref.fReconstructor)
- ,fPlane(ref.fPlane)
+ ,fClusterIter(0x0)
+ ,fClusterIdx(0)
+ ,fDet(ref.fDet)
,fMom(ref.fMom)
,fSnp(ref.fSnp)
,fTgl(ref.fTgl)
,fdX(ref.fdX)
+ ,fXref(ref.fXref)
+ ,fExB(ref.fExB)
{
//
// Copy Constructor performing a deep copy
SetBit(kOwner, kFALSE);
for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
+ memcpy(fRefCov, ref.fRefCov, 3*sizeof(Double_t));
+ memcpy(fCov, ref.fCov, 3*sizeof(Double_t));
}
if(this != &ref){
ref.Copy(*this);
}
+ SetBit(kOwner, kFALSE);
+
return *this;
}
//AliInfo("");
AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
- target.fPlane = fPlane;
+ target.fClusterIter = 0x0;
+ target.fClusterIdx = 0;
+ target.fDet = fDet;
target.fMom = fMom;
target.fSnp = fSnp;
target.fTgl = fTgl;
target.fdX = fdX;
+ target.fXref = fXref;
+ target.fExB = fExB;
target.fReconstructor = fReconstructor;
for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
+ memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
+ memcpy(target.fCov, fCov, 3*sizeof(Double_t));
AliTRDseed::Copy(target);
}
fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
fZref[0] = z;
fZref[1] = track->GetTgl();
+
+ const Double_t *cov = track->GetCovariance();
+ fRefCov[0] = cov[0]; // Var(y)
+ fRefCov[1] = cov[1]; // Cov(yz)
+ fRefCov[2] = cov[5]; // Var(z)
//printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
return kTRUE;
nclusters[slice]++;
} // End of loop over clusters
- if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
+ //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
+ if(nslices == AliTRDReconstructor::kLQslices){
// calculate mean charge per slice (only LQ PID)
for(int is=0; is<nslices; is++){
if(nclusters[is]) fdEdx[is] /= nclusters[is];
}
}
+//____________________________________________________________________
+void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
+{
+// Return corrected position of the cluster taking into
+// account variation of the drift velocity with drift length.
+
+
+ // drift velocity correction TODO to be moved to the clusterizer
+ const Float_t cx[] = {
+ -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
+ -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
+ 3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
+ 3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
+ -1.2050e-01
+ };
+
+ // PRF correction TODO to be replaced by the gaussian
+ // approximation with full error parametrization and // moved to the clusterizer
+ const Float_t cy[AliTRDgeometry::kNlayer][3] = {
+ { 4.014e-04, 8.605e-03, -6.880e+00},
+ {-3.061e-04, 9.663e-03, -6.789e+00},
+ { 1.124e-03, 1.105e-02, -6.825e+00},
+ {-1.527e-03, 1.231e-02, -6.777e+00},
+ { 2.150e-03, 1.387e-02, -6.783e+00},
+ {-1.296e-03, 1.486e-02, -6.825e+00}
+ };
+
+ Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
+ x = c->GetX() - cx[c->GetLocalTimeBin()];
+ y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
+ return;
+}
//____________________________________________________________________
Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
// Sets the a priori probabilities
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
- fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane);
+ fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
}
return &fProb[0];
}
//____________________________________________________________________
-void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
+void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
{
-// Computes covariance in the y-z plane at radial point x
+// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
+// and returns the results in the preallocated array cov[3] as :
+// cov[0] = Var(y)
+// cov[1] = Cov(yz)
+// cov[2] = Var(z)
+//
+// Details
+//
+// For the linear transformation
+// BEGIN_LATEX
+// Y = T_{x} X^{T}
+// END_LATEX
+// The error propagation has the general form
+// BEGIN_LATEX
+// C_{Y} = T_{x} C_{X} T_{x}^{T}
+// END_LATEX
+// We apply this formula 2 times. First to calculate the covariance of the tracklet
+// at point x we consider:
+// BEGIN_LATEX
+// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
+// END_LATEX
+// and secondly to take into account the tilt angle
+// BEGIN_LATEX
+// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
+// END_LATEX
+//
+// using simple trigonometrics one can write for this last case
+// BEGIN_LATEX
+// C_{Y}=#frac{1}{1+tg^{2}#alpha} #(){#splitline{(#sigma_{y}^{2}+tg^{2}#alpha#sigma_{z}^{2}) __ tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2})}{tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2}) __ (#sigma_{z}^{2}+tg^{2}#alpha#sigma_{y}^{2})}}
+// END_LATEX
+// which can be aproximated for small alphas (2 deg) with
+// BEGIN_LATEX
+// C_{Y}=#(){#splitline{#sigma_{y}^{2} __ (#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha}{((#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha __ #sigma_{z}^{2}}}
+// END_LATEX
+//
+// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
+// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
+// account for extra misalignment/miscalibration uncertainties.
+//
+// Author :
+// Alex Bercuci <A.Bercuci@gsi.de>
+// Date : Jan 8th 2009
+//
+ Double_t xr = fX0-x;
+ Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
+ Double_t sz2 = fPadLength*fPadLength/12.;
- Int_t ic = 0; while (!fClusters[ic]) ic++;
- AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
- Double_t exB = fCalib->GetOmegaTau(fCalib->GetVdriftAverage(fClusters[ic]->GetDetector()), -AliTracker::GetBz()*0.1);
+ // insert systematic uncertainties
+ Double_t sys[15];
+ fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
+ sy2 += sys[0];
+ sz2 += sys[1];
+
+ // rotate covariance matrix
+ Double_t t2 = fTilt*fTilt;
+ Double_t correction = 1./(1. + t2);
+ cov[0] = (sy2+t2*sz2)*correction;
+ cov[1] = fTilt*(sz2 - sy2)*correction;
+ cov[2] = (t2*sy2+sz2)*correction;
+}
- Double_t sy2 = fSigmaY2*fSigmaY2 + .2*(fYfit[1]-exB)*(fYfit[1]-exB);
- Double_t sz2 = fPadLength/12.;
+//____________________________________________________________________
+void AliTRDseedV1::SetExB()
+{
+// Retrive the tg(a_L) from OCDB. The following information are used
+// - detector index
+// - column and row position of first attached cluster.
+//
+// If no clusters are attached to the tracklet a random central chamber
+// position (c=70, r=7) will be used to retrieve the drift velocity.
+//
+// Author :
+// Alex Bercuci <A.Bercuci@gsi.de>
+// Date : Jan 8th 2009
+//
- //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
+ AliCDBManager *cdb = AliCDBManager::Instance();
+ if(cdb->GetRun() < 0){
+ AliError("OCDB manager not properly initialized");
+ return;
+ }
- cov[0] = sy2;
- cov[1] = fTilt*(sy2-sz2);
- cov[2] = sz2;
+ AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
+ AliTRDCalROC *fCalVdriftROC = fCalib->GetVdriftROC(fDet);
+ const AliTRDCalDet *fCalVdriftDet = fCalib->GetVdriftDet();
+
+ Int_t col = 70, row = 7;
+ AliTRDcluster **c = &fClusters[0];
+ if(fN){
+ Int_t ic = 0;
+ while (ic<AliTRDseed::knTimebins && !(*c)){ic++; c++;}
+ if(*c){
+ col = (*c)->GetPadCol();
+ row = (*c)->GetPadRow();
+ }
+ }
- // insert systematic uncertainties calibration and misalignment
- Double_t sys[15];
- fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
- cov[0] += (sys[0]*sys[0]);
- cov[2] += (sys[1]*sys[1]);
+ Double_t vd = fCalVdriftDet->GetValue(fDet) * fCalVdriftROC->GetValue(col, row);
+ fExB = fCalib->GetOmegaTau(vd, -0.1*AliTracker::GetBz());
}
-
//____________________________________________________________________
void AliTRDseedV1::SetOwner()
{
}
AliTRDchamberTimeBin *layer = 0x0;
- if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7 && c){
- TClonesArray clusters("AliTRDcluster", 24);
- clusters.SetOwner(kTRUE);
- AliTRDcluster *cc = 0x0;
- Int_t det=-1, ncl, ncls = 0;
- for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
- if(!(layer = chamber->GetTB(iTime))) continue;
- if(!(ncl = Int_t(*layer))) continue;
- for(int ic=0; ic<ncl; ic++){
- cc = (*layer)[ic];
- det = cc->GetDetector();
- new(clusters[ncls++]) AliTRDcluster(*cc);
- }
- }
- AliInfo(Form("N clusters[%d] = %d", fPlane, ncls));
-
- Int_t ref = c ? 1 : 0;
- TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
+ AliTRDtrackingChamber ch(*chamber);
+ ch.SetOwner();
+ TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
cstreamer << "AttachClustersIter"
- << "det=" << det
- << "ref=" << ref
- << "clusters.=" << &clusters
+ << "chamber.=" << &ch
<< "tracklet.=" << this
- << "cl.=" << c
<< "\n";
}
fZ[iTime] = cl->GetZ();
ncl++;
}
- if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
if(ncl>1){
// calculate length of the time bin (calibration aware)
for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
if(!fClusters[jTime]) continue;
fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
+ break;
}
- break;
}
}
// update x reference positions (calibration/alignment aware)
for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!fClusters[iTime]) continue;
- fX[iTime] = fClusters[iTime]->GetX() - fX0;
+ fX[iTime] = fX0 - fClusters[iTime]->GetX();
}
AliTRDseed::Update();
}
- if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
if(IsOK()){
tquality = GetQuality(kZcorr);
} // Loop: iter
if (!IsOK()) return kFALSE;
- CookLabels();
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
+
+ // refit tracklet with errors
+ //SetExB(); Fit(kFALSE, 2);
+
UpdateUsed();
return kTRUE;
}
return kTRUE;
}
+//____________________________________________________________
+void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
+{
+// Fill in all derived information. It has to be called after recovery from file or HLT.
+// The primitive data are
+// - list of clusters
+// - detector (as the detector will be removed from clusters)
+// - position of anode wire (fX0) - temporary
+// - track reference position and direction
+// - momentum of the track
+// - time bin length [cm]
+//
+// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
+//
+ fReconstructor = rec;
+ AliTRDgeometry g;
+ AliTRDpadPlane *pp = g.GetPadPlane(fDet);
+ fTilt = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
+ fPadLength = pp->GetLengthIPad();
+ fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
+ fTgl = fZref[1];
+ fN = 0; fN2 = 0; fMPads = 0.;
+ AliTRDcluster **cit = &fClusters[0];
+ for(Int_t ic = knTimebins; ic--; cit++){
+ if(!(*cit)) return;
+ fN++; fN2++;
+ fX[ic] = (*cit)->GetX() - fX0;
+ fY[ic] = (*cit)->GetY();
+ fZ[ic] = (*cit)->GetZ();
+ }
+ Update(); // Fit();
+ CookLabels();
+ GetProbability();
+}
+
+
//____________________________________________________________________
-Bool_t AliTRDseedV1::Fit()
+Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
{
//
// Linear fit of the tracklet
//
const Int_t kClmin = 8;
- const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
- AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
+
+ // cluster error parametrization parameters
+ // 1. sy total charge
+ const Float_t sq0inv = 0.019962; // [1/q0]
+ const Float_t sqb = 1.0281564; //[cm]
+ // 2. sy for the PRF
+ const Float_t scy[AliTRDgeometry::kNlayer][4] = {
+ {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
+ {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
+ {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
+ {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
+ {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
+ {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
+ };
+ // 3. sy parallel to the track
+ const Float_t sy0 = 2.649e-02; // [cm]
+ const Float_t sya = -8.864e-04; // [cm]
+ const Float_t syb = -2.435e-01; // [cm]
+
+ // 4. sx parallel to the track
+ const Float_t sxgc = 5.427e-02;
+ const Float_t sxgm = 7.783e-01;
+ const Float_t sxgs = 2.743e-01;
+ const Float_t sxe0 =-2.065e+00;
+ const Float_t sxe1 =-2.978e-02;
+
+ // 5. sx perpendicular to the track
+// const Float_t sxd0 = 1.881e-02;
+// const Float_t sxd1 =-4.101e-01;
+// const Float_t sxd2 = 1.572e+00;
+
+ // 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;
+
+ const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
+ AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
+ TLinearFitter fitterY(1, "pol1");
// convertion factor from square to gauss distribution for sigma
Double_t convert = 1./TMath::Sqrt(12.);
-
+
// book cluster information
- Double_t xc[knTimebins+1], yc[knTimebins], zc[knTimebins+1], sy[knTimebins], sz[knTimebins+1];
+ Double_t q, xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
Int_t zRow[knTimebins];
- AliTRDcluster *c = 0x0;
- Int_t nc = 0;
- for (Int_t ic=0; ic<kNtb; ic++) {
+
+ Int_t ily = AliTRDgeometry::GetLayer(fDet);
+ fN = 0; fXref = 0.; Double_t ssx = 0.;
+ AliTRDcluster *c=0x0, **jc = &fClusters[0];
+ for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
zRow[ic] = -1;
xc[ic] = -1.;
yc[ic] = 999.;
zc[ic] = 999.;
sy[ic] = 0.;
sz[ic] = 0.;
- if(!(c = fClusters[ic])) continue;
+ if(!(c = (*jc))) continue;
if(!c->IsInChamber()) continue;
+
Float_t w = 1.;
if(c->GetNPads()>4) w = .5;
if(c->GetNPads()>5) w = .2;
- zRow[nc] = c->GetPadRow();
- xc[nc] = fX0 - c->GetX();
- yc[nc] = c->GetY();
- zc[nc] = c->GetZ();
- sy[ic] = w; // all clusters have the same sigma
- sz[ic] = fPadLength*convert;
- fitterZ.AddPoint(&xc[ic], zc[ic], sz[ic]);
- nc++;
+
+ zRow[fN] = c->GetPadRow();
+ // correct cluster position for PRF and v drift
+ Double_t x, y; GetClusterXY(c, x, y);
+ xc[fN] = fX0 - x;
+ yc[fN] = y;
+ zc[fN] = c->GetZ();
+
+ // extrapolated y value for the track
+ yt = y0 - xc[fN]*dydx;
+ // extrapolated z value for the track
+ zt = z0 - xc[fN]*dzdx;
+ // tilt correction
+ if(tilt) yc[fN] -= fTilt*(zc[fN] - zt);
+
+ // ELABORATE CLUSTER ERROR
+ // TODO to be moved to AliTRDcluster
+ q = TMath::Abs(c->GetQ());
+ Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
+ // basic y error (|| to track).
+ sy[fN] = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
+ //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN, sy[fN]*1.e4);
+ // y error due to total charge
+ sy[fN] += sqb*(1./q - sq0inv);
+ //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
+ // y error due to PRF
+ sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
+ //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
+
+ sy[fN] *= sy[fN];
+
+ // ADD ERROR ON x
+ // error of drift length parallel to the track
+ Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
+ //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
+ // error of drift length perpendicular to the track
+ //sx += sxd0 + sxd1*d + sxd2*d*d;
+ sx *= sx; // square sx
+ // update xref
+ fXref += xc[fN]/sx; ssx+=1./sx;
+
+ // add error from ExB
+ if(errors>0) sy[fN] += fExB*fExB*sx;
+ //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
+
+ // global radial error due to misalignment/miscalibration
+ Double_t sx0 = 0.; sx0 *= sx0;
+ // add sx contribution to sy due to track angle
+ if(errors>1) sy[fN] += tgg*(sx+sx0);
+ // TODO we should add tilt pad correction here
+ //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
+ c->SetSigmaY2(sy[fN]);
+
+ sy[fN] = TMath::Sqrt(sy[fN]);
+ fitterY.AddPoint(&xc[fN], yc[fN]/*-yt*/, sy[fN]);
+
+ sz[fN] = fPadLength*convert;
+ fitterZ.AddPoint(&xc[fN], zc[fN], sz[fN]);
+ fN++;
}
// to few clusters
- if (nc < kClmin) return kFALSE;
-
+ if (fN < kClmin) return kFALSE;
- Int_t zN[2*35];
- Int_t nz = AliTRDtrackerV1::Freq(nc, zRow, zN, kFALSE);
+ // fit XY
+ fitterY.Eval();
+ fYfit[0] = fitterY.GetParameter(0);
+ fYfit[1] = -fitterY.GetParameter(1);
+ // store covariance
+ Double_t *p = fitterY.GetCovarianceMatrix();
+ fCov[0] = p[0]; // variance of y0
+ fCov[1] = p[1]; // covariance of y0, dydx
+ fCov[2] = p[3]; // variance of dydx
+ // store ref radial position.
+ fXref /= ssx; fXref = fX0 - fXref;
+
+ // check par row crossing
+ Int_t zN[2*AliTRDseed::knTimebins];
+ Int_t nz = AliTRDtrackerV1::Freq(fN, zRow, zN, kFALSE);
// more than one pad row crossing
if(nz>2) return kFALSE;
-
- // estimate reference parameter at average x
- Double_t y0 = fYref[0];
- Double_t dydx = fYref[1];
- Double_t dzdx = fZref[1];
- zc[nc] = fZref[0];
+
// determine z offset of the fit
+ Float_t zslope = 0.;
Int_t nchanges = 0, nCross = 0;
if(nz==2){ // tracklet is crossing pad row
// Find the break time allowing one chage on pad-rows
// with maximal number of accepted clusters
Int_t padRef = zRow[0];
- for (Int_t ic=1; ic<nc; ic++) {
+ for (Int_t ic=1; ic<fN; ic++) {
if(zRow[ic] == padRef) continue;
// debug
// evaluate parameters of the crossing point
Float_t sx = (xc[ic-1] - xc[ic])*convert;
- xc[nc] = .5 * (xc[ic-1] + xc[ic]);
- zc[nc] = .5 * (zc[ic-1] + zc[ic]);
- sz[nc] = TMath::Max(dzdx * sx, .01);
- dzdx = zc[ic-1] > zc[ic] ? 1. : -1.;
- padRef = zRow[ic];
- nCross = ic;
+ fCross[0] = .5 * (xc[ic-1] + xc[ic]);
+ fCross[2] = .5 * (zc[ic-1] + zc[ic]);
+ fCross[3] = TMath::Max(dzdx * sx, .01);
+ zslope = zc[ic-1] > zc[ic] ? 1. : -1.;
+ padRef = zRow[ic];
+ nCross = ic;
nchanges++;
}
}
// condition on nCross and reset nchanges TODO
if(nchanges==1){
- if(dzdx * fZref[1] < 0.){
+ if(dzdx * zslope < 0.){
AliInfo("tracklet direction does not correspond to the track direction. TODO.");
}
SetBit(kRowCross, kTRUE); // mark pad row crossing
- fCross[0] = xc[nc]; fCross[2] = zc[nc]; fCross[3] = sz[nc];
- fitterZ.AddPoint(&xc[nc], zc[nc], sz[nc]);
+ fitterZ.AddPoint(&fCross[0], fCross[2], fCross[3]);
fitterZ.Eval();
- dzdx = fZref[1]; // we don't trust Parameter[1] ??;
- zc[nc] = fitterZ.GetFunctionParameter(0);
+ //zc[nc] = fitterZ.GetFunctionParameter(0);
+ fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
+ fCross[0] = fX0 - fCross[0];
} else if(nchanges > 1){ // debug
- AliInfo("ERROR in n changes!!!");
+ AliError("N pad row crossing > 1.");
return kFALSE;
}
-
- // estimate deviation from reference direction
- dzdx *= fTilt;
- for (Int_t ic=0; ic<nc; ic++) {
- yc[ic] -= y0 + xc[ic]*(dydx + dzdx) + fTilt * (zc[ic] - zc[nc]);
- fitterY.AddPoint(&xc[ic], yc[ic], sy[ic]);
- }
- fitterY.Eval();
- fYfit[0] = y0+fitterY.GetFunctionParameter(0);
- fYfit[1] = dydx+fitterY.GetFunctionParameter(1);
- if(nchanges) fCross[1] = fYfit[0] + fCross[0] * fYfit[1];
-
-// printf("\nnz = %d\n", nz);
-// for(int ic=0; ic<35; ic++) printf("%d row[%d]\n", ic, zRow[ic]);
-//
-// for(int ic=0; ic<nz; ic++) printf("%d n[%d]\n", ic, zN[ic]);
+ UpdateUsed();
return kTRUE;
}
-//___________________________________________________________________
-void AliTRDseedV1::Draw(Option_t*)
-{
-}
//___________________________________________________________________
-void AliTRDseedV1::Print(Option_t*) const
+void AliTRDseedV1::Print(Option_t *o) const
{
//
// Printing the seedstatus
//
- printf("Seed status :\n");
- printf(" fTilt = %f\n", fTilt);
- printf(" fPadLength = %f\n", fPadLength);
- printf(" fX0 = %f\n", fX0);
- for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++) {
- const Char_t *isUsable = fUsable[ic]?"Yes":"No";
- printf(" %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
- , ic
- , fX[ic]
- , fY[ic]
- , fZ[ic]
- , fIndexes[ic]
- , ((void*) fClusters[ic])
- , isUsable);
- }
+ AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
+ AliInfo(Form("Nattach[%2d] Nfit[%2d] Nuse[%2d] pads[%f]", fN, fN2, fNUsed, fMPads));
+ AliInfo(Form("x[%7.2f] y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fX0, fYfit[0], fZfit[0], fYfit[1], fZfit[1]));
+ AliInfo(Form("Ref y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
- printf(" fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
- printf(" fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
- printf(" fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
- printf(" fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
- printf(" fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
- printf(" fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
- printf(" fSigmaY =%f\n", fSigmaY);
- printf(" fSigmaY2=%f\n", fSigmaY2);
- printf(" fMeanz =%f\n", fMeanz);
- printf(" fZProb =%f\n", fZProb);
- printf(" fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
- printf(" fN =%d\n", fN);
- printf(" fN2 =%d (>8 isOK)\n",fN2);
- printf(" fNUsed =%d\n", fNUsed);
- printf(" fFreq =%d\n", fFreq);
- printf(" fNChange=%d\n", fNChange);
- printf(" fMPads =%f\n", fMPads);
-
- printf(" fC =%f\n", fC);
- printf(" fCC =%f\n",fCC);
- printf(" fChi2 =%f\n", fChi2);
- printf(" fChi2Z =%f\n", fChi2Z);
+
+ if(strcmp(o, "a")!=0) return;
+
+ AliTRDcluster* const* jc = &fClusters[0];
+ for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
+ if(!(*jc)) continue;
+ (*jc)->Print(o);
+ }
}
+
+//___________________________________________________________________
+Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
+{
+ // Checks if current instance of the class has the same essential members
+ // as the given one
+
+ if(!o) return kFALSE;
+ const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
+ if(!inTracklet) return kFALSE;
+
+ for (Int_t i = 0; i < 2; i++){
+ if ( fYref[i] != inTracklet->GetYref(i) ) return kFALSE;
+ if ( fZref[i] != inTracklet->GetZref(i) ) return kFALSE;
+ }
+
+ if ( fSigmaY != inTracklet->GetSigmaY() ) return kFALSE;
+ if ( fSigmaY2 != inTracklet->GetSigmaY2() ) return kFALSE;
+ if ( fTilt != inTracklet->GetTilt() ) return kFALSE;
+ if ( fPadLength != inTracklet->GetPadLength() ) return kFALSE;
+
+ for (Int_t i = 0; i < knTimebins; i++){
+ if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
+ if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
+ if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
+ if ( fIndexes[i] != inTracklet->GetIndexes(i) ) return kFALSE;
+ if ( fUsable[i] != inTracklet->IsUsable(i) ) return kFALSE;
+ }
+
+ for (Int_t i=0; i < 2; i++){
+ if ( fYfit[i] != inTracklet->GetYfit(i) ) return kFALSE;
+ if ( fZfit[i] != inTracklet->GetZfit(i) ) return kFALSE;
+ if ( fYfitR[i] != inTracklet->GetYfitR(i) ) return kFALSE;
+ if ( fZfitR[i] != inTracklet->GetZfitR(i) ) return kFALSE;
+ if ( fLabels[i] != inTracklet->GetLabels(i) ) return kFALSE;
+ }
+
+ if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
+ if ( fZProb != inTracklet->GetZProb() ) return kFALSE;
+ if ( fN2 != inTracklet->GetN2() ) return kFALSE;
+ if ( fNUsed != inTracklet->GetNUsed() ) return kFALSE;
+ if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
+ if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
+ if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
+
+ if ( fC != inTracklet->GetC() ) return kFALSE;
+ if ( fCC != inTracklet->GetCC() ) return kFALSE;
+ if ( fChi2 != inTracklet->GetChi2() ) return kFALSE;
+ // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
+
+ if ( fDet != inTracklet->GetDetector() ) return kFALSE;
+ if ( fMom != inTracklet->GetMomentum() ) return kFALSE;
+ if ( fdX != inTracklet->GetdX() ) return kFALSE;
+
+ for (Int_t iCluster = 0; iCluster < knTimebins; iCluster++){
+ AliTRDcluster *curCluster = fClusters[iCluster];
+ AliTRDcluster *inCluster = inTracklet->GetClusters(iCluster);
+ if (curCluster && inCluster){
+ if (! curCluster->IsEqual(inCluster) ) {
+ curCluster->Print();
+ inCluster->Print();
+ return kFALSE;
+ }
+ } else {
+ // if one cluster exists, and corresponding
+ // in other tracklet doesn't - return kFALSE
+ if(curCluster || inCluster) return kFALSE;
+ }
+ }
+ return kTRUE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff