modifications for PropagateBack
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 Jan 2013 19:52:57 +0000 (19:52 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 Jan 2013 19:52:57 +0000 (19:52 +0000)
ITS/UPGRADE/AliITSURecoDet.cxx
ITS/UPGRADE/AliITSURecoDet.h
ITS/UPGRADE/AliITSURecoLayer.cxx
ITS/UPGRADE/AliITSURecoLayer.h
ITS/UPGRADE/AliITSUTrackHyp.cxx
ITS/UPGRADE/AliITSUTrackHyp.h
ITS/UPGRADE/AliITSUTrackerGlo.cxx
ITS/UPGRADE/AliITSUTrackerGlo.h

index 5e4dd90..8149b68 100644 (file)
@@ -134,3 +134,53 @@ void AliITSURecoDet::IndexLayers()
     SetRMax(GetLayer(fNLayers-1)->GetRMax()+kRMargin);
   }
 }
+
+//______________________________________________________
+Int_t AliITSURecoDet::FindLastLayerID(Double_t r, int dir) const
+{
+  // find the last layer which the particle moving in direction dir (1:outward,-1:inward) 
+  // will traverse on its way to radius r 
+  int ilr;
+  //
+  if (dir>0) {
+    for (ilr=0;ilr<fNLayers;ilr++) {
+      AliITSURecoLayer* lr = GetLayer(ilr);
+      if ( r<lr->GetR(-dir) ) break;  // this layer at least entered
+    }
+    return ilr--;  // -1 will correspond to point below the smalles layer
+  }
+  else {
+    for (ilr=fNLayers;ilr--;) {
+      AliITSURecoLayer* lr = GetLayer(ilr);
+      if ( r>lr->GetR(-dir) ) break; // this layer at least entered
+    }
+    ilr++;
+    return ilr<fNLayers ? ilr:-1; // -1 will correspond to point above outer layer
+  }
+  //
+}
+
+//______________________________________________________
+Int_t AliITSURecoDet::FindFirstLayerID(Double_t r, int dir) const
+{
+  // find the first layer which the particle moving in direction dir (1:outward,-1:inward) 
+  // will traverse starting from radius r 
+  int ilr;
+  //
+  if (dir>0) {
+    for (ilr=0;ilr<fNLayers;ilr++) {
+      AliITSURecoLayer* lr = GetLayer(ilr);
+      if ( r<lr->GetR(dir) ) break;  // this layer at least entered
+    }
+    return ilr<fNLayers ? ilr:-1;  // -1 will correspond to point above outer leayer
+  }
+  else {
+    for (ilr=fNLayers;ilr--;) {
+      AliITSURecoLayer* lr = GetLayer(ilr);
+      if ( r>lr->GetR(dir) ) break; // this layer at least entered
+      break;
+    }
+    return ilr; // -1 will correspond to point below inner layer
+  }
+  //
+}
index bbbcb4a..85cb552 100644 (file)
@@ -26,6 +26,8 @@ class AliITSURecoDet : public TNamed
   Int_t              GetNLayers()                  const {return fNLayers;}
   Int_t              GetNLayersActive()            const {return fNLayersActive;}
   Int_t              GetLrIDActive(Int_t lrActID)  const;
+  Int_t              FindLastLayerID(Double_t r, int dir)  const;
+  Int_t              FindFirstLayerID(Double_t r, int dir) const;
   AliITSURecoLayer*  GetLayer(Int_t i)             const;
   AliITSURecoLayer*  GetLayerActive(Int_t i)       const;
   AliITSUGeomTGeo*   GetGeom()                     const {return fGeom;}
@@ -75,14 +77,14 @@ inline Int_t AliITSURecoDet::GetLrIDActive(Int_t lrActID) const
 inline AliITSURecoLayer* AliITSURecoDet::GetLayer(Int_t i) const 
 {
   // get layer with global id=i
-  return i<fNLayers ? (AliITSURecoLayer*)fLayers.UncheckedAt(i):0;
+  return i>=0&&i<fNLayers ? (AliITSURecoLayer*)fLayers.UncheckedAt(i):0;
 }
 
 //_____________________________________________________________
 inline AliITSURecoLayer* AliITSURecoDet::GetLayerActive(Int_t i) const
 {
   // get layer with activeID=i
-  return i<fNLayersActive ? (AliITSURecoLayer*)fLayersActive.UncheckedAt(i):0;
+  return i>=0&&i<fNLayersActive ? (AliITSURecoLayer*)fLayersActive.UncheckedAt(i):0;
 }
 
 //______________________________________________________
