reco update: methods for smoothing, track hypotheses container added
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Jan 2013 18:56:06 +0000 (18:56 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Jan 2013 18:56:06 +0000 (18:56 +0000)
ITS/UPGRADE/AliITSURecoParam.cxx
ITS/UPGRADE/AliITSUSeed.cxx
ITS/UPGRADE/AliITSUSeed.h
ITS/UPGRADE/AliITSUTrackHyp.cxx [new file with mode: 0644]
ITS/UPGRADE/AliITSUTrackHyp.h [new file with mode: 0644]
ITS/UPGRADE/AliITSUTrackerGlo.cxx
ITS/UPGRADE/AliITSUTrackerGlo.h
ITS/UPGRADE/CMakelibITSUpgradeRec.pkg
ITS/UPGRADE/ITSUpgradeRecLinkDef.h

index 175b796..be2da0d 100644 (file)
@@ -36,7 +36,7 @@ const Double_t AliITSURecoParam::fgkSigmaRoadY                    = 1.;//1000e-4
 const Double_t AliITSURecoParam::fgkSigmaRoadZ                    = 1.;//1000e-4;
 const Double_t AliITSURecoParam::fgkMaxTr2ClChi2                  = 15.;
 const Double_t AliITSURecoParam::fgkTanLorentzAngle               = 0;
-const Double_t AliITSURecoParam::fgkMissPenalty                   = 3.0;
+const Double_t AliITSURecoParam::fgkMissPenalty                   = 2.0;
 //
 // hardwired params for TPC-ITS border layer
 const Double_t AliITSURecoParam::fgkTPCITSWallRMin                = 50.;
index 2907b20..66b3bad 100644 (file)
@@ -1,6 +1,7 @@
 #include <TString.h>
 #include <TMath.h>
 #include "AliITSUSeed.h"
+#include "AliLog.h"
 using namespace TMath;
 
 ClassImp(AliITSUSeed)
@@ -11,6 +12,7 @@ AliITSUSeed::AliITSUSeed()
   ,fClID(0)
   ,fChi2Glo(0)
   ,fChi2Cl(0)
+  ,fChi2Penalty(0)
   ,fParent(0)
 {
   // def c-tor
@@ -20,7 +22,7 @@ AliITSUSeed::AliITSUSeed()
 //_________________________________________________________________________
 AliITSUSeed::~AliITSUSeed()
 {
-  // d-rot
+  // d-tor
 }
 
 //_________________________________________________________________________
@@ -30,11 +32,13 @@ AliITSUSeed::AliITSUSeed(const AliITSUSeed& src)
   ,fClID(src.fClID)
   ,fChi2Glo(src.fChi2Glo)
   ,fChi2Cl(src.fChi2Cl)
+  ,fChi2Penalty(src.fChi2Penalty)
   ,fParent(src.fParent) 
 {
   // def c-tor
   for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
-  for (int i=kNBElem;i--;) fBMatrix[i] = src.fBMatrix[i];
+  for (int i=kNKElem;i--;) fKMatrix[i] = src.fKMatrix[i];
+  for (int i=kNRElem;i--;) fRMatrix[i] = src.fRMatrix[i];
   fResid[0]=src.fResid[0];
   fResid[1]=src.fResid[1];  
   fCovIYZ[0]=src.fCovIYZ[0];
@@ -58,11 +62,21 @@ void AliITSUSeed::Print(Option_t* opt) const
 {
   // print seed info
   int lr,cl = GetLrCluster(lr);
-  printf("%cLr%d Cl:%4d Chi2Glo:%7.2f(%7.2f) Chi2Cl:",IsKilled() ? '-':' ',
-        lr,cl,GetChi2Glo(),GetChi2GloNrm());
-  cl<0 ? printf("   NA  ") : printf("%7.2f",GetChi2Cl());
+  printf("%cLr%d Cl:%4d Chi2Glo:%7.2f(%7.2f) Chi2Cl:%7.2f Penalty: %7.2f",IsKilled() ? '-':' ',
+        lr,cl,GetChi2Glo(),GetChi2GloNrm(),GetChi2Cl(), GetChi2Penalty());
   printf(" |"); 
-  for (int i=0;i<=12;i++) printf("%c",HasClusterOnLayer(i) ? '+':'-'); printf("|\n");
+  int lrc=0;
+  const AliITSUSeed *sdc = this;
+  while(1) {
+    if (lrc<lr) printf(".");
+    else {
+      sdc = sdc->GetParent(lrc);
+      if (!sdc) break;
+      printf("%c",sdc->GetClusterID()<0 ? '.': (sdc->IsFake() ? '-':'+')); 
+    }
+    lrc++;
+  }
+  printf("|\n");
   TString opts = opt; opts.ToLower();
   if (opts.Contains("etp")) AliExternalTrackParam::Print();
   if (opts.Contains("parent") && GetParent()) GetParent()->Print(opt);
@@ -72,7 +86,7 @@ void AliITSUSeed::Print(Option_t* opt) const
 Float_t AliITSUSeed::GetChi2GloNrm() const
 {
   int ndf = 2*GetNLayersHit() - 5;
-  return ndf>0 ? fChi2Glo/ndf : fChi2Glo;
+  return (ndf>0 ? fChi2Glo/ndf : fChi2Glo) + fChi2Penalty;
 }
 
 
@@ -85,8 +99,11 @@ Int_t AliITSUSeed::Compare(const TObject* obj)  const
   if (!IsKilled() && sd->IsKilled()) return -1;
   if ( IsKilled() &&!sd->IsKilled()) return  1;
   //
-  if      (GetChi2Glo()+kTol<sd->GetChi2Glo()) return -1;
-  else if (GetChi2Glo()-kTol>sd->GetChi2Glo()) return  1;
+  float chi2This  = GetChi2GloNrm();
+  float chi2Other = sd->GetChi2GloNrm();
+
+  if      (chi2This+kTol<chi2Other) return -1;
+  else if (chi2This-kTol>chi2Other) return  1;
   return 0;
 }
 
@@ -97,7 +114,7 @@ Bool_t AliITSUSeed::IsEqual(const TObject* obj)  const
   const AliITSUSeed* sd = (const AliITSUSeed*)obj;
   const Float_t kTol = 1e-5;
   if (IsKilled() != sd->IsKilled()) return kFALSE;
-  return Abs(GetChi2Glo() - sd->GetChi2Glo())<kTol;
+  return Abs(GetChi2GloNrm() - sd->GetChi2GloNrm())<kTol;
 }
 
 //______________________________________________________________________________
@@ -207,6 +224,79 @@ Bool_t AliITSUSeed::PropagateToX(Double_t xk, Double_t b)
 }
 
 //______________________________________________________________________________
+Bool_t AliITSUSeed::RotateToAlpha(Double_t alpha) 
+{
+  // Transform this track to the local coord. system rotated
+  // by angle "alpha" (rad) with respect to the global coord. system. 
+  //
+  if (TMath::Abs(fP[2]) >= kAlmost1) {
+     AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2])); 
+     return kFALSE;
+  }
+  //
+  if      (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
+  else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
+  //
+  Double_t &fP0=fP[0];
+  Double_t &fP2=fP[2];
+  Double_t &fC00=fC[0];
+  Double_t &fC10=fC[1];
+  Double_t &fC20=fC[3];
+  Double_t &fC21=fC[4];
+  Double_t &fC22=fC[5];
+  Double_t &fC30=fC[6];
+  Double_t &fC32=fC[8];
+  Double_t &fC40=fC[10];
+  Double_t &fC42=fC[12];
+  //
+  Double_t x=fX;
+  Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
+  Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
+  // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
+  // direction in local frame is along the X axis
+  if ((cf*ca+sf*sa)<0) {
+    AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
+    return kFALSE;
+  }
+  //
+  Double_t tmp=sf*ca - cf*sa;
+
+  if (TMath::Abs(tmp) >= kAlmost1) {
+     if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))  
+        AliWarning(Form("Rotation failed ! %.10e",tmp));
+     return kFALSE;
+  }
+  fAlpha = alpha;
+  fX =  x*ca + fP0*sa;
+  fP0= -x*sa + fP0*ca;
+  fP2=  tmp;
+
+  if (TMath::Abs(cf)<kAlmost0) {
+    AliError(Form("Too small cosine value %f",cf)); 
+    cf = kAlmost0;
+  } 
+
+  Double_t rr=(ca+sf/cf*sa);  
+
+  fC00 *= (ca*ca);
+  fC10 *= ca;
+  fC20 *= ca*rr;
+  fC21 *= rr;
+  fC22 *= rr*rr;
+  fC30 *= ca;
+  fC32 *= rr;
+  fC40 *= ca;
+  fC42 *= rr;
+  //
+  fRMatrix[kR00] = ca;
+  fRMatrix[kR22] = rr; 
+  //
+  CheckCovariance();
+
+  return kTRUE;
+}
+
+//______________________________________________________________________________
 Bool_t AliITSUSeed::GetTrackingXAtXAlpha(double xOther, double alpOther, double bz, double &xdst)
 {
   // calculate X and Y in the tracking frame of the track, corresponding to other X,Alpha tracking
@@ -214,9 +304,9 @@ Bool_t AliITSUSeed::GetTrackingXAtXAlpha(double xOther, double alpOther, double
   double &y=fP[0], &sf=fP[2], cf=Sqrt((1.-sf)*(1.+sf));
   double eta = xOther - fX*ca - y*sa;
   double xi  = sf*ca - cf*sa;
-  if (xi>= kAlmost1) return kFALSE;
+  if (Abs(xi)>= kAlmost1) return kFALSE;
   double nu  = xi + GetC(bz)*eta;
-  if (nu>= kAlmost1) return kFALSE;
+  if (Abs(nu)>= kAlmost1) return kFALSE;
   xdst = xOther*ca - sa*( y*ca-fX*sa + eta*(xi+nu)/(Sqrt((1.-xi)*(1.+xi)) + Sqrt((1.-nu)*(1.+nu))) );
   return kTRUE;
 }
@@ -247,28 +337,27 @@ Bool_t AliITSUSeed::Update()
 {
   // Update the track parameters with the measurement stored during GetPredictedChi2
   //
-  Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
-  Double_t 
-  &fC00=fC[0],
-  &fC10=fC[1],   &fC11=fC[2],  
-  &fC20=fC[3],   &fC21=fC[4],   &fC22=fC[5],
-  &fC30=fC[6],   &fC31=fC[7],   &fC32=fC[8],   &fC33=fC[9],  
-  &fC40=fC[10],  &fC41=fC[11],  &fC42=fC[12],  &fC43=fC[13], &fC44=fC[14];
+  Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4],
+    &fC00=fC[kS00],
+    &fC10=fC[kS10],  &fC11=fC[kS11],  
+    &fC20=fC[kS20],  &fC21=fC[kS21],  &fC22=fC[kS22],
+    &fC30=fC[kS30],  &fC31=fC[kS31],  &fC32=fC[kS32],  &fC33=fC[kS33],  
+    &fC40=fC[kS40],  &fC41=fC[kS41],  &fC42=fC[kS42],  &fC43=fC[kS43], &fC44=fC[kS44];
   //
   double &r00=fCovIYZ[0],&r01=fCovIYZ[1],&r11=fCovIYZ[2];
   double &dy=fResid[0], &dz=fResid[1];
   //
-  // store info needed for smoothing in the fBMatrix
-  double &k00 = fBMatrix[kB00] = fC00*r00+fC10*r01;
-  double &k01 = fBMatrix[kB01] = fC00*r01+fC10*r11;
-  double &k10 = fBMatrix[kB10] = fC10*r00+fC11*r01;
-  double &k11 = fBMatrix[kB11] = fC10*r01+fC11*r11;  
-  double &k20 = fBMatrix[kB20] = fC20*r00+fC21*r01;
-  double &k21 = fBMatrix[kB21] = fC20*r01+fC21*r11;
-  double &k30 = fBMatrix[kB30] = fC30*r00+fC31*r01;
-  double &k31 = fBMatrix[kB31] = fC30*r01+fC31*r11;
-  double &k40 = fBMatrix[kB40] = fC40*r00+fC41*r01;
-  double &k41 = fBMatrix[kB41] = fC40*r01+fC41*r11;
+  // store info needed for smoothing in the fKMatrix
+  double &k00 = fKMatrix[kK00] = fC00*r00+fC10*r01;
+  double &k01 = fKMatrix[kK01] = fC00*r01+fC10*r11;
+  double &k10 = fKMatrix[kK10] = fC10*r00+fC11*r01;
+  double &k11 = fKMatrix[kK11] = fC10*r01+fC11*r11;  
+  double &k20 = fKMatrix[kK20] = fC20*r00+fC21*r01;
+  double &k21 = fKMatrix[kK21] = fC20*r01+fC21*r11;
+  double &k30 = fKMatrix[kK30] = fC30*r00+fC31*r01;
+  double &k31 = fKMatrix[kK31] = fC30*r01+fC31*r11;
+  double &k40 = fKMatrix[kK40] = fC40*r00+fC41*r01;
+  double &k41 = fKMatrix[kK41] = fC40*r01+fC41*r11;
   //
   Double_t sf=fP2 + k20*dy + k21*dz;
   if (TMath::Abs(sf) > kAlmost1) return kFALSE;  
@@ -298,10 +387,316 @@ Bool_t AliITSUSeed::Update()
   
   fC44-=k40*c04+k41*c14; 
   //
-  k00 -= 1;
-  k11 -= 1;
-  //
   CheckCovariance();
+  //
+  return kTRUE;
+}
+
+
+//____________________________________________________________________
+Bool_t AliITSUSeed::Smooth(Double_t vecL[5],Double_t matL[15]) 
+{
+  // Prepare MBF smoothing auxiliary params for smoothing at prev. point:
+  // \hat{l_N} = 0
+  // \hat{L_N} = 0
+  // \tilde{l_j} = -H_j^T N_{j}^{-1} z_j + B_{j}^T \hat{l_j}
+  // \tilde{L_j} =  H_j^T N_{j}^{-1} H_j + B_j^T \hat{L_j} B_j
+  // \hat{l_j} = F_j^T \tilde{l_{j+1}}
+  // \hat{L_j} = F_j^T \tilde{L_{j+1}} F_j
+  //
+  // P_{j/N} = P_{j/j} - P_{j/j} \hat{L_j} P_{j/j}
+  // \hat{x_{j/N}} = \hat{x_{j/j}} - P_{j/j} \hat{l_j}
+  //
+  // N^-1 = fCovIYZ
+  // z = fResid
+  // B = I - K H
+  // H = {{1,0,0,0,0},{0,1,0,0,0}}
+  // 
+  // calc. \tilde{l_j} 
+  //
+  if (GetClusterID()<0) return kTRUE;
+  //
+
+  double 
+    &k00=fKMatrix[kK00],&k01=fKMatrix[kK01],
+    &k10=fKMatrix[kK10],&k11=fKMatrix[kK11],
+    &k20=fKMatrix[kK20],&k21=fKMatrix[kK21],
+    &k30=fKMatrix[kK30],&k31=fKMatrix[kK31],
+    &k40=fKMatrix[kK40],&k41=fKMatrix[kK41];
+  double 
+    &l00=matL[kS00],
+    &l10=matL[kS10], &l11=matL[kS11],  
+    &l20=matL[kS20], &l21=matL[kS21], &l22=matL[kS22],
+    &l30=matL[kS30], &l31=matL[kS31], &l32=matL[kS32], &l33=matL[kS33],  
+    &l40=matL[kS40], &l41=matL[kS41], &l42=matL[kS42], &l43=matL[kS43], &l44=matL[kS44];
+  //
+  // calculate correction
+  double corrVec[5]={0},corrMat[15]={0};
+  corrVec[0] = fC[kS00]*vecL[0] + fC[kS10]*vecL[1] + fC[kS20]*vecL[2] + fC[kS30]*vecL[3] + fC[kS40]*vecL[4]; 
+  corrVec[1] = fC[kS10]*vecL[0] + fC[kS11]*vecL[1] + fC[kS21]*vecL[2] + fC[kS31]*vecL[3] + fC[kS41]*vecL[4]; 
+  corrVec[2] = fC[kS20]*vecL[0] + fC[kS21]*vecL[1] + fC[kS22]*vecL[2] + fC[kS32]*vecL[3] + fC[kS42]*vecL[4]; 
+  corrVec[3] = fC[kS30]*vecL[0] + fC[kS31]*vecL[1] + fC[kS32]*vecL[2] + fC[kS33]*vecL[3] + fC[kS43]*vecL[4]; 
+  corrVec[4] = fC[kS40]*vecL[0] + fC[kS41]*vecL[1] + fC[kS42]*vecL[2] + fC[kS43]*vecL[3] + fC[kS44]*vecL[4]; 
+  //
+  double *crm = ProdABA(fC,matL);
+  for (int i=0;i<15;i++) corrMat[i] = crm[i];
+
+  double vcL0 = vecL[0], vcL1 = vecL[1];
+  vecL[0] -= k00*vcL0+k10*vcL1+k20*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
+  vecL[1] -= k01*vcL0+k11*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
+
+  /*
+  double vcL0 = vecL[0], vcL1 = vecL[1];
+  vecL[0] -= k00*vcL0+k10*vcL1+fKMatrix[kK20]*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
+  vecL[1] -= k01*vcL0+fKMatrix[kK11]*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
+  vecL[3] += fFMatrix[kF13]*vecL[1]; 
+  vecL[4]  = fFMatrix[kF04]*vecL[0] + fFMatrix[kF14]*vecL[1] + fFMatrix[kF24]*vecL[2] + fFMatrix[kF44]*vecL[4];
+  vecL[2] += fFMatrix[kF02]*vecL[0] + fFMatrix[kF12]*vecL[1];
+  //
+  */
+  // and \hat{l_j} in one go
+
+  // L = H^T * sg * H + (I-KH)^T * L * (I - KH)
+  double v00 =  k00*l00+k10*l10+k20*l20+k30*l30+k40*l40;
+  double v10 =  k00*l10+k10*l11+k20*l21+k30*l31+k40*l41;
+  double v20 =  k00*l20+k10*l21+k20*l22+k30*l32+k40*l42;
+  double v30 =  k00*l30+k10*l31+k20*l32+k30*l33+k40*l43;
+  double v40 =  k00*l40+k10*l41+k20*l42+k30*l43+k40*l44;
+  //
+  double v01 =  k01*l00+k11*l10+k21*l20+k31*l30+k41*l40;
+  double v11 =  k01*l10+k11*l11+k21*l21+k31*l31+k41*l41;
+  double v21 =  k01*l20+k11*l21+k21*l22+k31*l32+k41*l42;
+  double v31 =  k01*l30+k11*l31+k21*l32+k31*l33+k41*l43;
+  double v41 =  k01*l40+k11*l41+k21*l42+k31*l43+k41*l44;
+  //
+  // (H^T * K^T * L * K * H) - (L * K * H) - (H^T * K^T * L) + (H^T*N^-1*H)
+  l00 += k00*v00 + k10*v10 + k20*v20 + k30*v30 + k40*v40 - v00 - v00 + fCovIYZ[0];
+  l10 += k01*v00 + k11*v10 + k21*v20 + k31*v30 + k41*v40 - v01 - v10 + fCovIYZ[1];
+  l11 += k01*v01 + k11*v11 + k21*v21 + k31*v31 + k41*v41 - v11 - v11 + fCovIYZ[2];
+  //
+  l20 -= v20;
+  l21 -= v21;
+  l30 -= v30;
+  l31 -= v31;
+  l40 -= v40;
+  l41 -= v41;
+  //
+  printf("CorrMt:\n");
+  printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
+        corrMat[kS00],corrMat[kS10],corrMat[kS11],corrMat[kS20],corrMat[kS21],corrMat[kS22],
+        corrMat[kS30],corrMat[kS31],corrMat[kS32],corrMat[kS33],
+        corrMat[kS40],corrMat[kS41],corrMat[kS42],corrMat[kS43],corrMat[kS44]);
+  
+  printf("SMcorr: %+e %+e %+e %+e %+e\n",corrVec[0],corrVec[1],corrVec[2],corrVec[3],corrVec[4]);
+
+  printf("State : "); this->AliExternalTrackParam::Print("");
+  //
+  printf("\nBefore transport back (RotElems: %+e %+e)\n",fRMatrix[kR00],fRMatrix[kR22]);
+  printf("Res: %+e %+e | Err: %+e %+e %+e\n",fResid[0],fResid[1],fCovIYZ[0],fCovIYZ[1],fCovIYZ[2]);
+  printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
+  //
+  printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
+        matL[kS00],matL[kS10],matL[kS11],matL[kS20],matL[kS21],matL[kS22],
+        matL[kS30],matL[kS31],matL[kS32],matL[kS33],matL[kS40],matL[kS41],matL[kS42],matL[kS43],matL[kS44]);
+  //
+  printf("F: "); for (int i=0;i<kNFElem;i++) printf("%+e ",fFMatrix[i]); printf("\n");  
+  printf("K: "); for (int i=0;i<kNKElem;i++) printf("%+e ",fKMatrix[i]); printf("\n");  
+  //
+  // apply rotation matrix (diagonal)
+  vecL[0] *= fRMatrix[kR00];
+  vecL[2] *= fRMatrix[kR22];
+  //
+  l00 *= fRMatrix[kR00]*fRMatrix[kR00];
+  l10 *= fRMatrix[kR00];
+  l20 *= fRMatrix[kR22]*fRMatrix[kR00];
+  l21 *= fRMatrix[kR22];
+  l22 *= fRMatrix[kR22]*fRMatrix[kR22];
+  l30 *= fRMatrix[kR00];
+  l32 *= fRMatrix[kR22];
+  l40 *= fRMatrix[kR00];
+  l42 *= fRMatrix[kR22];
+  //
+  // Apply translation matrix F^T. Note, that fFMatrix keeps non-trivial elems of F-1 = f, except the e-loss coeff f44
+  // We need F^T*L* F = L + (L*f) + (L*f)^T + f^T * (L*f)
+  //
+  double 
+    &f02=fFMatrix[kF02],&f04=fFMatrix[kF04],
+    &f12=fFMatrix[kF12],&f13=fFMatrix[kF13],&f14=fFMatrix[kF14],
+    &f24=fFMatrix[kF24],
+    f44 =fFMatrix[kF44];
+  //
+  vecL[4]  = f04*vecL[0]+f14*vecL[1]+f24*vecL[2]+f44*vecL[4];
+  vecL[3] += f13*vecL[1];
+  vecL[2] += f02*vecL[0]+f12*vecL[1];
+  //
+  f44 -= 1.0; // !!!!!
+  //
+  //b = L*f
+  Double_t b02=l00*f02+l10*f12, b03=l10*f13, b04=l00*f04+l10*f14+l20*f24+l40*f44;
+  Double_t b12=l10*f02+l11*f12, b13=l11*f13, b14=l10*f04+l11*f14+l21*f24+l41*f44;
+  Double_t b22=l20*f02+l21*f12, b23=l21*f13, b24=l20*f04+l21*f14+l22*f24+l42*f44;
+  Double_t b32=l30*f02+l31*f12, b33=l31*f13, b34=l30*f04+l31*f14+l32*f24+l43*f44;
+  Double_t b42=l40*f02+l41*f12, b43=l41*f13, b44=l40*f04+l41*f14+l42*f24+l44*f44;
+  //
+  //a = f^T * b = f^T * L * f, profit from symmetry
+  Double_t a22=f02*b02+f12*b12, a33=f13*b13, a44=f04*b04+f14*b14+f24*b24+f44*b44,
+    a32=f13*b12, //= a23=f02*b03+f12*b13, 
+    a42=f02*b04+f12*b14, //f04*b02+f14*b12+f24*b22+f44*b42 = a24
+    a43=f13*b14;         //f04*b03+f14*b13+f24*b23+f44*b43 = a34
+  //
+  // F^T*L* F = L + (b + b^T + a)
+  l44 += b44 + b44 + a44;
+  l43 += b43 + b34 + a43;
+  l42 += b42 + b24 + a42;
+  l41 += b14;
+  l40 += b04;
+  l33 += b33 + b33 + a33;
+  l32 += b32 + b23 + a32;
+  l31 += b13;
+  l30 += b03;
+  l22 += b22 + b23 + a22;
+  l21 += b12;
+  l20 += b02;
+  //
+  printf("After transport back\n");
+  printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
+  //
+  printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
+        matL[kS00],matL[kS10],matL[kS11],matL[kS20],matL[kS21],matL[kS22],
+        matL[kS30],matL[kS31],matL[kS32],matL[kS33],matL[kS40],matL[kS41],matL[kS42],matL[kS43],matL[kS44]);
+
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Double_t* AliITSUSeed::ProdABA(const double a[15],const double b[15]) const
+{
+  // product of symmetric matrices A*B*A
+  //
+  const Short_t knd[5][5] = {
+    {kS00,kS10,kS20,kS30,kS40},
+    {kS10,kS11,kS21,kS31,kS41},
+    {kS20,kS21,kS22,kS32,kS42},
+    {kS30,kS31,kS32,kS33,kS43},
+    {kS40,kS41,kS42,kS43,kS44}
+  };
+  //
+  static double aba[15];
+  // 1) ba = B*A
+  double ba[5][5];
+  for (int i=5;i--;) for (int j=5;j--;) {
+      ba[i][j] = 0;
+      for (int k=5;k--;) ba[i][j] += b[knd[i][k]]*a[knd[k][j]];
+    }
+  //
+  // 2) A * ba, lower triangle only
+  for (int i=5;i--;) for (int j=i+1;j--;) {
+      aba[knd[i][j]] = 0;
+      for (int k=5;k--;) aba[knd[i][j]] += a[knd[i][k]]*ba[k][j];
+    }
+  //
+  return &aba[0];
+}
 
