X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDseedV1.cxx;h=6819af98502983eaea2cdccd527e40e94d41b74f;hb=8fc736d78f1e6acbc7afb601e3ea60663d281a03;hp=d276e1b936fc19427ece5b313ad0369743d81751;hpb=5f1ae1e77d9e02701773b26c5ce9e080dccc0bc9;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx index d276e1b936f..6819af98502 100644 --- a/TRD/AliTRDseedV1.cxx +++ b/TRD/AliTRDseedV1.cxx @@ -36,8 +36,6 @@ //////////////////////////////////////////////////////////////////////////// #include "TMath.h" -#include "TLinearFitter.h" -#include "TClonesArray.h" // tmp #include #include "AliLog.h" @@ -62,9 +60,6 @@ ClassImp(AliTRDseedV1) -TLinearFitter *AliTRDseedV1::fgFitterY = NULL; -TLinearFitter *AliTRDseedV1::fgFitterZ = NULL; - //____________________________________________________________________ AliTRDseedV1::AliTRDseedV1(Int_t det) :AliTRDtrackletBase() @@ -77,6 +72,7 @@ AliTRDseedV1::AliTRDseedV1(Int_t det) ,fDiffL(0.) ,fDiffT(0.) ,fClusterIdx(0) + ,fErrorMsg(0) ,fN(0) ,fDet(det) ,fPt(0.) @@ -87,7 +83,6 @@ AliTRDseedV1::AliTRDseedV1(Int_t det) ,fZ(0.) ,fS2Y(0.) ,fS2Z(0.) - ,fC(0.) ,fChi2(0.) { // @@ -95,7 +90,7 @@ AliTRDseedV1::AliTRDseedV1(Int_t det) // 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.; @@ -105,6 +100,8 @@ AliTRDseedV1::AliTRDseedV1(Int_t det) 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; @@ -123,6 +120,7 @@ AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref) ,fDiffL(0.) ,fDiffT(0.) ,fClusterIdx(0) + ,fErrorMsg(0) ,fN(0) ,fDet(-1) ,fPt(0.) @@ -133,7 +131,6 @@ AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref) ,fZ(0.) ,fS2Y(0.) ,fS2Z(0.) - ,fC(0.) ,fChi2(0.) { // @@ -200,6 +197,7 @@ void AliTRDseedV1::Copy(TObject &ref) const target.fDiffL = fDiffL; target.fDiffT = fDiffT; target.fClusterIdx = 0; + target.fErrorMsg = fErrorMsg; target.fN = fN; target.fDet = fDet; target.fPt = fPt; @@ -210,12 +208,11 @@ void AliTRDseedV1::Copy(TObject &ref) const 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]; @@ -224,6 +221,7 @@ void AliTRDseedV1::Copy(TObject &ref) const 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); @@ -255,24 +253,28 @@ Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track) //_____________________________________________________________________________ -void AliTRDseedV1::Reset() +void AliTRDseedV1::Reset(Option_t *opt) { - // - // Reset seed - // +// +// Reset seed. If option opt="c" is given only cluster arrays are cleared. +// + for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1; + memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*)); + fN=0; SetBit(kRowCross, kFALSE); + if(strcmp(opt, "c")==0) return; + fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.; fDiffL=0.;fDiffT=0.; fClusterIdx=0; - fN=0; + fErrorMsg = 0; fDet=-1; 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.; - for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1; - 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.; @@ -295,7 +297,7 @@ void AliTRDseedV1::Update(const AliTRDtrackV1 *trk) 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()); @@ -386,13 +388,10 @@ void AliTRDseedV1::CookdEdx(Int_t nslices) // 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; icGetX()); @@ -414,16 +413,7 @@ void AliTRDseedV1::CookdEdx(Int_t nslices) //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 1) && (out[3] > 1)) fLabels[1] = out[2]; } +//____________________________________________________________ +Float_t AliTRDseedV1::GetAnodeWireOffset(Float_t zt) +{ + 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::GetdQdl(Int_t ic, Float_t *dl) const @@ -506,7 +509,8 @@ Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const } 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.; } //____________________________________________________________ @@ -541,6 +545,20 @@ Float_t AliTRDseedV1::GetMomentum(Float_t *err) const 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; icGetRecoParam() ->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])); + Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/ TMath::Sqrt((1.0 - GetSnp()*GetSnp()) / (1.0 + GetTgl()*GetTgl())); //calculate dE/dx - CookdEdx(fkReconstructor->GetNdEdxSlices()); - + 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; ispecGetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane()); + fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices()); return kTRUE; } @@ -677,18 +695,23 @@ void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const 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. @@ -707,6 +730,7 @@ Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d) // Author A.Bercuci // 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 : @@ -715,26 +739,33 @@ Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d) // 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)c[2]?-1.:1.)); + l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.)); + if(l[0]GetVolumeId() : 0; -} - -//____________________________________________________________________ -TLinearFitter* AliTRDseedV1::GetFitterY() -{ - if(!fgFitterY) fgFitterY = new TLinearFitter(1, "pol1"); - fgFitterY->ClearPoints(); - return fgFitterY; + for(Int_t ic(0);icGetVolumeId(); + } + return 0; } -//____________________________________________________________________ -TLinearFitter* AliTRDseedV1::GetFitterZ() -{ - if(!fgFitterZ) fgFitterZ = new TLinearFitter(1, "pol1"); - fgFitterZ->ClearPoints(); - return fgFitterZ; -} //____________________________________________________________________ void AliTRDseedV1::Calibrate() @@ -831,6 +848,9 @@ 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); } @@ -851,10 +871,10 @@ void AliTRDseedV1::SetOwner() void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *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(); } @@ -888,32 +908,34 @@ Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t // Author : Alexandru Bercuci // Debug : level >3 - Bool_t kPRINT = kFALSE; - if(!fkReconstructor->GetRecoParam() ){ - AliError("Seed can not be used without a valid RecoParam."); + 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); - if(fkReconstructor->IsHLT())cp.SetRPhiMethod(AliTRDcluster::kCOG); - Calibrate(); + if(fkReconstructor->IsHLT()) cp.SetRPhiMethod(AliTRDcluster::kCOG); + if(!IsCalibrated()) Calibrate(); - if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady); + AliDebug(4, ""); + AliDebug(4, Form("syKalman[%f] rY[%f] rZ[%f]", syRef, kroady, kroadz)); // working variables const Int_t kNrows = 16; @@ -944,7 +966,9 @@ Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t // get estimated road kroady = 3.*TMath::Sqrt(12.*(s2yTrk + s2yCl)); - if(kPRINT) printf(" %2d dx[%f] yt[%f] zt[%f] sT[um]=%6.2f sy[um]=%6.2f syTilt[um]=%6.2f yRoad[mm]=%f\n", it, dx, yt, zt, 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()), 1.e4*TMath::Sqrt(s2yCl), 1.e1*kroady); + AliDebug(5, Form(" %2d x[%f] yt[%f] zt[%f]", it, dx, yt, zt)); + + AliDebug(5, Form(" syTrk[um]=%6.2f syCl[um]=%6.2f syClTlt[um]=%6.2f Ry[mm]=%f", 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()), 1.e4*TMath::Sqrt(s2yCl), 1.e1*kroady)); // select clusters cond[0] = yt; cond[2] = kroady; @@ -961,7 +985,7 @@ Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t continue; }*/ Int_t r = c->GetPadRow(); - if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r); + AliDebug(5, Form(" -> dy[%f] yc[%f] r[%d]", TMath::Abs(dy), c->GetY(), r)); clst[r][ncl[r]] = c; blst[r][ncl[r]] = kTRUE; idxs[r][ncl[r]] = idx[ic]; @@ -969,116 +993,177 @@ Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t ncl[r]++; ncls++; if(ncl[r] >= kNcls) { - AliWarning(Form("Cluster candidates reached buffer limit %d. Some may be lost.", kNcls)); + AliWarning(Form("Cluster candidates row[%d] reached buffer limit[%d]. Some may be lost.", r, kNcls)); kBUFFER = kTRUE; break; } } if(kBUFFER) break; } - if(kPRINT) printf("Found %d clusters\n", ncls); - if(ncls0 && lr-ir != 1){ - if(kPRINT) printf("W - gap in rows attached !!\n"); + if(lr>0 && ir-lr != 1){ + AliDebug(2, "Rows attached not continuous. Turn on selection."); + kRowSelection=kTRUE; } - if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]); - // Evaluate truncated mean on the y direction - if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8)); - else { - mean = 0.; syDis = 0.; - continue; - } - if(fkReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){ - TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::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"; - } + AliDebug(5, Form(" r[%d] n[%d]", ir, ncl[ir])); + // Evaluate truncated mean on the y direction + if(ncl[ir] < 4) continue; + AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean[nr], syDis[nr], Int_t(ncl[ir]*.8)); // 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 + AliDebug(4, Form(" m_%d[%+5.3f (%5.3fs)] s[%f]", nr, mean[nr], TMath::Abs(mean[nr]/syDis[nr]), syDis[nr])); + // remove outliers based on a 3 sigmaDistr level Bool_t kFOUND = kFALSE; for(Int_t ic = ncl[ir]; ic--;){ - if(yres[ir][ic] - mean > 3. * syDis){ + if(yres[ir][ic] - mean[nr] > 3. * syDis[nr]){ blst[ir][ic] = kFALSE; continue; } - nrow[nr]++; kFOUND = kTRUE; + nrow[nr]++; rowId[nr]=ir; kFOUND = kTRUE; + } + if(kFOUND){ + vdy[nr].Use(nrow[nr], yres[ir]); + nr++; } - // exit loop - if(kFOUND) nr++; lr = ir; if(nr>=3) break; } - if(kPRINT) printf("lr[%d] nr[%d] nrow[0]=%d nrow[1]=%d nrow[2]=%d\n", lr, nr, nrow[0], nrow[1], nrow[2]); - - // classify cluster rows - Int_t row = -1; - switch(nr){ - case 1: - row = lr; - break; - case 2: - SetBit(kRowCross, kTRUE); // mark pad row crossing - if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;} - else{ - row = lr; lr = 1; - nrow[2] = nrow[1]; - nrow[1] = nrow[0]; - nrow[0] = nrow[2]; + if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){ + TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker); + UChar_t stat(0); + if(IsKink()) SETBIT(stat, 1); + if(IsStandAlone()) SETBIT(stat, 2); + cstreamer << "AttachClusters" + << "stat=" << stat + << "det=" << fDet + << "pt=" << fPt + << "s2y=" << s2yTrk + << "r0=" << rowId[0] + << "dy0=" << &vdy[0] + << "m0=" << mean[0] + << "s0=" << syDis[0] + << "r1=" << rowId[1] + << "dy1=" << &vdy[1] + << "m1=" << mean[1] + << "s1=" << syDis[1] + << "r2=" << rowId[2] + << "dy2=" << &vdy[2] + << "m2=" << mean[2] + << "s2=" << syDis[2] + << "\n"; + } + + + // analyze gap in rows attached + if(kRowSelection){ + SetErrorMsg(kAttachRowGap); + Int_t rowRemove(-1); + if(nr==2){ // select based on minimum distance to track projection + if(TMath::Abs(mean[0])nrow[0]) AliDebug(2, Form("Conflicting mean[%f < %f] but ncl[%d < %d].", TMath::Abs(mean[0]), TMath::Abs(mean[1]), nrow[0], nrow[1])); + }else{ + if(nrow[1] %f] but ncl[%d > %d].", TMath::Abs(mean[0]), TMath::Abs(mean[1]), nrow[0], nrow[1])); + Swap(nrow[0],nrow[1]); Swap(rowId[0],rowId[1]); + Swap(mean[0],mean[1]); Swap(syDis[0],syDis[1]); + } + rowRemove=1; nr=1; + } else if(nr==3){ // select based on 2 consecutive rows + if(rowId[1]==rowId[0]+1 && rowId[1]!=rowId[2]-1){ + nr=2;rowRemove=2; + } else if(rowId[1]!=rowId[0]+1 && rowId[1]==rowId[2]-1){ + Swap(nrow[0],nrow[2]); Swap(rowId[0],rowId[2]); + Swap(mean[0],mean[2]); Swap(syDis[0],syDis[2]); + nr=2; rowRemove=2; + } } - break; - case 3: + if(rowRemove>0){nrow[rowRemove]=0; rowId[rowRemove]=-1;} + } + AliDebug(4, Form(" Ncl[%d[%d] + %d[%d] + %d[%d]]", nrow[0], rowId[0], nrow[1], rowId[1], nrow[2], rowId[2])); + + if(nr==3){ SetBit(kRowCross, kTRUE); // mark pad row crossing - break; + SetErrorMsg(kAttachRow); + const Float_t am[]={TMath::Abs(mean[0]), TMath::Abs(mean[1]), TMath::Abs(mean[2])}; + AliDebug(4, Form("complex row configuration\n" + " r[%d] n[%d] m[%6.3f] s[%6.3f]\n" + " r[%d] n[%d] m[%6.3f] s[%6.3f]\n" + " r[%d] n[%d] m[%6.3f] s[%6.3f]\n" + , rowId[0], nrow[0], am[0], syDis[0] + , rowId[1], nrow[1], am[1], syDis[1] + , rowId[2], nrow[2], am[2], syDis[2])); + Int_t id[]={0,1,2}; TMath::Sort(3, am, id, kFALSE); + // backup + Int_t rnn[3]; memcpy(rnn, nrow, 3*sizeof(Int_t)); + Int_t rid[3]; memcpy(rid, rowId, 3*sizeof(Int_t)); + Double_t rm[3]; memcpy(rm, mean, 3*sizeof(Double_t)); + Double_t rs[3]; memcpy(rs, syDis, 3*sizeof(Double_t)); + nrow[0]=rnn[id[0]]; rowId[0]=rid[id[0]]; mean[0]=rm[id[0]]; syDis[0]=rs[id[0]]; + nrow[1]=rnn[id[1]]; rowId[1]=rid[id[1]]; mean[1]=rm[id[1]]; syDis[1]=rs[id[1]]; + nrow[2]=0; rowId[2]=-1; mean[2] = 1.e3; syDis[2] = 1.e3; + AliDebug(4, Form("solved configuration\n" + " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n" + " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n" + " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n" + , rowId[0], nrow[0], mean[0], syDis[0] + , rowId[1], nrow[1], mean[1], syDis[1] + , rowId[2], nrow[2], mean[2], syDis[2])); + nr=2; + } else if(nr==2) { + SetBit(kRowCross, kTRUE); // mark pad row crossing + if(nrow[1] > nrow[0]){ // swap row order + Swap(nrow[0],nrow[1]); Swap(rowId[0],rowId[1]); + Swap(mean[0],mean[1]); Swap(syDis[0],syDis[1]); + } } - if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]); - if(row<0) return kFALSE; // Select and store clusters // We should consider here : // 1. How far is the chamber boundary // 2. How big is the mean - Int_t n = 0; + Int_t n(0); Float_t dyc[kNclusters]; memset(dyc,0,kNclusters*sizeof(Float_t)); for (Int_t ir = 0; ir < nr; ir++) { - Int_t jr = row + ir*lr; - if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr); + Int_t jr(rowId[ir]); + AliDebug(4, Form(" Attaching Ncl[%d]=%d ...", jr, ncl[jr])); for (Int_t ic = 0; ic < ncl[jr]; ic++) { if(!blst[jr][ic])continue; c = clst[jr][ic]; - Int_t it = c->GetPadTime(); + Int_t it(c->GetPadTime()); + Int_t idx(it+kNtb*ir); + if(fClusters[idx]){ + AliDebug(4, Form("Many cluster candidates on row[%2d] tb[%2d].", jr, it)); + // TODO should save also the information on where the multiplicity happened and its size + SetErrorMsg(kAttachMultipleCl); + // TODO should also compare with mean and sigma for this row + if(yres[jr][ic] > dyc[idx]) continue; + } + // TODO proper indexing of clusters !! - fIndexes[it+kNtb*ir] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]); - fClusters[it+kNtb*ir] = c; - - //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]); - + fIndexes[idx] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]); + fClusters[idx] = c; + dyc[idx] = yres[jr][ic]; n++; } - } + } + SetN(n); // number of minimum numbers of clusters expected for the tracklet - if (n < kClmin){ - //AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", n, kClmin)); + if (GetN() < kClmin){ + AliDebug(1, Form("NOT ENOUGH CLUSTERS %d ATTACHED TO THE TRACKLET [min %d] FROM FOUND %d.", GetN(), kClmin, n)); + SetErrorMsg(kAttachClAttach); return kFALSE; } - SetN(n); // Load calibration parameters for this tracklet Calibrate(); @@ -1112,10 +1197,8 @@ void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec) // 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; @@ -1134,16 +1217,17 @@ void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec) //____________________________________________________________________ -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 // @@ -1228,16 +1312,23 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) // Author // A.Bercuci + 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; @@ -1245,15 +1336,39 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) // 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; icGetRecoParam(); //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; icIsInChamber()) continue; + // compute pseudo tilt correction + if(kDZDX){ + fZfit[0] = c->GetZ(); + if(rc){ + for(Int_t kc=AliTRDseedV1::kNtb; kcIsInChamber()) 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; @@ -1263,50 +1378,63 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) 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 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]; @@ -1322,37 +1450,26 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) 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--; - } - // fit XZ - 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 + // fit QZ + if(opt!=1 && IsRowCross()){ + if(!fitterZ.Eval()) SetErrorMsg(kFitFailedZ); + if(!HasError(kFitFailedZ) && fitterZ.GetFunctionParameter(1)!=0.){ + // 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 } - - 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 { - 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; } @@ -1619,13 +1736,17 @@ void AliTRDseedV1::Print(Option_t *o) const 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(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |"); 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(strcmp(o, "a")!=0) return;