index 02c6919..17cc909 100644 (file)
@@ -1,9 +1,10 @@
 #include <TClonesArray.h>
 #include "AliITSURecoLayer.h"
-#include "AliITSUGeomTGeo.h"
 #include "AliITSsegmentation.h"
 #include "AliITSUAux.h"
 #include "AliITSUClusterPix.h"
+#include "AliITSUGeomTGeo.h"
+#include "AliLog.h"
 
 using namespace AliITSUAux;
 using namespace TMath;
@@ -312,3 +313,12 @@ Int_t  AliITSURecoLayer::Compare(const TObject* obj) const
   return dr>0 ? 1:-1;
   //      
 }
+
+//_________________________________________________________________
+AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const 
+{
+  // get sensor from its global id
+  i -= fITSGeom->GetFirstModIndex(fActiveID);
+  if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstModIndex(fActiveID),fActiveID));
+  return (AliITSURecoSens*)fSensors[i];
+}
index 343d2b3..677ea6f 100644 (file)
@@ -5,6 +5,7 @@
 #include <TObjArray.h>
 #include <TClonesArray.h>
 #include "AliITSURecoSens.h"
+
 class AliITSUGeomTGeo;
 class AliITSsegmentation;
 class AliCluster;
@@ -57,6 +58,7 @@ class AliITSURecoLayer : public TNamed
   //
   AliITSURecoSens*   GetSensor(Int_t i)            const {return i<0 ? 0:(AliITSURecoSens*)fSensors[i];}
   AliITSURecoSens*   GetSensor(Int_t ld,Int_t is)  const {return GetSensor(ld*fNSensInLadder+is);}
+  AliITSURecoSens*   GetSensorFromID(Int_t i)      const;
   TClonesArray*      GetClusters()                 const {return (TClonesArray*)fClusters;}
   Int_t              GetNClusters()                const {return fClusters ? fClusters->GetEntriesFast() : 0;}
   AliCluster*        GetCluster(Int_t icl)         const {return (AliCluster*)fClusters->UncheckedAt(icl);}
index 28ac92d..b8d3de5 100644 (file)
@@ -81,11 +81,12 @@ void AliITSUTrackHyp::DefineWinner(int lr, int id)
 }
 
 //__________________________________________________________________
-Double_t AliITSUTrackHyp::GetPredictedChi2(const AliCluster */*c*/) const
+Double_t AliITSUTrackHyp::GetPredictedChi2(const AliCluster *cl) const
 {
-  // NA
-  AliFatal("Not to be used");
-  return 0;
+  // calculate chi2 to cluster
+  Double_t p[2]={cl->GetY(), cl->GetZ()};
+  Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
+  return AliExternalTrackParam::GetPredictedChi2(p,cov);
 }
 
 //__________________________________________________________________
@@ -103,3 +104,34 @@ Bool_t AliITSUTrackHyp::Update(const AliCluster* /*c*/, Double_t /*chi2*/, Int_t
   AliFatal("Not to be used");
   return kFALSE;
 }
+
+//__________________________________________________________________
+Bool_t AliITSUTrackHyp::Update(const AliCluster* cl)
+{
+  // update with cluster
+  Double_t p[2]={cl->GetY(), cl->GetZ()};
+  Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
+  double chi2 = AliExternalTrackParam::GetPredictedChi2(p,cov);
+  if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
+  SetChi2(GetChi2()+chi2);
+  return kTRUE;
+}
+
+//__________________________________________________________________
+void AliITSUTrackHyp::FetchClusterInfo(Int_t *clIDarr) const
+{
+  // fill cl.id's in the array. The clusters of layer L will be set at slots
+  // clID[2L] (and clID[2L+1] if there is an extra cluster).
+  for (int i=fNLayers<<1;i--;) clIDarr[i]=-1;
+  AliITSUSeed* seed = GetWinner();
+  Int_t lr;
+  while(seed) {
+    int clID = seed->GetLrCluster(lr);
+    if (clID>=0) {
+      int slotLr = lr<<1;
+      clIDarr[ clIDarr[slotLr]<0 ? slotLr : slotLr+1 ] = clID;
+    }
+    seed = (AliITSUSeed*)seed->GetParent();
+  }
+}
+
index bf9144a..7a3b2dd 100644 (file)
@@ -27,8 +27,10 @@ class AliITSUTrackHyp: public AliKalmanTrack
   const TObjArray*   GetLayerSeeds(Int_t lr) const {return lr<fNLayers ? &fLayerSeeds[lr] : 0;}
   void               AddSeed(AliITSUSeed* seed, Int_t lr) {fLayerSeeds[lr].AddLast(seed);}
   void               SetESDTrack(AliESDtrack* esdtr) {fESDTrack = esdtr;}