+/*
+//____________________________________________________________________
+Bool_t AliITSUSeed::Smooth(Double_t vecL[5],Double_t matL[15]) 
+{
+  // Prepare MBF smoothing auxiliary params for smoothing at prev. point:
+  // \hat{l_N} = 0
+  // \hat{L_N} = 0
+  // \tilde{l_j} = -H_j^T N_{j}^{-1} z_j + B_{j}^T \hat{l_j}
+  // \tilde{L_j} =  H_j^T N_{j}^{-1} H_j + B_j^T \hat{L_j} B_j
+  // \hat{l_j} = F_j^T \tilde{l_{j+1}}
+  // \hat{L_j} = F_j^T \tilde{L_{j+1}} F_j
+  //
+  // P_{j/N} = P_{j/j} - P_{j/j} \hat{L_j} P_{j/j}
+  // \hat{x_{j/N}} = \hat{x_{j/j}} - P_{j/j} \hat{l_j}
+  //
+  // N^-1 = fCovIYZ
+  // z = fResid
+  // B = I - K H
+  // H = {{1,0,0,0,0},{0,1,0,0,0}}
+  // 
+  // calc. \tilde{l_j} and \hat{l_j} in one go
+  //
+  if (GetClusterID()<0) return kTRUE;
+  //
+  double 
+    &k00=fKMatrix[kK00],&k01=fKMatrix[kK01],
+    &k10=fKMatrix[kK10],&k11=fKMatrix[kK11],
+    &k20=fKMatrix[kK20],&k21=fKMatrix[kK21],
+    &k30=fKMatrix[kK30],&k31=fKMatrix[kK31],
+    &k40=fKMatrix[kK40],&k41=fKMatrix[kK41];
+  double 
+    &matL00=matL[kS00],
+    &matL10=matL[kS01],  &matL11=matL[kS11],  
+    &matL20=matL[kS20],  &matL21=matL[kS21],  &matL22=matL[kS22],
+    &matL30=matL[kS30],  &matL31=matL[kS31],  &matL32=matL[kS32],  &matL33=matL[kS33],  
+    &matL40=matL[kS40],  &matL41=matL[kS41],  &matL42=matL[kS42],  &matL43=matL[kS43], &matL44=matL[kS44];
+  //
+  double vcL0 = vecL[0], vcL1 = vecL[1];
+  vecL[0] -= k00*vcL0+k10*vcL1+fKMatrix[kK20]*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
+  vecL[1] -= k01*vcL0+fKMatrix[kK11]*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
+  vecL[3] += fFMatrix[kF13]*vecL[1]; 
+  vecL[4] =  fFMatrix[kF04]*vecL[0] + fFMatrix[kF14]*vecL[1] + fFMatrix[kF24]*vecL[2] + fFMatrix[kF44]*vecL[4];
+  vecL[2] += fFMatrix[kF02]*vecL[0] + fFMatrix[kF12]*vecL[1];
+  //
+
+  // L = H^T * sg * H + (I-KH)^T * L * (I - KH)
+  double v00 =  k00*matL00+k10*matL10+k20*matL20+k30*matL30+k40*matL40;
+  double v10 =  k00*matL10+k10*matL11+k20*matL21+k30*matL31+k40*matL41;
+  double v20 =  k00*matL20+k10*matL12+k20*matL22+k30*matL32+k40*matL42;
+  double v30 =  k00*matL30+k10*matL13+k20*matL23+k30*matL33+k40*matL43;
+  double v40 =  k00*matL40+k10*matL14+k20*matL24+k30*matL34+k40*matL44;
+  //
+  double v01 =  k01*matL00+k11*matL10+k21*matL20+k31*matL30+k41*matL40;
+  double v11 =  k01*matL01+k11*matL11+k21*matL21+k31*matL31+k41*matL41;
+  double v21 =  k01*matL02+k11*matL12+k21*matL22+k31*matL32+k41*matL42;
+  double v31 =  k01*matL03+k11*matL13+k21*matL23+k31*matL33+k41*matL43;
+  double v41 =  k01*matL04+k11*matL14+k21*matL24+k31*matL34+k41*matL44;
+  //
+  double t00 =  k00*matL00+k10*matL01+k20*matL02+k30*matL03+k40*matL04;
+  double t10 =  k00*matL10+k10*matL11+k20*matL12+k30*matL13+k40*matL14;
+  double t20 =  k00*matL20+k10*matL21+k20*matL22+k30*matL23+k40*matL24;
+  double t30 =  k00*matL30+k10*matL31+k20*matL32+k30*matL33+k40*matL34;
+  double t40 =  k00*matL40+k10*matL41+k20*matL42+k30*matL43+k40*matL44;
+  //
+  double t01 =  k01*matL00+k11*matL01+k21*matL02+k31*matL03+k41*matL04;
+  double t11 =  k01*matL10+k11*matL11+k21*matL12+k31*matL13+k41*matL14;
+  double t21 =  k01*matL20+k11*matL21+k21*matL22+k31*matL23+k41*matL24;
+  double t31 =  k01*matL30+k11*matL31+k21*matL32+k31*matL33+k41*matL34;
+  double t41 =  k01*matL40+k11*matL41+k21*matL42+k31*matL43+k41*matL44;
+  //
+  // (H^T * K^T * L * K * H) - (L * K * H) - (H^T * K^T * L) + (H^T*N^-1*H)
+  matL00 += k00*v00+k10*v10+k20*v20*k30*v30+k40*v40 - t00 - v00 + fCovIYZ[0];
+  matL01 += k01*v00+k11*v10+k21*v20*k31*v30+k41*v40 - t01 - v10 + fCovIYZ[1];
+  matL10 += k00*v01+k10*v11+k20*v21*k30*v31+k40*v41 - t10 - v01 + fCovIYZ[1];
+  matL11 += k01*v01+k11*v11+k21*v21*k31*v31+k41*v41 - t11 - v11 + fCovIYZ[2];
+  //
+  matL20 -= t20;
+  matL21 -= t21;
+  matL30 -= t30;
+  matL31 -= t31;
+  matL40 -= t40;
+  matL41 -= t41;
+  //
+  matL02 -= v20;
+  matL03 -= v30;
+  matL04 -= v40;
+  matL12 -= v21;
+  matL13 -= v31;
+  matL14 -= v41;
+  //
+  printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
+  printf("F: "); for (int i=0;i<kNFElem;i++) printf("%+e ",fFMatrix[i]); printf("\n");  
+  printf("K: "); for (int i=0;i<kNKElem;i++) printf("%+e ",fKMatrix[i]); printf("\n");  
+  //
+  for (int j=0;j<5;j++) {
+    for (int i=0;i<5;i++) printf("%+e ",matL[j][i]); printf("\n");  
+  }
+  //
   return kTRUE;
 }
+
+ */
index 5a9ee1d..c5b1450 100644 (file)
@@ -9,9 +9,11 @@ using namespace AliITSUAux;
 class AliITSUSeed: public AliExternalTrackParam
 {
  public:
-  enum {kKilled=BIT(14)};
+  enum {kKilled=BIT(14),kFake=BIT(15)};
   enum {kF02,kF04,kF12,kF13,kF14,kF24, kF44,kNFElem}; // non-trivial elems of propagation matrix
-  enum {kB00,kB01,kB10,kB11,kB20,kB21,kB30,kB31,kB40,kB41, kNBElem}; // non-trivial elems of B matrix (I - K*H)
+  enum {kK00,kK01,kK10,kK11,kK20,kK21,kK30,kK31,kK40,kK41, kNKElem}; // non-trivial elems of gain matrix
+  enum {kS00,kS10,kS11,kS20,kS21,kS22,kS30,kS31,kS32,kS33,kS40,kS41,kS42,kS43,kS44,kNSElem}; // elements of 5x5 sym matrix
+  enum {kR00,kR22,kNRElem};                        // non trivial elements of rotation matrix
   //
   AliITSUSeed();
   AliITSUSeed(const AliITSUSeed& src);
@@ -23,8 +25,9 @@ class AliITSUSeed: public AliExternalTrackParam
   void            SetLr(Int_t lr)                        {SetLrClusterID(lr,-1);} // lr w/o cluster
   void            SetLrClusterID(UInt_t id)              {fClID = id;}
   void            SetParent(TObject* par)                {fParent = par;}
-  void            SetChi2Cl(Double_t v)                  {fChi2Glo += fChi2Cl= v;}
+  void            SetChi2Cl(Double_t v)                  {fChi2Cl= v; v>0 ? fChi2Glo+=v : fChi2Penalty -= v;}
   void            Kill(Bool_t v=kTRUE)                   {SetBit(kKilled, v);}
+  void            SetFake(Bool_t v=kTRUE)                {SetBit(kFake, v);}
   //
   UInt_t          GetLrClusterID()                 const {return fClID;}
   Int_t           GetLrCluster(Int_t &lr)          const {return UnpackCluster(fClID,lr);}
@@ -35,10 +38,13 @@ class AliITSUSeed: public AliExternalTrackParam
   UShort_t        GetHitsPattern()                 const {return fHitsPattern;}
   Float_t         GetChi2Cl()                      const {return fChi2Cl;}
   Float_t         GetChi2Glo()                     const {return fChi2Glo;}
+  Float_t         GetChi2Penalty()                 const {return fChi2Penalty;}
   Float_t         GetChi2GloNrm()                  const;
   Bool_t          IsKilled()                       const {return TestBit(kKilled);}
+  Bool_t          IsFake()                         const {return TestBit(kFake);}
   //
   TObject*        GetParent()                      const {return fParent;}
+  const AliITSUSeed* GetParent(Int_t lr)              const;
   //
   virtual Bool_t  IsSortable()                     const {return kTRUE;}
   virtual Bool_t  IsEqual(const TObject* obj)      const;
@@ -49,20 +55,25 @@ class AliITSUSeed: public AliExternalTrackParam
   void            ApplyELoss2FMatrix(Double_t frac, Bool_t beforeProp);
   Bool_t          ApplyMaterialCorrection(Double_t xOverX0, Double_t xTimesRho, Double_t mass, Bool_t beforeProp);
   Bool_t          PropagateToX(Double_t xk, Double_t b);
+  Bool_t          RotateToAlpha(Double_t alpha);
   Bool_t          GetTrackingXAtXAlpha(double xOther,double alpOther,double bz, double &x);
   Double_t        GetPredictedChi2(Double_t p[2],Double_t cov[3]);
   Bool_t          Update();
+  Bool_t          Smooth(Double_t vecL[5],Double_t matL[15]);
+  Double_t*       ProdABA(const double a[15],const double b[15]) const;
   //
  protected:
   //
   Double_t              fFMatrix[kNFElem];  // matrix of propagation from prev layer (non-trivial elements)
+  Double_t              fKMatrix[kNKElem];  // Gain matrix non-trivial elements (note: standard MBF formula uses I-K*H)
+  Double_t              fRMatrix[kNRElem];  // rotation matrix non-trivial elements
   Double_t              fCovIYZ[3];         // inverted matrix of propagation + meas errors = [Hi * Pi|i-1 * Hi^T + Ri]^-1
   Double_t              fResid[2];          // residuals vector
-  Double_t              fBMatrix[kNBElem];  // K*H - I matix non-trivial elements (note: standard MBF formula uses I-K*H)
   UShort_t              fHitsPattern;       // bit pattern of hits
   UInt_t                fClID;              // packed cluster info (see AliITSUAux::PackCluster)
-  Float_t               fChi2Glo;           // current chi2 global
-  Float_t               fChi2Cl;            // track-cluster chi2
+  Float_t               fChi2Glo;           // current chi2 global (sum of track-cluster chi2's on layers with hit)
+  Float_t               fChi2Cl;            // track-cluster chi2 (if >0) or penalty for missing cluster (if < 0)
+  Float_t               fChi2Penalty;       // total penalty (e.g. for missing clusters)
   TObject*              fParent;            // parent track (in higher tree hierarchy)
   
   ClassDef(AliITSUSeed,1)
@@ -110,5 +121,17 @@ inline void AliITSUSeed::ApplyELoss2FMatrix(Double_t frac, Bool_t beforeProp)
   }
 }
 
