--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+#include "AliTPCpolyTrack.h"
+#include "TMath.h"
+
+ClassImp(AliTPCpolyTrack)
+
+
+AliTPCpolyTrack::AliTPCpolyTrack()
+{
+ Reset();
+}
+
+void AliTPCpolyTrack::Reset()
+{
+ //
+ // reset track
+ fSumX = fSumX2= fSumX3=fSumX4 = fSumY=fSumYX=fSumYX2=fSumZ=fSumZX=fSumZX2=fSumW =0;
+ fNPoints = 0;
+}
+
+void AliTPCpolyTrack::AddPoint(Double_t x, Double_t y, Double_t z,Double_t sy, Double_t sz)
+{
+ //
+ //
+ if (fNPoints==0){
+ fMaxX = x;
+ fMinX = x;
+ }else{
+ if (x>fMaxX) fMaxX=x;
+ if (x<fMinX) fMinX=x;
+ }
+
+ Double_t x2 = x*x;
+ Double_t w = 2./(sy+sz);
+ fSumW += w;
+ //
+ fSumX += x*w;
+ fSumX2 += x2*w;
+ fSumX3 += x2*x*w;
+ fSumX4 += x2*x2*w;
+ //
+ fSumY +=y*w;
+ fSumYX +=y*x*w;
+ fSumYX2 +=y*x2*w;
+ //
+ fSumZ +=z*w;
+ fSumZX +=z*x*w;
+ fSumZX2 +=z*x2*w;
+ //
+ fX[fNPoints] = x;
+ fY[fNPoints] = y;
+ fZ[fNPoints] = z;
+ fSY[fNPoints] = sy;
+ fSZ[fNPoints] = sz;
+
+ fNPoints++;
+
+}
+
+void AliTPCpolyTrack::UpdateParameters()
+{
+ //
+ //
+ //Update fit parameters
+ if (fNPoints>4){
+ Fit2(fSumY,fSumYX,fSumYX2,fSumX,fSumX2,fSumX3,fSumX4,fSumW,fA,fB,fC);
+ // Fit2(fSumZ,fSumZX,fSumZX2,fSumX,fSumX2,fSumX3,fSumX4,fNPoints,fD,fE,fF);
+ Fit1(fSumZ,fSumZX,fSumX,fSumX2,fSumW,fD,fE,fF);
+ }
+ else
+ {
+ Fit1(fSumY,fSumYX,fSumX,fSumX2,fSumW,fA,fB,fC);
+ Fit1(fSumZ,fSumZX,fSumX,fSumX2,fSumW,fD,fE,fF);
+ }
+}
+
+
+void AliTPCpolyTrack::GetFitPoint(Double_t x, Double_t &y, Double_t &z)
+{
+ y = fA+fB*x+fC*x*x;
+ z = fD+fE*x+fF*x*x;
+}
+
+
+void AliTPCpolyTrack::Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
+ Double_t fSumX, Double_t fSumX2, Double_t fSumX3,
+ Double_t fSumX4, Double_t fSumW,
+ Double_t &a, Double_t &b, Double_t &c)
+{
+ //fit of second order
+ Double_t det =
+ fSumW* (fSumX2*fSumX4-fSumX3*fSumX3) -
+ fSumX* (fSumX*fSumX4-fSumX3*fSumX2)+
+ fSumX2* (fSumX*fSumX3-fSumX2*fSumX2);
+
+ if (TMath::Abs(det)> 0.000000000000001) {
+ a =
+ (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
+ fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
+ fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
+ b=
+ (fSumW*(fSumYX*fSumX4-fSumX3*fSumYX2)-
+ fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
+ fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
+ c=
+ (fSumW*(fSumX2*fSumYX2-fSumYX*fSumX3)-
+ fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
+ fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
+ }
+}
+
+void AliTPCpolyTrack::Fit1(Double_t fSumY, Double_t fSumYX,
+ Double_t fSumX, Double_t fSumX2,
+ Double_t fSumW, Double_t &a, Double_t &b, Double_t &c)
+{
+ //
+ //
+ //
+ Double_t det = fSumW*fSumX2-fSumX*fSumX;
+ if (TMath::Abs(det)> 0.000000000000001) {
+ b = (fSumW*fSumYX-fSumX*fSumY)/det;
+ a = (fSumX2*fSumY-fSumX*fSumYX)/det;
+ c = 0;
+ }else{
+ a =fSumYX/fSumX;
+ b =0;
+ c =0;
+ }
+
+}
#include "AliTPCParam.h"
#include "AliTPCClustersRow.h"
#include "AliComplexCluster.h"
+#include "AliTPCpolyTrack.h"
#include "TStopwatch.h"
}
+void AliTPCseed::Reset()
+{
+ //
+ //PH SetN(0);
+ fNFoundable = 0;
+ ResetCovariance();
+ SetChi2(0);
+ for (Int_t i=0;i<200;i++) fClusterIndex[i]=-1;
+}
+
Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
{
Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
- if (TMath::Abs(cur*fX-eta) >= 0.99999) {
- Int_t n=GetNumberOfClusters();
- if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
+ if (TMath::Abs(cur*fX-eta) >= 0.9) {
+ // Int_t n=GetNumberOfClusters();
+ //if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
return 0;
}
// if (row < nr) return 1; // don't prolongate if not information until now -
//
Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
- if (!t.PropagateTo(x)) return 0;
+ // if (t.GetRadius()>x+10 ) return 0;
+
+ if (!t.PropagateTo(x)) {
+ t.fStopped = kTRUE;
+ return 0;
+ }
// update current
t.fCurrentSigmaY = GetSigmaY(&t);
t.fCurrentSigmaZ = GetSigmaZ(&t);
}
else
{
- t.fNFoundable++;
+ if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
+ else
+ return 0;
}
//calculate
Float_t maxdistance = roady*roady + roadz*roadz;
Int_t row = fSectors->GetRowNumber(xt)-1;
if (row < nr) return 1; // don't prolongate if not information until now -
Double_t x=fSectors->GetX(nr);
- if (!t.PropagateTo(x)) return 0;
+ // if (t.fStopped) return 0;
+ // if (t.GetRadius()>x+10 ) return 0;
+ if (!t.PropagateTo(x)){
+ t.fStopped =kTRUE;
+ return 0;
+ }
// update current
t.fCurrentSigmaY = GetSigmaY(&t);
t.fCurrentSigmaZ = GetSigmaZ(&t);
}
else
{
- t.fNFoundable++;
+ if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
+ else
+ return 0;
}
if (t.fCurrentCluster) {
// printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
+ if ( (rdistancey>1) || (rdistancez>1)) return 0;
if (rdistance>4) return 0;
if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
- if (t.GetSnp()>0.2)
+ if (t.GetSnp()<0.8)
FollowToNext(t,nr);
}
return 1;
}
}
- Float_t summin = TMath::Min(sum1,sum2);
- Float_t ratio = sum/Float_t(summin);
+ Float_t summin = TMath::Min(sum1+1,sum2+1);
+ Float_t ratio = (sum+1)/Float_t(summin);
return ratio;
}
if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
Int_t sum1,sum2;
Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
+ //if (sum1==0) {
+ // pt->Desactivate(removalindex); // arr->RemoveAt(i);
+ // break;
+ //}
if (ratio>factor){
// if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex); // arr->RemoveAt(j);
Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
Int_t nseed = arr->GetEntriesFast();
for (Int_t i=0;i<nseed;i++)
fSeeds->AddLast(arr->RemoveAt(i));
- }
- // fSeeds = MakeSeedsSectors(0,fkNOS);
+ }
+ // fSeeds = MakeSeedsSectors(0,fkNOS);
}
TObjArray * AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
MakeSeeds(arr, sec, nup-1, nup-1-gap);
MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
}
+ gap = Int_t(0.3* nrows);
+ for (Int_t sec=sec1; sec<sec2;sec++){
+ //find secondaries
+ MakeSeeds2(arr, sec, nup-1, nup-1-gap);
+ MakeSeeds2(arr, sec, nup-1-shift, nup-1-shift-gap);
+ MakeSeeds2(arr, sec, nup-1-2*shift, nup-1-2*shift-gap);
+ //MakeSeeds2(arr, sec, nup-1-3*shift, nup-1-3*shift-gap);
+ MakeSeeds2(arr, sec, 30, 0);
+ }
+
Int_t nseed=arr->GetEntriesFast();
gap=Int_t(0.3*nrows);
// continue seeds
continue;
}
delete arr->RemoveAt(i);
- }
+ }
+
//
//remove seeds which overlaps
RemoveOverlap(arr,0.6,1);
//, &t=*pt;
if (!pt) continue;
if ((pt->IsActive()) && pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
- //FollowBackProlongation(t,nup);
+ //pt->Reset();
+ //FollowBackProlongation(*pt,nup-1);
+ //if ( pt->GetNumberOfClusters() < pt->fNFoundable*0.5 || pt->GetNumberOfClusters()<10 )
+ //delete arr->RemoveAt(i);
+ //else
+ // pt->Reset();
continue;
}
delete arr->RemoveAt(i);
- }
+ }
+ //RemoveOverlap(arr,0.6,1);
return arr;
}
}
}
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
+ //-----------------------------------------------------------------
+ // This function creates track seeds - without vertex constraint
+ //-----------------------------------------------------------------
+
+ Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
+ // Double_t cs=cos(alpha), sn=sin(alpha);
+ Int_t row0 = (i1+i2)/2;
+ Int_t drow = (i1-i2)/2;
+ const AliTPCRow& kr0=fSectors[sec][row0];
+ const AliTPCRow& krm=fSectors[sec][row0-1];
+ const AliTPCRow& krp=fSectors[sec][row0+1];
+ AliTPCRow * kr=0;
+
+ AliTPCpolyTrack polytrack;
+ Int_t nclusters=fSectors[sec][row0];
+
+ for (Int_t is=0; is < nclusters; is++) {
+ const AliTPCclusterMI * cl= kr0[is];
+ Double_t x = kr0.GetX();
+
+ // Initialization of the polytrack
+ polytrack.Reset();
+
+ Double_t y0= cl->GetY();
+ Double_t z0= cl->GetZ();
+ polytrack.AddPoint(x,y0,z0);
+ Float_t roady = 5*TMath::Sqrt(cl->GetSigmaY2()+0.2);
+ Float_t roadz = 5*TMath::Sqrt(cl->GetSigmaZ2()+0.2);
+ //
+ x = krm.GetX();
+ cl = krm.FindNearest(y0,z0,roady,roadz);
+ if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),roady,roadz);
+ //
+ x = krp.GetX();
+ cl = krp.FindNearest(y0,z0,roady,roadz);
+ if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),cl->GetSigmaY2()+0.05,cl->GetSigmaZ2()+0.05);
+ //
+ polytrack.UpdateParameters();
+ // follow polytrack
+ roadz = 0.6;
+ roady = 0.6;
+ //
+ Double_t yn,zn;
+ Int_t nfoundable = polytrack.GetN();
+ Int_t nfound = nfoundable;
+ for (Int_t ddrow = 2; ddrow<drow;ddrow++){
+ for (Int_t delta = -1;delta<=1;delta+=2){
+ Int_t row = row0+ddrow*delta;
+ kr = &(fSectors[sec][row]);
+ Double_t xn = kr->GetX();
+ Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone;
+ polytrack.GetFitPoint(xn,yn,zn);
+ if (TMath::Abs(yn)>ymax) continue;
+ nfoundable++;
+ AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
+ if (cln) {
+ polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),cln->GetSigmaY2()+0.05,cln->GetSigmaZ2()+0.05);
+ nfound++;
+ }
+ }
+ polytrack.UpdateParameters();
+ }
+ if ((nfound>0.5*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
+ // add polytrack candidate
+ Double_t x[5], c[15];
+ Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
+ polytrack.GetBoundaries(x3,x1);
+ x2 = (x1+x3)/2.;
+ polytrack.GetFitPoint(x1,y1,z1);
+ polytrack.GetFitPoint(x2,y2,z2);
+ polytrack.GetFitPoint(x3,y3,z3);
+ //
+ //is track pointing to the vertex ?
+ Double_t x0,y0,z0;
+ x0=0;
+ polytrack.GetFitPoint(x0,y0,z0);
+ if ( (TMath::Abs(z0-GetZ())<10) && (TMath::Abs(y0-GetY())<5)){ //if yes apply vertex constraint
+ x3 = 0;
+ y3 = GetY();
+ z3 = GetZ();
+ }
+
+ x[0]=y1;
+ x[1]=z1;
+ x[4]=f1(x1,y1,x2,y2,x3,y3);
+ if (TMath::Abs(x[4]) >= 0.0066) continue;
+ x[2]=f2(x1,y1,x2,y2,x3,y3);
+ //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
+ x[3]=f3(x1,y1,x2,y2,z1,z2);
+ if (TMath::Abs(x[3]) > 1.2) continue;
+ if (TMath::Abs(x[2]) > 0.99) continue;
+ // Double_t a=asin(x[2]);
+
+
+ Double_t sy=1.5, sz=1.5;
+ Double_t sy1=1.5, sz1=1.5;
+ Double_t sy2=1.3, sz2=1.3;
+ Double_t sy3=1.5;
+ //sz3=1.5;
+ //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
+ // Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
+ //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
+
+ Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
+ Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
+ Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
+ Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
+ Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
+ Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
+ Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
+ Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
+ Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
+ Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
+
+ c[0]=sy1;
+ c[1]=0.; c[2]=sz1;
+ c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
+ c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
+ c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
+ c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
+ c[13]=f30*sy1*f40+f32*sy2*f42;
+ c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
+
+ UInt_t index=0;
+ //kr0.GetIndex(is);
+ AliTPCseed *track=new AliTPCseed(index, x, c, x1, sec*alpha+shift);
+ track->fStopped =kFALSE;
+ track->fIsSeeding = kTRUE;
+ Int_t rc=FollowProlongation(*track, i2);
+ track->fLastPoint = i1; // first cluster in track position
+ if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
+ else arr->AddLast(track);
+ }
+ }
+
+
+
+}
+
+
+
+
+
//_____________________________________________________________________________
Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
//-----------------------------------------------------------------
Int_t gap=Int_t(0.3*nrows);
Int_t i;
+ //RemoveOverlap(fSeeds,0.6,2);
Int_t nseed=fSeeds->GetEntriesFast();
// outer sectors parallel tracking
ParallelTracking(fSectors->GetNRows()-gap-1,0);
+ //ParallelTracking(fSectors->GetNRows()-1,0);
+ //RemoveOverlap(fSeeds, 0.6,3);
+ // ParallelTracking(49,0);
printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
RemoveOverlap(fSeeds, 0.6,3);
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
if (!pt) continue;
Int_t nc=t.GetNumberOfClusters();
+ if (nc<20) continue;
t.CookdEdx(0.02,0.6);
CookLabel(pt,0.1); //For comparison only
// if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
- if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (nc>Int_t(0.3*nrows)))){
+ if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (t.fNFoundable>Int_t(0.3*nrows)))){
iotrack=pt;
tracktree.Fill();
cerr<<found++<<'\r';
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
if (!pt) continue;
if (!pt->IsActive()) continue;
+ if (pt->fRelativeSector>17) {
+ continue;
+ }
UpdateClusters(t,i,nr);
}
// prolonagate to the nearest cluster - if founded
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
if (!pt) continue;
- if (!pt->IsActive()) continue;
+ if (!pt->IsActive()) continue;
+ if (pt->fRelativeSector>17) {
+ continue;
+ }
FollowToNextCluster(i,nr);
}
// for (Int_t i= 0;i<fN;i++)
+//___________________________________________________________________
+AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
+ //-----------------------------------------------------------------------
+ // Return the index of the nearest cluster in z y
+ //-----------------------------------------------------------------------
+ Float_t maxdistance = roady*roady + roadz*roadz;
+
+ AliTPCclusterMI *cl =0;
+ for (Int_t i=Find(z-roadz); i<fN; i++) {
+ AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
+ if (c->GetZ() > z+roadz) break;
+ if ( (c->GetY()-y) > roady ) continue;
+ Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
+ if (maxdistance>distance) {
+ maxdistance = distance;
+ cl=c;
+ }
+ }
+ return cl;
+}
+
+
+
+
+
AliTPCseed::AliTPCseed():AliTPCtrack(){
//
fRow=0;
// TMath::Abs( Int_t(iz) - iz + 0.5);
//ampc *= 1.15*(1-0.3*dy);
//ampc *= 1.15*(1-0.3*dz);
-
+ // Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
+ //ampc *=zfactor;
}
else{
ampc = 1.0*point->GetCPoint().GetMax();
//ampc *= 1.15*(1-0.3*dy);
//ampc *= 1.15*(1-0.3*dz);
+ // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
+ //ampc *=zfactor;
}
ampc *= 2.0; // put mean value to channel 50
//ampc *= 0.58; // put mean value to channel 50
Float_t w = 1.;
// if (type>0) w = 1./(type/2.-0.5);
+ Float_t z = TMath::Abs(point->GetCPoint().GetZ());
if (i<64) {
ampc /= 0.6;
- }
- if (i>128){
- ampc /=1.5;
- }
+ //ampc /= (1+0.0008*z);
+ } else
+ if (i>128){
+ ampc /=1.5;
+ //ampc /= (1+0.0008*z);
+ }else{
+ //ampc /= (1+0.0008*z);
+ }
+
if (type<0) { //amp at the border - lower weight
// w*= 2.;
+
continue;
}
if (rsigma>1.5) ampc/=1.3; // if big backround