+  void               FetchClusterInfo(Int_t* clIDarr) const;
   //
   void               SetChi2(Double_t chi2) {fChi2 = chi2;}
+  Bool_t             Update(const AliCluster* c);
   //
   virtual Double_t   GetPredictedChi2(const AliCluster *c) const;
   virtual Bool_t     PropagateTo(Double_t xr, Double_t x0, Double_t rho);
index 8020cfd..fa0621b 100644 (file)
@@ -54,6 +54,7 @@ AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
   ,fSeedsPool("AliITSUSeed",0)
   ,fTrCond()
   ,fTrackPhase(-1)
+  ,fClInfo(0)
 {
   // Default constructor
   if (rec) Init(rec);
@@ -65,6 +66,7 @@ AliITSUTrackerGlo::~AliITSUTrackerGlo()
  // Default destructor
  //  
   delete fITS;
+  delete[] fClInfo;
   //
 }
 
@@ -74,7 +76,9 @@ void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
   // init with external reconstructor
   //
   fITS = new AliITSURecoDet(rec->GetGeom(),"ITSURecoInterface");
-  for (int ilr=fITS->GetNLayersActive();ilr--;) {
+  int nLr = fITS->GetNLayersActive();
+  fClInfo = new Int_t[nLr<<1];
+  for (int ilr=nLr;ilr--;) {
     fITS->GetLayerActive(ilr)->SetClusters(rec->GetClusters(ilr));
   }
   //
@@ -313,7 +317,8 @@ void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
        // since the transport matrix should be defined in this frame.
        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 (!PropagateSeed(&seedT,xs,fCurrMass)) continue;
+       //      if (!seedT.PropagateToX(xs,bz)) continue;
        //      if (!seedT.Rotate(sens->GetPhiTF())) continue;
        if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
        //
@@ -387,6 +392,9 @@ Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t
   //  
   const double kToler = 1e-6; // tolerance for layer on-surface check
   //
+  //
+  if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
+  //
   int dir = lTo > lFrom ? 1:-1;
   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
   Bool_t checkFirst = kTRUE;
@@ -394,7 +402,7 @@ Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t
     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
     if (lrFr) {
       Bool_t doLayer = kTRUE;
-      double xToGo = dir>0 ? lrFr->GetRMax() : lrFr->GetRMin();
+      double xToGo = lrFr->GetR(dir);
       if (checkFirst) { // do we need to track till the surface of the current layer ?
        checkFirst = kFALSE;
        if      (dir>0) { if (curR2-xToGo*xToGo>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
@@ -410,7 +418,7 @@ Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t
     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
     //
     // go the entrance of the layer, assuming no materials in between
-    double xToGo = dir>0 ? lrTo->GetRMin() : lrTo->GetRMax();
+    double xToGo = lrTo->GetR(-dir);
     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) return kFALSE;
     lrFr = lrTo;
@@ -426,6 +434,8 @@ Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lF
   //  
   const double kToler = 1e-6; // tolerance for layer on-surface check
   //
+  if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
+  //
   int dir = lTo > lFrom ? 1:-1;
   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
   Bool_t checkFirst = kTRUE;
@@ -433,7 +443,7 @@ Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lF
     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
     if (lrFr) {
       Bool_t doLayer = kTRUE;
-      double xToGo = dir>0 ? lrFr->GetRMax() : lrFr->GetRMin();
+      double xToGo = lrFr->GetR(dir);
       if (checkFirst) { // do we need to track till the surface of the current layer ?
        checkFirst = kFALSE;
        if      (dir>0) { if (curR2-xToGo*xToGo>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
@@ -449,7 +459,7 @@ Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lF
     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
     //
     // go the entrance of the layer, assuming no materials in between
-    double xToGo = dir>0 ? lrTo->GetRMin() : lrTo->GetRMax();
+    double xToGo = lrTo->GetR(-dir);
     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) return kFALSE;
     lrFr = lrTo;
@@ -524,7 +534,9 @@ Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
   Bool_t goodCl = kFALSE;
   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;}}
+  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
@@ -609,22 +621,34 @@ Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Doubl
   Int_t dir         = (xpos<xToGo) ? 1:-1;
   Double_t xyz0[3],xyz1[3],param[7];
   //
+  Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
   if (matCorr) seed->GetXYZ(xyz1);   //starting global position
   while ( (xToGo-xpos)*dir > kEpsilon){
     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
     Double_t x    = xpos+step;
     Double_t bz=GetBz();   // getting the local Bz
     if (!seed->PropagateToX(x,bz))  return kFALSE;
-    if (matCorr) {
+    double ds = 0;
+    if (matCorr || updTime) {
       xyz0[0]=xyz1[0]; // global pos at the beginning of step
       xyz0[1]=xyz1[1];
       xyz0[2]=xyz1[2];
       seed->GetXYZ(xyz1);    //  // global pos at the end of step
-      MeanMaterialBudget(xyz0,xyz1,param);     
-      Double_t xrho=param[0]*param[4], xx0=param[1];
-      if (dir>0) xrho = -xrho; // outward should be negative
-      if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) return kFALSE;
+      if (matCorr) {
+       MeanMaterialBudget(xyz0,xyz1,param);    
+       Double_t xrho=param[0]*param[4], xx0=param[1];
+       if (dir>0) xrho = -xrho; // outward should be negative
+       if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) return kFALSE;
+       ds = param[4];
+      }
+       else { // matCorr is not requested but time integral is
+       double d0 = xyz1[0]-xyz0[0];
+       double d1 = xyz1[1]-xyz0[1];
+       double d2 = xyz1[2]-xyz0[2];    
+       ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
+      }     
     }
+    if (updTime) seed->AddTimeStep(ds);
     xpos = seed->GetX();
   }
   return kTRUE;
@@ -717,24 +741,60 @@ void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp)
 Bool_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest)
 {
   // refit track till radius rDest
+  AliITSUTrackHyp tmpTr;
+  //
   double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
   int dir,lrStart,lrStop;
   //
-  if (rCurr<rDest) {
-    dir = 1;
-    lrStart = 0;
-    lrStop  = fITS->GetNLayers()-1;
-  }
-  else {
-    dir = 1;
-    lrStart = fITS->GetNLayers()-1;
-    lrStop  = 0;    
-  }
+  dir = rCurr<rDest ? 1 : -1;
+  lrStart = fITS->FindFirstLayerID(rCurr,dir);
+  lrStop  = fITS->FindLastLayerID(rDest,dir);
+  if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
+  //
+  trc->FetchClusterInfo(fClInfo);
+  fCurrMass = trc->GetMass();
+  tmpTr.AliKalmanTrack::operator=(*trc);
+  tmpTr.SetChi2(0);
+  int iclLr[2],nclLr;
   //
   for (int ilr=lrStart;ilr!=lrStop;ilr+=dir) {
     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
     if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
+    int ilrA2,ilrA = lr->GetActiveID();
+    // passive layer or active w/o hits will be traversed on the way to next cluster
+    if (!lr->IsActive() || fClInfo[ilrA2=(ilrA<<1)]<0) continue; 
+    //
+    if (!TransportToLayer(&tmpTr,lrStart,ilr)) return kFALSE; 
+    lrStart = ilr;
+    //
+    // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
+    nclLr=0;
+    if (dir>0) { // clusters are stored in increasing radius order
+      iclLr[nclLr++]=fClInfo[ilrA2++];
+      if (fClInfo[ilrA2]>=0) iclLr[nclLr++]=fClInfo[ilrA2];
+    }
+    else {
+      if ( fClInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=fClInfo[ilrA2+1];
+      iclLr[nclLr++]=fClInfo[ilrA2];
+    }
+    //
+    for (int icl=0;icl<nclLr;icl++) {
+      AliITSUClusterPix* clus =  (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
+      AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
+      if (!tmpTr.Rotate(sens->GetPhiTF())) return kFALSE;
+      if (!PropagateSeed(&tmpTr,sens->GetXTF()+clus->GetX(),fCurrMass)) return kFALSE;
+      if (!tmpTr.Update(clus)) return kFALSE;
+    }
     //
   }
-  
+  // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
+  // Still, try to go as close as possible to rDest.
+  //
+  if (lrStart!=lrStop) {
+    if (!TransportToLayer(&tmpTr,lrStart,lrStop)) return kTRUE;
+    // go to the exit from layer
+    // TODO
+  }
+  trc->AliKalmanTrack::operator=(tmpTr);
+  return kTRUE;
 }
index 211cc07..d1931a3 100644 (file)
@@ -94,6 +94,7 @@ class AliITSUTrackerGlo : public AliTracker {
   //
   AliITSUTrackCond                fTrCond;         // tmp, to be moved to recoparam
   Int_t                           fTrackPhase;     // tracking phase
+  Int_t*                          fClInfo;         //! auxiliary track cluster info
   //
   ClassDef(AliITSUTrackerGlo,1)   //ITS upgrade tracker