+//_________________________________________________________________________
+inline const AliITSUSeed* AliITSUSeed::GetParent(Int_t lr) const
+{
+  // get parent at given layer
+  const AliITSUSeed* par=this;
+  int lrt;
+  while( par && (lrt=par->GetLayerID())>=0 ) {
+    if (lrt==lr) break;
+    par = dynamic_cast<const AliITSUSeed*>(par->GetParent());
+  }
+  return par;
+}
 
 #endif
diff --git a/ITS/UPGRADE/AliITSUTrackHyp.cxx b/ITS/UPGRADE/AliITSUTrackHyp.cxx
new file mode 100644 (file)
index 0000000..80bbe99
--- /dev/null
@@ -0,0 +1,62 @@
+#include "AliITSUTrackHyp.h"
+
+ClassImp(AliITSUTrackHyp)
+
+
+
+//__________________________________________________________________
+AliITSUTrackHyp::AliITSUTrackHyp(Int_t nlr) 
+: fNLayers(nlr)
+  ,fESDSeed(0)
+  ,fLayerSeeds(0)
+{
+  // def. c-tor
+  if (fNLayers>0) fLayerSeeds = new TObjArray[fNLayers];
+}
+
+//__________________________________________________________________
+AliITSUTrackHyp::~AliITSUTrackHyp() 
+{
+  // d-tor
+  delete[] fLayerSeeds;
+}
+
+//__________________________________________________________________
+AliITSUTrackHyp::AliITSUTrackHyp(const AliITSUTrackHyp &src)
+  : TObject(src)
+  , fNLayers(src.fNLayers)
+  , fESDSeed(src.fESDSeed)
+  , fLayerSeeds(0)
+{
+  // copy c-tor
+  if (fNLayers>0) {
+    fLayerSeeds = new TObjArray[fNLayers];
+    for (int ilr=fNLayers;ilr--;) {
+      int ns = src.GetNSeeds(ilr);
+      for (int isd=0;isd<ns;isd++) {
+       AliITSUSeed* sd = src.GetSeed(ilr,isd);
+       if (sd->IsKilled()) continue;
+       AddSeed(sd,ilr);
+      }      
+    }
+  }
+  //
+}
+
+//__________________________________________________________________
+AliITSUTrackHyp &AliITSUTrackHyp::operator=(const AliITSUTrackHyp &src)
+{
+  // copy 
+  if (this == &src) return *this;
+  this->~AliITSUTrackHyp();
+  new(this) AliITSUTrackHyp(src);
+  return *this;
+  //
+}
+
+//__________________________________________________________________
+void AliITSUTrackHyp::Print(Option_t* ) const
+{
+  printf("Track Hyp.#%4d. NSeeds:",GetUniqueID());
+  for (int i=0;i<fNLayers;i++) printf(" (%d) %3d",i,GetNSeeds(i)); printf("\n");
+}
diff --git a/ITS/UPGRADE/AliITSUTrackHyp.h b/ITS/UPGRADE/AliITSUTrackHyp.h
new file mode 100644 (file)
index 0000000..ffeb2d0
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef ALIITSUTRACKHYP_H
+#define ALIITSUTRACKHYP_H
+
+#include <TObject.h>
+#include <TObjArray.h>
+#include "AliITSUSeed.h"
+
+
+// Container for track hypotheses
+
+class AliITSUTrackHyp: public TObject
+{
+ public:
+  AliITSUTrackHyp(Int_t nlr=0);
+  AliITSUTrackHyp(const AliITSUTrackHyp& src);
+  AliITSUTrackHyp &operator=(const AliITSUTrackHyp &src);
+  virtual ~AliITSUTrackHyp();
+  //
+  Int_t              GetNLayers()        const {return fNLayers;}
+  Int_t              GetNSeeds(Int_t lr) const {return lr<fNLayers ? fLayerSeeds[lr].GetEntriesFast() : 1;}
+  AliITSUSeed*       GetSeed(Int_t lr, Int_t id) const;
+  AliITSUSeed*       GetESDSeed()        const {return fESDSeed;}
+  const TObjArray*   GetLayerSeeds(Int_t lr) const {return lr<fNLayers ? &fLayerSeeds[lr] : 0;}
+  void               AddSeed(AliITSUSeed* seed, Int_t lr);
+  void               SetESDSeed(AliITSUSeed* seed) {fESDSeed = seed;}
+  //
+  
+  //
+  virtual void       Print(Option_t* option = "") const;
+  //
+ protected:
+  UChar_t          fNLayers;               // number of layers
+  AliITSUSeed*     fESDSeed;               // bare esd (TPC) seed
+  TObjArray*       fLayerSeeds;            // seeds of given layer
+  //
+  ClassDef(AliITSUTrackHyp,1)
+};
+
+//___________________________________________________________________
+inline AliITSUSeed* AliITSUTrackHyp::GetSeed(Int_t lr, Int_t id) const
+{
+  //return requested seed of given layer, no check is done on seed index
+  return lr<fNLayers ? (AliITSUSeed*)fLayerSeeds[lr].UncheckedAt(id) : fESDSeed;
+}
+
+//___________________________________________________________________
+inline void AliITSUTrackHyp::AddSeed(AliITSUSeed* seed, Int_t lr)
+{
+  //return requested seed of given layer, no check is done
+  if (lr<fNLayers) fLayerSeeds[lr].AddLast(seed);
+  else             fESDSeed = seed;
+}
+
+#endif
index ed51e0d..7f83805 100644 (file)
@@ -41,8 +41,6 @@ using namespace TMath;
 
 
 //----------------- tmp stuff -----------------
