ITS/UPGRADE/AliITSUTrackerGlo.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17   TO CHECK
  18   GetBz usage
  19
  20  */
  21
  22 //-------------------------------------------------------------------------
  23 //  Implementation of the ITS Upgrade tracker mother class.
  24 //-------------------------------------------------------------------------
  25 #include <TTree.h>
  26 #include <Riostream.h>
  27 #include <TMath.h>
  28
  29 #include "AliITSUTrackerGlo.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDtrack.h"
  32 #include "AliITSURecoDet.h"
  33 #include "AliITSURecoSens.h"
  34 #include "AliITSUReconstructor.h"
  35 #include "AliITSReconstructor.h"
  36 #include "AliITSUSeed.h"
  37 #include "AliITSUAux.h"
  38 #include "AliITSUClusterPix.h"
  39 using namespace AliITSUAux;
  40 using namespace TMath;
  41
  42
  43 //----------------- tmp stuff -----------------
  44
  45 ClassImp(AliITSUTrackerGlo)
  46 //_________________________________________________________________________
  47 AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
  48 :  fReconstructor(rec)
  49   ,fITS(0)
  50   ,fCurrESDtrack(0)
  51   ,fCurrMass(kPionMass)
  52   ,fHypStore(100)
  53   ,fCurrHyp(0)
  54   ,fSeedsPool("AliITSUSeed",0)
  55   ,fTrCond()
  56   ,fTrackPhase(-1)
  57 {
  58   // Default constructor
  59   if (rec) Init(rec);
  60 }
  61
  62 //_________________________________________________________________________
  63 AliITSUTrackerGlo::~AliITSUTrackerGlo()
  64 {
  65  // Default destructor
  66  //
  67   delete fITS;
  68   //
  69 }
  70
  71 //_________________________________________________________________________
  72 void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
  73 {
  74   // init with external reconstructor
  75   //
  76   fITS = new AliITSURecoDet(rec->GetGeom(),"ITSURecoInterface");
  77   for (int ilr=fITS->GetNLayersActive();ilr--;) {
  78     fITS->GetLayerActive(ilr)->SetClusters(rec->GetClusters(ilr));
  79   }
  80   //
  81   fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
  82   //
  83   fTrCond.SetNLayers(fITS->GetNLayersActive());
  84   fTrCond.AddNewCondition(5);
  85   fTrCond.AddGroupPattern( (0x1<<0)|(0x1<<1) );
  86   fTrCond.AddGroupPattern( (0x1<<3)|(0x1<<4) );
  87   fTrCond.AddGroupPattern( (0x1<<5)|(0x1<<6) );
  88   //
  89   fTrCond.AddNewCondition(5);
  90   fTrCond.AddGroupPattern( (0x1<<0)|(0x1<<2) );
  91   fTrCond.AddGroupPattern( (0x1<<3)|(0x1<<4) );
  92   fTrCond.AddGroupPattern( (0x1<<5)|(0x1<<6) );
  93   //
  94   fTrCond.AddNewCondition(5);
  95   fTrCond.AddGroupPattern( (0x1<<1)|(0x1<<2) );
  96   fTrCond.AddGroupPattern( (0x1<<3)|(0x1<<4) );
  97   fTrCond.AddGroupPattern( (0x1<<5)|(0x1<<6) );
  98   //
  99   printf("Tracking Conditions: ");
 100   fTrCond.Print();
 101 }
 102
 103 //_________________________________________________________________________
 104 Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
 105 {
 106   //
 107   //
 108   fTrackPhase = kClus2Tracks;
 109   int nTrESD = esdEv->GetNumberOfTracks();
 110   AliInfo(Form("Will try to find prolongations for %d tracks",nTrESD));
 111   AliITSUReconstructor::GetRecoParam()->Print();
 112   fHypStore.Delete();
 113   if (fHypStore.GetSize()<nTrESD) fHypStore.Expand(nTrESD+100);
 114   //
 115   fITS->ProcessClusters();
 116   // select ESD tracks to propagate
 117   for (int itr=0;itr<nTrESD;itr++) {
 118     AliESDtrack *esdTr = esdEv->GetTrack(itr);
 119     AliInfo(Form("Processing track %d | MCLabel: %d",itr,esdTr->GetTPCLabel()));
 120     if (!NeedToProlong(esdTr)) continue;  // are we interested in this track?
 121     FindTrack(esdTr, itr);
 122   }
 123   //
 124   printf("Hypotheses for current event (N seeds in pool: %d, size: %d)\n",fSeedsPool.GetEntriesFast(),fSeedsPool.GetSize());
 125   fHypStore.Print();
 126   FinalizeHypotheses();
 127   //
 128   return 0;
 129 }
 130
 131 //_________________________________________________________________________
 132 Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent *esdEv)
 133 {
 134   //
 135   // Do outward fits in ITS
 136   //
 137   fTrackPhase = kPropBack;
 138   int nTrESD = esdEv->GetNumberOfTracks();
 139   AliInfo(Form("Will fit %d tracks",nTrESD));
 140   //
 141   double bz0 = GetBz();
 142   Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
 143   AliITSUTrackHyp dummyTr,*currTr=0;
 144   const double kWatchStep=10.; // for larger steps watch arc vs segment difference
 145   Double_t times[AliPID::kSPECIES];
 146   //
 147   //
 148   for (int itr=0;itr<nTrESD;itr++) {
 149     AliESDtrack *esdTr = esdEv->GetTrack(itr);
 150     // Start time integral and add distance from current position to vertex
 151     if (esdTr->IsOn(AliESDtrack::kITSout)) continue; // already done
 152     //
 153     esdTr->GetXYZ(xyzTrk);
 154     Double_t dst = 0.;     // set initial track lenght, tof
 155     {
 156       double dxs = xyzTrk[0] - xyzVtx[0];
 157       double dys = xyzTrk[1] - xyzVtx[1];
 158       double dzs = xyzTrk[2] - xyzVtx[2];
 159       // RS: for large segment steps d use approximation of cicrular arc s by
 160       // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
 161       // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
 162       // Hence s^2/d^2 = (1+1/6 p^2)^2
 163       dst = dxs*dxs + dys*dys;
 164       if (dst > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
 165         double crv = Abs(esdTr->GetC(bz0));
 166         double fcarc = 1.+crv*crv*dst/6.;
 167         dst *= fcarc*fcarc;
 168       }
 169       dst += dzs*dzs;
 170       dst = Sqrt(dst);
 171     }
 172     //
 173     esdTr->SetStatus(AliESDtrack::kTIME);
 174     //
 175     if (!esdTr->IsOn(AliESDtrack::kITSin)) { // case of tracks w/o ITS prolongation: just set the time integration
 176       dummyTr.StartTimeIntegral();
 177       dummyTr.AddTimeStep(TMath::Sqrt(dst));
 178       dummyTr.GetIntegratedTimes(times);
 179       esdTr->SetIntegratedTimes(times);
 180       esdTr->SetIntegratedLength(dummyTr.GetIntegratedLength());
 181       continue;
 182     }
 183     //
 184     currTr = GetTrackHyp(itr);
 185     currTr->StartTimeIntegral();
 186     currTr->AddTimeStep(dst);
 187     currTr->ResetCovariance(10.);
 188     RefitTrack(currTr,fITS->GetRMax()); // propagate to exit from the ITS/TPC screen
 189     //
 190   }
 191   //
 192   return 0;
 193 }
 194
 195 //_________________________________________________________________________
 196 Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent * /*event*/)
 197 {
 198   //
 199   // To be implemented
 200   //
 201
 202   Info("RefitInward","To be implemented");
 203   return 0;
 204 }
 205
 206 //_________________________________________________________________________
 207 Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
 208 {
 209   // read from tree (if pointer provided) or directly from the ITS reco interface
 210   //
 211   return fReconstructor->LoadClusters(treeRP);
 212 }
 213
 214 //_________________________________________________________________________
 215 void AliITSUTrackerGlo::UnloadClusters()
 216 {
 217   //
 218   // To be implemented
 219   //
 220
 221   Info("UnloadClusters","To be implemented");
 222 }
 223 //_________________________________________________________________________
 224 AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
 225 {
 226   //
 227   // To be implemented
 228   //
 229
 230   Info("GetCluster","To be implemented");
 231   return 0x0;
 232 }
 233
 234 //_________________________________________________________________________
 235 Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
 236 {
 237   // do we need to match this track to ITS?
 238   //
 239   static double bz = GetBz();
 240   if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
 241       esdTr->IsOn(AliESDtrack::kTPCout) ||
 242       esdTr->IsOn(AliESDtrack::kITSin)  ||
 243       esdTr->GetKinkIndex(0)>0) return kFALSE;
 244   //
 245   if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
 246   //
 247   float dtz[2];
 248   esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz);
 249   // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
 250   if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation())
 251        && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
 252            Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
 253   //
 254   return kTRUE;
 255 }
 256
 257 //_________________________________________________________________________
 258 void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
 259 {
 260   // find prolongaion candidates finding for single seed
 261   //
 262   AliITSUSeed seedUC;  // copy of the seed from the upper layer
 263   AliITSUSeed seedT;   // transient seed between the seedUC and new prolongation hypothesis
 264   //
 265   if (!InitHypothesis(esdTr,esdID)) return;  // initialize prolongations hypotheses tree
 266   //
 267   AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
 268   //
 269   TObjArray clArr; // container for transfer of clusters matching to seed
 270   //
 271   int nLrActive = fITS->GetNLayersActive();
 272   for (int ila=nLrActive;ila--;) {
 273     int ilaUp = ila+1;                         // prolong seeds from layer above
 274     //
 275     // for the outermost layer the seed is created from the ESD track
 276     int nSeedsUp = (ilaUp==nLrActive) ? 1 : fCurrHyp->GetNSeeds(ilaUp);
 277     //
 278     for (int isd=0;isd<nSeedsUp;isd++) {
 279       AliITSUSeed* seedU;
 280       if (ilaUp==nLrActive) {
 281         seedU = 0;
 282         seedUC.InitFromESDTrack(esdTr);
 283       }
 284       else {
 285         seedU = fCurrHyp->GetSeed(ilaUp,isd);  // seed on prev.active layer to prolong
 286         seedUC = *seedU;                       // its copy will be prolonged
 287         seedUC.SetParent(seedU);
 288       }
 289       seedUC.ResetFMatrix();                    // reset the matrix for propagation to next layer
 290       // go till next active layer
 291       AliInfo(Form("working on Lr:%d Seed:%d of %d",ila,isd,nSeedsUp));
 292       if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) {
 293         //
 294         AliInfo("Transport failed");
 295         // Check if the seed satisfies to track definition
 296         if (NeedToKill(&seedUC,kTransportFailed) && seedU) seedU->Kill();
 297         continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
 298       }
 299       AliITSURecoLayer* lrA = fITS->GetLayerActive(ila);
 300       if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
 301         if (NeedToKill(&seedUC,kRWCheckFailed) && seedU) seedU->Kill();
 302         continue;
 303       }
 304       int nsens = lrA->FindSensors(&fTrImpData[kTrPhi0], hitSens);  // find detectors which may be hit by the track
 305       AliInfo(Form("Will check %d sensors on lr:%d ",nsens,ila));
 306       //
 307       double bz = GetBz();
 308       for (int isn=nsens;isn--;) {
 309         seedT = seedUC;
 310         AliITSURecoSens* sens = hitSens[isn];
 311         //
 312         // We need to propagate the seed to sensor on lrA staying the frame of the sensor from prev.layer,
 313         // since the transport matrix should be defined in this frame.
 314         double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
 315         if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) continue;
 316         if (!seedT.PropagateToX(xs,bz)) continue;
 317         //      if (!seedT.Rotate(sens->GetPhiTF())) continue;
 318         if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
 319         //
 320         int clID0 = sens->GetFirstClusterId();
 321         for (int icl=sens->GetNClusters();icl--;) {
 322           int res = CheckCluster(&seedT,ila,clID0+icl);
 323           //
 324           if (res==kStopSearchOnSensor) break;     // stop looking on this sensor
 325           if (res==kClusterNotMatching) continue;  // cluster does not match
 326           // cluster is matching and it was added to the hypotheses tree
 327         }
 328       }
 329       // cluster search is done. Do we need ta have a version of this seed skipping current layer
 330       seedT.SetLr(ila);
 331       if (!NeedToKill(&seedT,kMissingCluster)) {
 332         AliITSUSeed* seedSkp = NewSeedFromPool(&seedT);
 333         double penalty = -AliITSUReconstructor::GetRecoParam()->GetMissPenalty(ila);
 334         // to do: make penalty to account for probability to miss the cluster for good reason
 335         seedSkp->SetChi2Cl(penalty);
 336         AddProlongationHypothesis(seedSkp,ila);
 337       }
 338     }
 339     ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
 340     printf(">>> All hypotheses on lr %d: \n",ila);
 341     for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
 342       printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();
 343       //
 344       /*
 345       if (ila!=0) continue;
 346       double vecL[5] = {0};
 347       double matL[15] = {0};
 348       AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
 349       while(sp->GetParent()) {
 350         sp->Smooth(vecL,matL);
 351         if (sp->GetLayerID()>=fITS->GetNLayersActive()-1) break;
 352         sp = (AliITSUSeed*)sp->GetParent();
 353       }
 354       */
 355     }
 356   }
 357   //
 358   SaveCurrentTrackHypotheses();
 359   //
 360 }
 361
 362 //_________________________________________________________________________
 363 Bool_t AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esdID)
 364 {
 365   // init prolongaion candidates finding for single seed
 366   fCurrHyp = GetTrackHyp(esdID);
 367   if (fCurrHyp) return kTRUE;
 368   //
 369   fCurrMass = esdTr->GetMass();
 370   fCurrESDtrack = esdTr;
 371   if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
 372   //
 373   fCurrHyp = new AliITSUTrackHyp(fITS->GetNLayersActive());
 374   fCurrHyp->SetESDTrack(esdTr);
 375   fCurrHyp->SetUniqueID(esdID);
 376   fCurrHyp->SetMass(fCurrMass);
 377   SetTrackHyp(fCurrHyp,esdID);
 378   //
 379   return kTRUE;
 380   // TO DO
 381 }
 382
 383 //_________________________________________________________________________
 384 Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
 385 {
 386   // transport seed from layerFrom to the entrance of layerTo
 387   //
 388   const double kToler = 1e-6; // tolerance for layer on-surface check
 389   //
 390   int dir = lTo > lFrom ? 1:-1;
 391   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 392   Bool_t checkFirst = kTRUE;
 393   while(lFrom!=lTo) {
 394     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 395     if (lrFr) {
 396       Bool_t doLayer = kTRUE;
 397       double xToGo = dir>0 ? lrFr->GetRMax() : lrFr->GetRMin();
 398       if (checkFirst) { // do we need to track till the surface of the current layer ?
 399         checkFirst = kFALSE;
 400         if      (dir>0) { if (curR2-xToGo*xToGo>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
 401         else if (dir<0) { if (xToGo*xToGo-curR2>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
 402       }
 403       if (doLayer) {
 404         if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 405         // go via layer to its boundary, applying material correction.
 406         if (!PropagateSeed(seed,xToGo,fCurrMass, lrFr->GetMaxStep())) return kFALSE;
 407       }
 408     }
 409     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 410     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 411     //
 412     // go the entrance of the layer, assuming no materials in between
 413     double xToGo = dir>0 ? lrTo->GetRMin() : lrTo->GetRMax();
 414     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 415     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) return kFALSE;
 416     lrFr = lrTo;
 417   }
 418   return kTRUE;
 419   //
 420 }
 421
 422 //_________________________________________________________________________
 423 Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo)
 424 {
 425   // transport track from layerFrom to the entrance of layerTo
 426   //
 427   const double kToler = 1e-6; // tolerance for layer on-surface check
 428   //
 429   int dir = lTo > lFrom ? 1:-1;
 430   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 431   Bool_t checkFirst = kTRUE;
 432   while(lFrom!=lTo) {
 433     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 434     if (lrFr) {
 435       Bool_t doLayer = kTRUE;
 436       double xToGo = dir>0 ? lrFr->GetRMax() : lrFr->GetRMin();
 437       if (checkFirst) { // do we need to track till the surface of the current layer ?
 438         checkFirst = kFALSE;
 439         if      (dir>0) { if (curR2-xToGo*xToGo>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
 440         else if (dir<0) { if (xToGo*xToGo-curR2>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
 441       }
 442       if (doLayer) {
 443         if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 444         // go via layer to its boundary, applying material correction.
 445         if (!PropagateSeed(seed,xToGo,fCurrMass, lrFr->GetMaxStep())) return kFALSE;
 446       }
 447     }
 448     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 449     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 450     //
 451     // go the entrance of the layer, assuming no materials in between
 452     double xToGo = dir>0 ? lrTo->GetRMin() : lrTo->GetRMax();
 453     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 454     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) return kFALSE;
 455     lrFr = lrTo;
 456   }
 457   return kTRUE;
 458   //
 459 }
 460
 461 //_________________________________________________________________________
 462 Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
 463 {
 464   // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
 465   // as well as some aux info
 466   double bz = GetBz();
 467   AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
 468   seed->GetXYZ(&fTrImpData[kTrXIn]);    // lab position at the entrance from above
 469   static AliExternalTrackParam sc;   // seed copy for manipalitions
 470   sc = *seed;
 471   //
 472   fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
 473   if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
 474   double dr  = lrA->GetDR();                              // approximate X dist at the inner radius
 475   if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
 476     // special case: track does not reach inner radius, might be tangential
 477     double r = sc.GetD(0,0,bz);
 478     double x;
 479     if (!sc.GetXatLabR(r,x,bz,-1)) {
 480       sc.Print();
 481       AliFatal(Form("This should not happen: r=%f",r));
 482     }
 483     dr = Abs(sc.GetX() - x);
 484     if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) {
 485       sc.Print();
 486       AliFatal(Form("This should not happen: x=%f",x));
 487     }
 488   }
 489   //
 490   fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
 491   double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
 492   double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
 493   sgy = Sqrt(sgy)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY();
 494   sgz = Sqrt(sgz)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ();
 495   fTrImpData[kTrPhi0] = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
 496   fTrImpData[kTrZ0]   = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
 497   fTrImpData[kTrDPhi] = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]) + sgy/lrA->GetR();
 498   fTrImpData[kTrDZ]   = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn])   + sgz;
 499   //
 500   return kTRUE;
 501 }
 502
 503 //_________________________________________________________________________
 504 AliITSUSeed* AliITSUTrackerGlo::NewSeedFromPool(const AliITSUSeed* src)
 505 {
 506   // create new seed, optionally copying from the source
 507   return src ?
 508     new(fSeedsPool[fSeedsPool.GetEntriesFast()]) AliITSUSeed(*src) :
 509     new(fSeedsPool[fSeedsPool.GetEntriesFast()]) AliITSUSeed();
 510 }
 511
 512 //_________________________________________________________________________
 513 Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
 514 {
 515   // Check if the cluster (in tracking frame!) is matching to track.
 516   // The track must be already propagated to sensor tracking frame.
 517   // Returns:  kStopSearchOnSensor if the search on given sensor should be stopped,
 518   //           kClusterMatching    if the cluster is matching
 519   //           kClusterMatching    otherwise
 520   //
 521   const double kTolerX = 5e-4;
 522   AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
 523   //
 524   Bool_t goodCl = kFALSE;
 525   int currLabel = Abs(fCurrESDtrack->GetTPCLabel());
 526   //
 527   if (cl->GetLabel(0)>=0) {for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}}
 528   else goodCl = kTRUE;
 529   //
 530   if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
 531     if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
 532       if (goodCl) {printf("Loose good cl: Failed propagation. |"); cl->Print();}
 533       return kStopSearchOnSensor; // propagation failed, seedT is intact
 534     }
 535   }
 536   double dy = cl->GetY()-track->GetY();
 537   double dz = cl->GetZ()-track->GetZ();
 538   //
 539   double dy2 = dy*dy;
 540   double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
 541     AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY(); // RS TOOPTIMIZE
 542   if (dy2>tol2) {                          // the clusters are sorted in Z(col) then in Y(row).
 543     if (goodCl) {printf("Loose good cl: dy2=%e > tol2=%e |",dy2,tol2); cl->Print();}
 544     if (dy>0) return kStopSearchOnSensor;  // No chance that other cluster of this sensor will match (all Y's will be even larger)
 545     else      return kClusterNotMatching;   // Other clusters may match
 546   }
 547   double dz2 = dz*dz;
 548   tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
 549     AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ(); // RS TOOPTIMIZE
 550   if (dz2>tol2) {
 551     if (goodCl) {printf("Loose good cl: dz2=%e > tol2=%e |",dz2,tol2); cl->Print();}
 552     return kClusterNotMatching; // Other clusters may match
 553   }
 554   //
 555   // check chi2
 556   Double_t p[2]={cl->GetY(), cl->GetZ()};
 557   Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
 558   double chi2 = track->GetPredictedChi2(p,cov);
 559   if (chi2>AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr)) {
 560     if (goodCl) {
 561       printf("Loose good cl: Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e\n",
 562              chi2,AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr),dy,dz);
 563       track->Print("etp");
 564       cl->Print("");
 565     }
 566     return kClusterNotMatching;
 567   }
 568   //
 569   track = NewSeedFromPool(track);  // input track will be reused, use its clone for updates
 570   if (!track->Update()) {
 571     if (goodCl) {printf("Loose good cl: Failed update |"); cl->Print();}
 572     return kClusterNotMatching;
 573   }
 574   track->SetChi2Cl(chi2);
 575   track->SetLrClusterID(lr,clID);
 576   cl->IncreaseClusterUsage();
 577   //
 578   track->SetFake(!goodCl);
 579   //
 580   AliInfo(Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)",
 581                goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
 582   //
 583   AddProlongationHypothesis(track,lr);
 584   //
 585   return kClusterMatching;
 586 }
 587
 588 //_________________________________________________________________________
 589 Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
 590 {
 591   // check if the seed should not be discarded
 592   const UShort_t kMask = 0xffff;
 593   if (flag==kMissingCluster) {
 594     int lastChecked = seed->GetLayerID();
 595     UShort_t patt = seed->GetHitsPattern();
 596     if (lastChecked) patt |= ~(kMask<<lastChecked); // not all layers were checked, complete unchecked once by potential hits
 597     Bool_t seedOK = fTrCond.CheckPattern(patt);
 598     return !seedOK;
 599   }
 600   return kTRUE;
 601 }
 602
 603 //______________________________________________________________________________
 604 Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
 605 {
 606   // propagate seed to given x applying material correction if requested
 607   const Double_t kEpsilon = 1e-5;
 608   Double_t xpos     = seed->GetX();
 609   Int_t dir         = (xpos<xToGo) ? 1:-1;
 610   Double_t xyz0[3],xyz1[3],param[7];
 611   //
 612   if (matCorr) seed->GetXYZ(xyz1);   //starting global position
 613   while ( (xToGo-xpos)*dir > kEpsilon){
 614     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 615     Double_t x    = xpos+step;
 616     Double_t bz=GetBz();   // getting the local Bz
 617     if (!seed->PropagateToX(x,bz))  return kFALSE;
 618     if (matCorr) {
 619       xyz0[0]=xyz1[0]; // global pos at the beginning of step
 620       xyz0[1]=xyz1[1];
 621       xyz0[2]=xyz1[2];
 622       seed->GetXYZ(xyz1);    //  // global pos at the end of step
 623       MeanMaterialBudget(xyz0,xyz1,param);
 624       Double_t xrho=param[0]*param[4], xx0=param[1];
 625       if (dir>0) xrho = -xrho; // outward should be negative
 626       if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) return kFALSE;
 627     }
 628     xpos = seed->GetX();
 629   }
 630   return kTRUE;
 631 }
 632
 633 //______________________________________________________________________________
 634 Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
 635 {
 636   // propagate seed to given x applying material correction if requested
 637   const Double_t kEpsilon = 1e-5;
 638   Double_t xpos     = seed->GetX();
 639   Int_t dir         = (xpos<xToGo) ? 1:-1;
 640   Double_t xyz0[3],xyz1[3],param[7];
 641   //
 642   Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
 643   if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
 644   while ( (xToGo-xpos)*dir > kEpsilon){
 645     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 646     Double_t x    = xpos+step;
 647     Double_t bz=GetBz();   // getting the local Bz
 648     if (!seed->PropagateTo(x,bz))  return kFALSE;
 649     double ds = 0;
 650     if (matCorr || updTime) {
 651       xyz0[0]=xyz1[0]; // global pos at the beginning of step
 652       xyz0[1]=xyz1[1];
 653       xyz0[2]=xyz1[2];
 654       seed->GetXYZ(xyz1);    //  // global pos at the end of step
 655       //
 656       if (matCorr) {
 657         MeanMaterialBudget(xyz0,xyz1,param);
 658         Double_t xrho=param[0]*param[4], xx0=param[1];
 659         if (dir>0) xrho = -xrho; // outward should be negative
 660         if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 661         ds = param[4];
 662       }
 663       else { // matCorr is not requested but time integral is
 664         double d0 = xyz1[0]-xyz0[0];
 665         double d1 = xyz1[1]-xyz0[1];
 666         double d2 = xyz1[2]-xyz0[2];
 667         ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
 668       }
 669     }
 670     if (updTime) seed->AddTimeStep(ds);
 671     //
 672     xpos = seed->GetX();
 673   }
 674   return kTRUE;
 675 }
 676
 677 //______________________________________________________________________________
 678 void AliITSUTrackerGlo::SaveCurrentTrackHypotheses()
 679 {
 680   // RS: shall we clean up killed seeds?
 681   fCurrHyp = 0;
 682   // TODO
 683
 684 }
 685
 686 //______________________________________________________________________________
 687 void AliITSUTrackerGlo::FinalizeHypotheses()
 688 {
 689   // select winner for each hypothesis, remove cl. sharing conflicts
 690   AliInfo("TODO");
 691   //
 692   int nh = fHypStore.GetEntriesFast();
 693   for (int ih=0;ih<nh;ih++) {
 694     AliITSUTrackHyp* hyp = (AliITSUTrackHyp*) fHypStore.UncheckedAt(ih);
 695     if (!hyp) continue;
 696     hyp->DefineWinner();  // TODO
 697     UpdateESDTrack(hyp);
 698   }
 699
 700 }
 701
 702 //______________________________________________________________________________
 703 void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp)
 704 {
 705   // update ESD track with current best hypothesis
 706   AliESDtrack* esdTr = hyp->GetESDTrack();
 707   if (!esdTr) return;
 708   AliITSUSeed* win = hyp->GetWinner();
 709   if (!win) return;
 710   esdTr->UpdateTrackParams(hyp,AliESDtrack::kITSin);
 711   //
 712   // transfer module indices
 713   // TODO
 714 }
 715
 716 //______________________________________________________________________________
 717 Bool_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest)
 718 {
 719   // refit track till radius rDest
 720   double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
 721   int dir,lrStart,lrStop;
 722   //
 723   if (rCurr<rDest) {
 724     dir = 1;
 725     lrStart = 0;
 726     lrStop  = fITS->GetNLayers()-1;
 727   }
 728   else {
 729     dir = 1;
 730     lrStart = fITS->GetNLayers()-1;
 731     lrStop  = 0;
 732   }
 733   //
 734   for (int ilr=lrStart;ilr!=lrStop;ilr+=dir) {
 735     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
 736     if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
 737     //
 738   }
 739
 740 }