]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCTracklet.cxx
bug fix in the vertex selection
[u/mrichter/AliRoot.git] / TPC / AliTPCTracklet.cxx
index 51a2b059ff7f5213cdeb2524ef893d962b27034a..9f9f4e86b834a3bea8a4b7ca646929f225f55038 100755 (executable)
@@ -182,6 +182,9 @@ AliTPCTracklet::~AliTPCTracklet() {
   delete fPrimary;
 }
 
+
+
+/*
 void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
   //
   // Fit using Kalman filter
@@ -192,18 +195,29 @@ void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
     return;
   }
   // fit from inner to outer row
+  Double_t covar[15];
+  for (Int_t i=0;i<15;i++) covar[i]=0;
+  covar[0]=10.*10.;
+  covar[2]=10.*10.;
+  covar[5]=10.*10./(64.*64.);
+  covar[9]=10.*10./(64.*64.);
+  covar[14]=1*1;
+
+  //
   AliTPCseed *outerSeed=new AliTPCseed(*t);
   Int_t n=0;
   for (Int_t i=0;i<160;++i) {
+    
     AliTPCclusterMI *c=t->GetClusterPointer(i);
     if (c && RejectCluster(c,outerSeed)) continue;
     if (c&&c->GetDetector()==sector) {
       if (n==1)        {
        outerSeed->ResetCovariance(100.);
+       outerSeed->AddCovariance(covar);
       }
       ++n;
       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
-      Double_t cov[3]={0.1,0.,0.1}; //TODO: correct error parametrisation
+      Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
       if (!outerSeed->PropagateTo(r[0]) ||
          !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) {
        delete outerSeed;
@@ -214,9 +228,13 @@ void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
   }
   if (outerSeed)
     fOuter=new AliExternalTrackParam(*outerSeed);
-  delete outerSeed;
   // fit from outer to inner rows
-  AliTPCseed *innerSeed=new AliTPCseed(*t);
+  //  AliTPCseed *innerSeed=new AliTPCseed(*t);
+  AliTPCseed *innerSeed=0;
+  if (fOuter) innerSeed=new AliTPCseed(*outerSeed);
+  if (!innerSeed) innerSeed=new AliTPCseed(*t);
+  delete outerSeed;
+
   n=0;
   for (Int_t i=159;i>=0;--i) {
     AliTPCclusterMI *c=t->GetClusterPointer(i);
@@ -224,10 +242,11 @@ void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
     if (c&&c->GetDetector()==sector) {
       if (n==1)        {
        innerSeed->ResetCovariance(100.);
+       innerSeed->AddCovariance(covar);
       }
       ++n;
       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
-      Double_t cov[3]={0.1,0.,0.1};
+      Double_t cov[3]={0.01,0.,0.01};
       if (!innerSeed->PropagateTo(r[0]) ||
          !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) {
        delete innerSeed;
@@ -257,6 +276,137 @@ void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
 
   delete t;
 }
+*/
+
+
+void AliTPCTracklet::FitKalman(const AliTPCseed *seed,Int_t sector) {
+  //
+  // Fit using Kalman filter
+  //
+  AliTPCseed *track=new AliTPCseed(*seed);
+  if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
+    delete track;
+    return;
+  }
+  // fit from inner to outer row
+  Double_t covar[15];
+  for (Int_t i=0;i<15;i++) covar[i]=0;
+  covar[0]=10.*10.;
+  covar[2]=10.*10.;
+  covar[5]=10.*10./(64.*64.);
+  covar[9]=10.*10./(64.*64.);
+  covar[14]=1*1;
+  Float_t xmin=1000, xmax=-10000;
+  Int_t imin=158, imax=0;
+  for (Int_t i=0;i<160;i++) {
+    AliTPCclusterMI *c=track->GetClusterPointer(i);
+    if (!c) continue;
+    if (c->GetDetector()!=sector)  continue;
+    if (c->GetX()<xmin) xmin=c->GetX();
+    if (c->GetX()>xmax) xmax=c->GetX();
+    if (i<imin) imin=i;
+    if (i>imax) imax=i;
+  }
+  if(imax-imin<10) {
+    delete track;
+    return;
+  }
+
+  for (Float_t x=track->GetX(); x<xmin; x++) track->PropagateTo(x);
+  track->AddCovariance(covar);
+  //
+  AliExternalTrackParam paramIn;
+  AliExternalTrackParam paramOut;
+  Bool_t isOK=kTRUE;
+  //
+  //
+  //
+  for (Int_t i=imin; i<=imax; i++){
+    AliTPCclusterMI *c=track->GetClusterPointer(i);
+    if (!c) continue;
+    Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
+    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
+    AliTPCseed::GetError(c, track,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    if (!track->PropagateTo(r[0])) {
+      isOK=kFALSE;
+      break;
+    }
+    if (RejectCluster(c,track)) continue;
+    if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
+  }
+  if (!isOK) { delete track; return;}
+  track->AddCovariance(covar);
+  //
+  //
+  //
+  for (Int_t i=imax; i>=imin; i--){
+    AliTPCclusterMI *c=track->GetClusterPointer(i);
+    if (!c) continue;
+    Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
+    Double_t cov[3]={0.01,0.,0.01}; 
+    AliTPCseed::GetError(c, track,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    if (!track->PropagateTo(r[0])) {
+      isOK=kFALSE;
+      break;
+    }
+    if (RejectCluster(c,track)) continue;
+    if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
+  }
+  if (!isOK) { delete track; return;}
+  paramIn = *track;
+  track->AddCovariance(covar);
+  //
+  //
+  for (Int_t i=imin; i<=imax; i++){
+    AliTPCclusterMI *c=track->GetClusterPointer(i);
+    if (!c) continue;
+    Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
+    Double_t cov[3]={0.01,0.,0.01}; 
+    AliTPCseed::GetError(c, track,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    if (!track->PropagateTo(r[0])) {
+      isOK=kFALSE;
+      break;
+    }
+    if (RejectCluster(c,track)) continue;
+    if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
+  }
+  if (!isOK) { delete track; return;}
+  paramOut=*track;
+  //
+  //
+  //
+  fOuter=new AliExternalTrackParam(paramOut);
+  fInner=new AliExternalTrackParam(paramIn);
+  //
+
+
+ //  // propagate to the primary vertex
+//   if (fInner) {
+//     AliExternalTrackParam  *param= new AliExternalTrackParam(*fInner);
+//     Double_t pos[]={0.,0.,0.};
+//     Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
+//     AliESDVertex vertex(pos,sigma);
+//     if (param->PropagateToVertex(&vertex))
+//       fPrimary=new AliExternalTrackParam(*param);
+//     delete param;
+//     // for better comparison one does not want to have alpha changed...
+//     if (fPrimary) if (!fPrimary->Rotate(fInner->GetAlpha())) {
+//       delete fPrimary;
+//       fPrimary=0;
+//     }
+//   }
+
+  delete track;
+}
+
+
+
 
 void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
                               TrackType type) {
@@ -467,7 +617,7 @@ Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
   Double_t dx=x0-xr0;
   if (dx*dx>=R*R) return kFALSE;
-  Double_t dy=TMath::Sqrt(R*R-dx*dx); //sign!!
+  Double_t dy=TMath::Sqrt((R-dx)*(R+dx)); //sign!!
   if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
     dy=-dy;
   Double_t y0=yr0+dy; 
@@ -569,19 +719,25 @@ Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
       t2m=new AliExternalTrackParam(*t2.GetOuter());
     }
     Double_t mx=.5*(t1m->GetX()+t2m->GetX());
-    Double_t b1,b2;
+    //Double_t b1,b2;
     Double_t xyz[3];
     t1m->GetXYZ(xyz);
-    b1=GetBz(xyz);
+    //b1=GetBz(xyz);
+    Double_t b1[3]; AliTracker::GetBxByBz(xyz,b1);
     t2m->GetXYZ(xyz);
-    b2=GetBz(xyz);
+    //b2=GetBz(xyz);
+    Double_t b2[3]; AliTracker::GetBxByBz(xyz,b2);
     if (t1m->Rotate(t2m->GetAlpha()) 
-       && t1m->PropagateTo(mx,b1) 
-       && t2m->PropagateTo(mx,b2));
+       //&& t1m->PropagateTo(mx,b1) 
+       //&& t2m->PropagateTo(mx,b2));
+       && t1m->PropagateToBxByBz(mx,b1) 
+       && t2m->PropagateToBxByBz(mx,b2));
     else
       if (t2m->Rotate(t1m->GetAlpha())
-         && t1m->PropagateTo(mx,b1) 
-         && t2m->PropagateTo(mx,b2));
+         //&& t1m->PropagateTo(mx,b1) 
+         //&& t2m->PropagateTo(mx,b2));
+         && t1m->PropagateToBxByBz(mx,b1) 
+         && t2m->PropagateToBxByBz(mx,b2));
       else {
        delete t1m;
        delete t2m;
@@ -594,11 +750,9 @@ Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
   return t1m&&t2m;
 }
 
-double AliTPCTracklet::GetBz(Double_t *xyz) {
-  if (AliTracker::UniformField())
-    return AliTracker::GetBz();
-  else
-    return AliTracker::GetBz(xyz);
+double AliTPCTracklet::GetBz(Double_t *xyz) 
+{
+  return AliTracker::GetBz(xyz);
 }
 
 void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,