-Int_t currLabel=-1;
-
 
 ClassImp(AliITSUTrackerGlo)
 //_________________________________________________________________________
@@ -51,7 +49,8 @@ AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
   ,fITS(0)
   ,fCurrESDtrack(0)
   ,fCurrMass(kPionMass)
-  ,fSeedsLr(0)
+  ,fHypStore(100)
+  ,fCurrHyp(0)
   ,fSeedsPool("AliITSUSeed",0)
   ,fTrCond()
 {
@@ -65,7 +64,6 @@ AliITSUTrackerGlo::~AliITSUTrackerGlo()
  // Default destructor
  //  
   delete fITS;
-  delete[] fSeedsLr;
   //
 }
 
@@ -80,8 +78,6 @@ void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
   }
   //
   fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
-  int n = fITS->GetNLayersActive()+1;
-  fSeedsLr = new TObjArray[n];
   //
   fTrCond.SetNLayers(fITS->GetNLayersActive());
   fTrCond.AddNewCondition(5);
@@ -109,17 +105,21 @@ Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
   //
   //
   AliITSUReconstructor::GetRecoParam()->Print();
-
+  int nTrESD = esdEv->GetNumberOfTracks();
+  fHypStore.Delete();
+  if (fHypStore.GetSize()<nTrESD) fHypStore.Expand(nTrESD+100);
+  //
   fITS->ProcessClusters();
   // select ESD tracks to propagate
-  int nTrESD = esdEv->GetNumberOfTracks();
   for (int itr=0;itr<nTrESD;itr++) {
     AliESDtrack *esdTr = esdEv->GetTrack(itr);
     AliInfo(Form("Processing track %d | MCLabel: %d",itr,esdTr->GetTPCLabel()));
-    currLabel = Abs(esdTr->GetTPCLabel());
-    FindTrack(esdTr);
+    FindTrack(esdTr, itr);
   }
-
+  //
+  printf("Hypotheses for current event (N seeds in pool: %d, size: %d)\n",fSeedsPool.GetEntriesFast(),fSeedsPool.GetSize());
+  fHypStore.Print();
+  //
   return 0;
 }
 
@@ -197,12 +197,12 @@ Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
 }
 
 //_________________________________________________________________________
-void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr)
+void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
 {
   // find prolongaion candidates finding for single seed
   //
   if (!NeedToProlong(esdTr)) return;  // are we interested in this track?
-  if (!InitSeed(esdTr))      return;  // initialize prolongations hypotheses tree
+  if (!InitSeed(esdTr,esdID)) return;  // initialize prolongations hypotheses tree
   //
   AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
   AliITSUSeed seedUC;  // copy of the seed from the upper layer
@@ -212,9 +212,9 @@ void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr)
   //
   for (int ila=fITS->GetNLayersActive();ila--;) {
     int ilaUp = ila+1;                         // prolong seeds from layer above
-    int nSeedsUp = GetNSeeds(ilaUp);
+    int nSeedsUp = fCurrHyp->GetNSeeds(ilaUp);
     for (int isd=0;isd<nSeedsUp;isd++) {
-      AliITSUSeed* seedU = GetSeed(ilaUp,isd);  // seed on prev.active layer to prolong
+      AliITSUSeed* seedU = fCurrHyp->GetSeed(ilaUp,isd);  // seed on prev.active layer to prolong
       seedUC = *seedU;                          // its copy will be prolonged
       seedUC.SetParent(seedU);
       seedUC.ResetFMatrix();                    // reset the matrix for propagation to next layer
@@ -245,7 +245,8 @@ void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr)
        double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
        if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) continue;
        if (!seedT.PropagateToX(xs,bz)) continue;
-       if (!seedT.Rotate(sens->GetPhiTF())) continue;
+       //      if (!seedT.Rotate(sens->GetPhiTF())) continue;
+       if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
        //
        int clID0 = sens->GetFirstClusterId();
        for (int icl=sens->GetNClusters();icl--;) {
@@ -258,43 +259,63 @@ void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr)
       }
       // cluster search is done. Do we need ta have a version of this seed skipping current layer
       seedT.SetLr(ila);
-      if (!NeedToKill(&seedT,kMissingCluster)) AddProlongationHypothesis(NewSeedFromPool(&seedT) ,ila);      
+      if (!NeedToKill(&seedT,kMissingCluster)) {
+       AliITSUSeed* seedSkp = NewSeedFromPool(&seedT);
+       double penalty = -AliITSUReconstructor::GetRecoParam()->GetMissPenalty(ila);
+       // to do: make penalty to account for probability to miss the cluster for good reason
+       seedSkp->SetChi2Cl(penalty);
+       AddProlongationHypothesis(seedSkp,ila);      
+      }
     }
+    ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
     printf(">>> All hypotheses on lr %d: \n",ila);
-    for (int ih=0;ih<GetNSeeds(ila);ih++) {
-      printf(" #%3d ",ih); GetSeed(ila,ih)->Print();
+    for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
+      printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();
+      //
+      /*
+      if (ila!=0) continue;
+      double vecL[5] = {0};
+      double matL[15] = {0};
+      AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
+      while(sp->GetParent()) {
+       sp->Smooth(vecL,matL);
+       if (sp->GetLayerID()>=fITS->GetNLayersActive()-1) break;
+       sp = (AliITSUSeed*)sp->GetParent();
+      }
+      */
     }
   }
   //
-  ResetSeedTree();
+  SaveCurrentTrackHypotheses();
+  //
 }
 
 //_________________________________________________________________________
-Bool_t AliITSUTrackerGlo::InitSeed(AliESDtrack *esdTr)
+Bool_t AliITSUTrackerGlo::InitSeed(AliESDtrack *esdTr, Int_t esdID)
 {
   // init prolongaion candidates finding for single seed
+  fCurrHyp = GetTrackHyp(esdID);
+  if (fCurrHyp) return kTRUE;
+  //
   fCurrMass = esdTr->GetMass();
   fCurrESDtrack = esdTr;
   if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
   //
-  //
   AliITSUSeed* seed = NewSeedFromPool();
   seed->SetLr(fITS->GetNLayersActive());   // fake layer
   seed->AliExternalTrackParam::operator=(*esdTr);
   seed->SetParent(esdTr);
-  AddProlongationHypothesis(seed,fITS->GetNLayersActive());
+  int nl = fITS->GetNLayersActive();
+  fCurrHyp = new AliITSUTrackHyp(nl);
+  fCurrHyp->SetESDSeed(seed);
+  AddProlongationHypothesis(seed,nl);
+  fCurrHyp->SetUniqueID(esdID);
+  SetTrackHyp(fCurrHyp,esdID);
   return kTRUE;
   // TO DO
 }
 
 //_________________________________________________________________________
-void AliITSUTrackerGlo::ResetSeedTree()
-{
-  // reset current hypotheses tree
-  for (int i=fITS->GetNLayersActive()+1;i--;) fSeedsLr[i].Clear();
-}
-
-//_________________________________________________________________________
 Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
 {
   // transport seed from layerFrom to the entrance of layerTo
@@ -397,7 +418,10 @@ Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
   AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
   //
   Bool_t goodCl = kFALSE;
-  for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) goodCl = kTRUE;
+  int currLabel = Abs(fCurrESDtrack->GetTPCLabel());
+  //
+  if (cl->GetLabel(0)>=0) {for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}}
+  else goodCl = kTRUE;
   //
   if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
     if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
@@ -430,7 +454,8 @@ Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
   double chi2 = track->GetPredictedChi2(p,cov);
   if (chi2>AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr)) {
     if (goodCl) {
-      printf("Loose good cl: Chi2=%e > Chi2Max=%e \n",chi2,AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr)); 
+      printf("Loose good cl: Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e\n",
+            chi2,AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr),dy,dz); 
       track->Print("etp");
       cl->Print("");
     }
@@ -446,7 +471,10 @@ Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
   track->SetLrClusterID(lr,clID);
   cl->IncreaseClusterUsage();
   //
-  AliInfo(Form("AddCl Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f",clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2));
+  track->SetFake(!goodCl);
+  //
+  AliInfo(Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)",
+              goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
   //
   AddProlongationHypothesis(track,lr);
   //
@@ -497,3 +525,12 @@ Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Doubl
   }
   return kTRUE;
 }
+
+//______________________________________________________________________________
+void AliITSUTrackerGlo::SaveCurrentTrackHypotheses()
+{
+  // RS: shall we clean up killed seeds?
+  fCurrHyp = 0;
+  // TODO
+  
+}
index 45995a1..8b87936 100644 (file)
@@ -11,6 +11,7 @@
 #include "AliESDEvent.h"
 #include "AliITSUSeed.h"
 #include "AliITSUTrackCond.h"
+#include "AliITSUTrackHyp.h"
 
 class AliITSUReconstructor;
 class AliITSURecoDet;
@@ -51,10 +52,8 @@ class AliITSUTrackerGlo : public AliTracker {
   //------------------------------------
   Bool_t                 NeedToProlong(AliESDtrack* estTr);
   void                   Init(AliITSUReconstructor* rec);
-  void                   FindTrack(AliESDtrack* esdTr);
-  Bool_t                 InitSeed(AliESDtrack *esdTr);
-  Int_t                  GetNSeeds(Int_t lr)              const {return fSeedsLr[lr].GetEntriesFast();} //RS TOCHECK
-  AliITSUSeed*           GetSeed(Int_t lr, Int_t sID)     const {return (AliITSUSeed*)fSeedsLr[lr].UncheckedAt(sID);} //RS TOCHECK
+  void                   FindTrack(AliESDtrack* esdTr, Int_t esdID);
+  Bool_t                 InitSeed(AliESDtrack *esdTr, Int_t esdID);
   Bool_t                 TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo);
   Bool_t                 PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep=1.0, Bool_t matCorr=kTRUE);
   //
@@ -64,9 +63,11 @@ class AliITSUTrackerGlo : public AliTracker {
   void                   AddProlongationHypothesis(AliITSUSeed* seed, Int_t lr);
   //
   AliITSUSeed*           NewSeedFromPool(const AliITSUSeed* src=0);
-  void                   DeleteLastSeedFromPool()               {fSeedsPool.RemoveLast();}
-  void                   ResetSeedTree();  // RS TOCHECK
-  //
+  AliITSUTrackHyp*       GetTrackHyp(Int_t id)               const  {return (AliITSUTrackHyp*)fHypStore.UncheckedAt(id);}
+  void                   SetTrackHyp(AliITSUTrackHyp* hyp,Int_t id) {fHypStore.AddAtAndExpand(hyp,id);}
+  void                   DeleteLastSeedFromPool()                   {fSeedsPool.RemoveLast();}
+  void                   SaveCurrentTrackHypotheses();
+ //
 
  private:
   
@@ -81,7 +82,8 @@ class AliITSUTrackerGlo : public AliTracker {
   Double_t                        fTrImpData[kNTrImpData];  // data on track impact on the layer
   //
   // the seeds management to be optimized
-  TObjArray*                      fSeedsLr;        // seeds at each layer
+  TObjArray                       fHypStore;       // storage for tracks hypotheses
+  AliITSUTrackHyp*                fCurrHyp;        // hypotheses container for current track
   TClonesArray                    fSeedsPool;      // pool for seeds
   //
   AliITSUTrackCond                fTrCond;         // tmp, to be moved to recoparam
@@ -94,7 +96,7 @@ class AliITSUTrackerGlo : public AliTracker {
 inline void AliITSUTrackerGlo::AddProlongationHypothesis(AliITSUSeed* seed, Int_t lr)
 {
   // add new seed prolongation hypothesis 
-  fSeedsLr[lr].AddLast(seed);
+  fCurrHyp->AddSeed(seed,lr);
   printf("*** Adding: "); seed->Print();
 }
 
index 3d5ca2c..6cd57d5 100644 (file)
@@ -37,6 +37,8 @@ set ( SRCS
     AliITSUSeed.cxx
     AliITSUTrackerGlo.cxx
     AliITSUTrackCond.cxx
+    AliITSUTrackHyp.cxx
+
 #
 #    v0/AliITSlayerUpgrade.cxx 
 #    v0/AliITStrackerUpgrade.cxx 
index 78a3c93..bcfd7df 100644 (file)
@@ -23,6 +23,8 @@
 #pragma link C++ class AliITSUSeed+;
 #pragma link C++ class AliITSUTrackerGlo+;
 #pragma link C++ class AliITSUTrackCond+;
+#pragma link C++ class AliITSUTrackHyp+;
+
 //
 
 // old v0