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
  44 //----------------- tmp stuff -----------------
  45
  46 ClassImp(AliITSUTrackerGlo)
  47
  48 const Double_t AliITSUTrackerGlo::fgkToler =  1e-6;// tolerance for layer on-surface check
  49
  50
  51 //_________________________________________________________________________
  52 AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
  53 :  fReconstructor(rec)
  54   ,fITS(0)
  55   ,fCurrESDtrack(0)
  56   ,fCurrMass(kPionMass)
  57   ,fHypStore(100)
  58   ,fCurrHyp(0)
  59   ,fSeedsPool("AliITSUSeed",0)
  60   ,fFreeSeedsID(0)
  61   ,fNFreeSeeds(0)
  62   ,fLastSeedID(0)
  63   ,fDefTrackConds(0)
  64   ,fCurrTrackCond(0)
  65   ,fTrackPhase(-1)
  66   ,fClInfo(0)
  67 {
  68   // Default constructor
  69   if (rec) Init(rec);
  70 }
  71
  72 //_________________________________________________________________________
  73 AliITSUTrackerGlo::~AliITSUTrackerGlo()
  74 {
  75  // Default destructor
  76  //
  77   delete fITS;
  78   delete[] fClInfo;
  79   //
  80 }
  81
  82 //_________________________________________________________________________
  83 void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
  84 {
  85   // init with external reconstructor
  86   //
  87   fITS = new AliITSURecoDet(rec->GetGeom(),"ITSURecoInterface");
  88   int nLr = fITS->GetNLayersActive();
  89   fClInfo = new Int_t[nLr<<1];
  90   for (int ilr=nLr;ilr--;) {
  91     fITS->GetLayerActive(ilr)->SetClusters(rec->GetClusters(ilr));
  92   }
  93   //
  94   fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
  95   fFreeSeedsID.Set(1000);
  96   //
  97   /*
  98   fTrCond.SetNLayers(fITS->GetNLayersActive());
  99   fTrCond.AddNewCondition(5);
 100   fTrCond.AddGroupPattern( (0x1<<0)|(0x1<<1) );
 101   fTrCond.AddGroupPattern( (0x1<<3)|(0x1<<4) );
 102   fTrCond.AddGroupPattern( (0x1<<5)|(0x1<<6) );
 103   //
 104   fTrCond.AddNewCondition(5);
 105   fTrCond.AddGroupPattern( (0x1<<0)|(0x1<<2) );
 106   fTrCond.AddGroupPattern( (0x1<<3)|(0x1<<4) );
 107   fTrCond.AddGroupPattern( (0x1<<5)|(0x1<<6) );
 108   //
 109   fTrCond.AddNewCondition(5);
 110   fTrCond.AddGroupPattern( (0x1<<1)|(0x1<<2) );
 111   fTrCond.AddGroupPattern( (0x1<<3)|(0x1<<4) );
 112   fTrCond.AddGroupPattern( (0x1<<5)|(0x1<<6) );
 113   //
 114   printf("Tracking Conditions: ");
 115   fTrCond.Print();
 116   */
 117 }
 118
 119 //_________________________________________________________________________
 120 void AliITSUTrackerGlo::CreateDefaultTrackCond()
 121 {
 122   // creates default tracking conditions to be used when recoparam does not provide them
 123   int nLr = fITS->GetNLayersActive();
 124   fClInfo = new Int_t[nLr<<1];
 125   //
 126   AliITSUTrackCond* cond = new AliITSUTrackCond();
 127   //
 128   cond->SetNLayers(fITS->GetNLayersActive());
 129   cond->AddNewCondition(nLr);
 130   cond->AddGroupPattern( 0xffff ); // require all layers hit
 131   //
 132   fDefTrackConds.AddLast(cond);
 133   //
 134   AliInfo("Created conditions: ");
 135   for (int i=0;i<fDefTrackConds.GetEntriesFast();i++) fDefTrackConds[i]->Print();
 136   //
 137 }
 138
 139
 140 //_________________________________________________________________________
 141 Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
 142 {
 143   //
 144   TObjArray *trackConds = 0;
 145   //
 146   SetTrackingPhase(kClus2Tracks);
 147   ResetSeedsPool();
 148   int nTrESD = esdEv->GetNumberOfTracks();
 149   AliInfo(Form("Will try to find prolongations for %d tracks",nTrESD));
 150   AliITSUReconstructor::GetRecoParam()->Print();
 151   fHypStore.Delete();
 152   if (fHypStore.GetSize()<nTrESD) fHypStore.Expand(nTrESD+100);
 153   //
 154   fITS->ProcessClusters();
 155   //
 156   int nTrackCond = AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();
 157   if (nTrackCond<1) {
 158     if (!fDefTrackConds.GetEntriesFast()) {
 159       AliInfo("No tracking conditions found in recoparams, creating default one requesting all layers hit");
 160       CreateDefaultTrackCond();
 161     }
 162     trackConds = &fDefTrackConds;
 163     nTrackCond = trackConds->GetEntriesFast();
 164   }
 165   else trackConds = AliITSUReconstructor::GetRecoParam()->GetTrackingConditions();
 166   //
 167   for (int icnd=0;icnd<nTrackCond;icnd++) {
 168     fCurrTrackCond = (AliITSUTrackCond*)trackConds->UncheckedAt(icnd);
 169     // select ESD tracks to propagate
 170     for (int itr=0;itr<nTrESD;itr++) {
 171       AliESDtrack *esdTr = esdEv->GetTrack(itr);
 172       AliInfo(Form("Processing track %d | MCLabel: %d",itr,esdTr->GetTPCLabel()));
 173       if (!NeedToProlong(esdTr)) continue;  // are we interested in this track?
 174       FindTrack(esdTr, itr);
 175     }
 176     //
 177     AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
 178     fHypStore.Print();
 179     FinalizeHypotheses();
 180   }
 181   //
 182   return 0;
 183 }
 184
 185 //_________________________________________________________________________
 186 Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent *esdEv)
 187 {
 188   //
 189   // Do outward fits in ITS
 190   //
 191   SetTrackingPhase(kPropBack);
 192   int nTrESD = esdEv->GetNumberOfTracks();
 193   AliInfo(Form("Will propagate back %d tracks",nTrESD));
 194   //
 195   double bz0 = GetBz();
 196   Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
 197   AliITSUTrackHyp dummyTr,*currTr=0;
 198   const double kWatchStep=10.; // for larger steps watch arc vs segment difference
 199   Double_t times[AliPID::kSPECIES];
 200   //
 201   //
 202   for (int itr=0;itr<nTrESD;itr++) {
 203     AliESDtrack *esdTr = esdEv->GetTrack(itr);
 204     // Start time integral and add distance from current position to vertex
 205     if (esdTr->IsOn(AliESDtrack::kITSout)) continue; // already done
 206     //
 207     esdTr->GetXYZ(xyzTrk);
 208     Double_t dst = 0.;     // set initial track lenght, tof
 209     {
 210       double dxs = xyzTrk[0] - xyzVtx[0];
 211       double dys = xyzTrk[1] - xyzVtx[1];
 212       double dzs = xyzTrk[2] - xyzVtx[2];
 213       // RS: for large segment steps d use approximation of cicrular arc s by
 214       // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
 215       // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
 216       // Hence s^2/d^2 = (1+1/6 p^2)^2
 217       dst = dxs*dxs + dys*dys;
 218       if (dst > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
 219         double crv = Abs(esdTr->GetC(bz0));
 220         double fcarc = 1.+crv*crv*dst/6.;
 221         dst *= fcarc*fcarc;
 222       }
 223       dst += dzs*dzs;
 224       dst = Sqrt(dst);
 225     }
 226     //
 227     esdTr->SetStatus(AliESDtrack::kTIME);
 228     //
 229     if (!esdTr->IsOn(AliESDtrack::kITSin)) { // case of tracks w/o ITS prolongation: just set the time integration
 230       dummyTr.AliExternalTrackParam::operator=(*esdTr);
 231       dummyTr.StartTimeIntegral();
 232       dummyTr.AddTimeStep(dst);
 233       dummyTr.GetIntegratedTimes(times);
 234       esdTr->SetIntegratedTimes(times);
 235       esdTr->SetIntegratedLength(dummyTr.GetIntegratedLength());
 236       continue;
 237     }
 238     //
 239     currTr = GetTrackHyp(itr);
 240     currTr->StartTimeIntegral();
 241     currTr->AddTimeStep(dst);
 242     printf("Before resetCov: "); currTr->AliExternalTrackParam::Print();
 243     currTr->ResetCovariance(10000);
 244     if (RefitTrack(currTr,fITS->GetRMax())) { // propagate to exit from the ITS/TPC screen
 245       UpdateESDTrack(currTr,AliESDtrack::kITSout);
 246     }
 247     else {
 248       AliInfo(Form("Refit Failed for track %d",itr));
 249     }
 250     //
 251   }
 252   //
 253   return 0;
 254 }
 255
 256 //_________________________________________________________________________
 257 Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent *esdEv)
 258 {
 259   //
 260   // refit the tracks inward, using current cov.matrix
 261   //
 262   SetTrackingPhase(kRefitInw);
 263   Int_t nTrESD = esdEv->GetNumberOfTracks();
 264   AliInfo(Form("Will refit inward %d tracks",nTrESD));
 265   AliITSUTrackHyp *currTr=0;
 266   //
 267   for (int itr=0;itr<nTrESD;itr++) {
 268     AliESDtrack *esdTr = esdEv->GetTrack(itr);
 269     // Start time integral and add distance from current position to vertex
 270     UInt_t trStat = esdTr->GetStatus();
 271     if ( !(trStat & AliESDtrack::kITSout) ) continue;
 272     if (   trStat & AliESDtrack::kITSrefit ) continue; // already done
 273     if (  (trStat & AliESDtrack::kTPCout) && !(trStat & AliESDtrack::kTPCrefit) ) continue;
 274     //
 275     currTr = GetTrackHyp(itr);
 276     currTr->AliExternalTrackParam::operator=(*esdTr);  // fetch current ESDtrack kinematics
 277     if (RefitTrack(currTr,fITS->GetRMin())) { // propagate up to inside radius of the beam pipe
 278       UpdateESDTrack(currTr,AliESDtrack::kITSrefit);
 279     }
 280     else {
 281       AliInfo(Form("Refit Failed for track %d",itr));
 282     }
 283   }
 284   //
 285   return 0;
 286 }
 287
 288 //_________________________________________________________________________
 289 Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
 290 {
 291   // read from tree (if pointer provided) or directly from the ITS reco interface
 292   //
 293   return fReconstructor->LoadClusters(treeRP);
 294 }
 295
 296 //_________________________________________________________________________
 297 void AliITSUTrackerGlo::UnloadClusters()
 298 {
 299   //
 300   // To be implemented
 301   //
 302
 303   Info("UnloadClusters","To be implemented");
 304 }
 305 //_________________________________________________________________________
 306 AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
 307 {
 308   //
 309   // To be implemented
 310   //
 311
 312   Info("GetCluster","To be implemented");
 313   return 0x0;
 314 }
 315
 316 //_________________________________________________________________________
 317 Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
 318 {
 319   // do we need to match this track to ITS?
 320   //
 321   static double bz = GetBz();
 322   if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
 323       esdTr->IsOn(AliESDtrack::kTPCout) ||
 324       esdTr->IsOn(AliESDtrack::kITSin)  ||
 325       esdTr->GetKinkIndex(0)>0) return kFALSE;
 326   //
 327   if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
 328   //
 329   float dtz[2];
 330   esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz);
 331   // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
 332   if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation())
 333        && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
 334            Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
 335   //
 336   return kTRUE;
 337 }
 338
 339 //_________________________________________________________________________
 340 void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
 341 {
 342   // find prolongaion candidates finding for single seed
 343   //
 344   AliITSUSeed seedUC;  // copy of the seed from the upper layer
 345   AliITSUSeed seedT;   // transient seed between the seedUC and new prolongation hypothesis
 346   //
 347   if (!InitHypothesis(esdTr,esdID)) return;  // initialize prolongations hypotheses tree
 348   //
 349   AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
 350   //
 351   TObjArray clArr; // container for transfer of clusters matching to seed
 352   //
 353   int nLrActive = fITS->GetNLayersActive();
 354   for (int ila=nLrActive;ila--;) {
 355     int ilaUp = ila+1;                         // prolong seeds from layer above
 356     //
 357     // for the outermost layer the seed is created from the ESD track
 358     int nSeedsUp = (ilaUp==nLrActive) ? 1 : fCurrHyp->GetNSeeds(ilaUp);
 359     //
 360     for (int isd=0;isd<nSeedsUp;isd++) {
 361       AliITSUSeed* seedU;
 362       if (ilaUp==nLrActive) {
 363         seedU = 0;
 364         seedUC.InitFromESDTrack(esdTr);
 365       }
 366       else {
 367         seedU = fCurrHyp->GetSeed(ilaUp,isd);  // seed on prev.active layer to prolong
 368         seedUC = *seedU;                       // its copy will be prolonged
 369         seedUC.SetParent(seedU);
 370       }
 371       seedUC.ResetFMatrix();                    // reset the matrix for propagation to next layer
 372       // go till next active layer
 373       AliInfo(Form("working on Lr:%d Seed:%d of %d",ila,isd,nSeedsUp));
 374       if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) {
 375         //
 376         AliInfo("Transport failed");
 377         // Check if the seed satisfies to track definition
 378         if (NeedToKill(&seedUC,kTransportFailed) && seedU) KillSeed(seedU,kTRUE);
 379         continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
 380       }
 381       AliITSURecoLayer* lrA = fITS->GetLayerActive(ila);
 382       if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
 383         if (NeedToKill(&seedUC,kRWCheckFailed) && seedU) KillSeed(seedU,kTRUE);
 384         continue;
 385       }
 386       int nsens = lrA->FindSensors(&fTrImpData[kTrPhi0], hitSens);  // find detectors which may be hit by the track
 387       AliInfo(Form("Will check %d sensors on lr:%d ",nsens,ila));
 388       //
 389       double bz = GetBz();
 390       for (int isn=nsens;isn--;) {
 391         seedT = seedUC;
 392         AliITSURecoSens* sens = hitSens[isn];
 393         //
 394         // We need to propagate the seed to sensor on lrA staying the frame of the sensor from prev.layer,
 395         // since the transport matrix should be defined in this frame.
 396         double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
 397         if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) continue;
 398         if (!PropagateSeed(&seedT,xs,fCurrMass)) continue;
 399         //      if (!seedT.PropagateToX(xs,bz)) continue;
 400         //      if (!seedT.Rotate(sens->GetPhiTF())) continue;
 401         if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
 402         //
 403         int clID0 = sens->GetFirstClusterId();
 404         for (int icl=sens->GetNClusters();icl--;) {
 405           int res = CheckCluster(&seedT,ila,clID0+icl);
 406           //
 407           if (res==kStopSearchOnSensor) break;     // stop looking on this sensor
 408           if (res==kClusterNotMatching) continue;  // cluster does not match
 409           // cluster is matching and it was added to the hypotheses tree
 410         }
 411       }
 412       // cluster search is done. Do we need ta have a version of this seed skipping current layer
 413       seedT.SetLr(ila);
 414       if (!NeedToKill(&seedT,kMissingCluster)) {
 415         AliITSUSeed* seedSkp = NewSeedFromPool(&seedT);
 416         double penalty = -AliITSUReconstructor::GetRecoParam()->GetMissPenalty(ila);
 417         // to do: make penalty to account for probability to miss the cluster for good reason
 418         seedSkp->SetChi2Cl(penalty);
 419         AddProlongationHypothesis(seedSkp,ila);
 420       }
 421     }
 422     ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
 423     printf(">>> All hypotheses on lr %d: \n",ila);
 424     for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
 425       printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();
 426       //
 427       /*
 428       if (ila!=0) continue;
 429       double vecL[5] = {0};
 430       double matL[15] = {0};
 431       AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
 432       while(sp->GetParent()) {
 433         sp->Smooth(vecL,matL);
 434         if (sp->GetLayerID()>=fITS->GetNLayersActive()-1) break;
 435         sp = (AliITSUSeed*)sp->GetParent();
 436       }
 437       */
 438     }
 439   }
 440   //
 441   SaveCurrentTrackHypotheses();
 442   //
 443 }
 444
 445 //_________________________________________________________________________
 446 Bool_t AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esdID)
 447 {
 448   // init prolongaion candidates finding for single seed
 449   fCurrHyp = GetTrackHyp(esdID);
 450   if (fCurrHyp) return kTRUE;
 451   //
 452   fCurrMass = esdTr->GetMass();
 453   fCurrESDtrack = esdTr;
 454   if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
 455   //
 456   fCurrHyp = new AliITSUTrackHyp(fITS->GetNLayersActive());
 457   fCurrHyp->SetESDTrack(esdTr);
 458   fCurrHyp->SetUniqueID(esdID);
 459   fCurrHyp->SetMass(fCurrMass);
 460   SetTrackHyp(fCurrHyp,esdID);
 461   //
 462   return kTRUE;
 463   // TO DO
 464 }
 465
 466 //_________________________________________________________________________
 467 Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
 468 {
 469   // transport seed from layerFrom to the entrance of layerTo
 470   //
 471   //
 472   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
 473   //
 474   int dir = lTo > lFrom ? 1:-1;
 475   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 476   Bool_t checkFirst = kTRUE;
 477   while(lFrom!=lTo) {
 478     if (lrFr) {
 479       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
 480       checkFirst = kFALSE;
 481     }
 482     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 483     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 484     //
 485     // go the entrance of the layer, assuming no materials in between
 486     double xToGo = lrTo->GetR(-dir);
 487     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 488     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) return kFALSE;
 489     lrFr = lrTo;
 490   }
 491   return kTRUE;
 492   //
 493 }
 494
 495 //_________________________________________________________________________
 496 Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo)
 497 {
 498   // transport track from layerFrom to the entrance of layerTo
 499   //
 500   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
 501   //
 502   int dir = lTo > lFrom ? 1:-1;
 503   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 504   Bool_t checkFirst = kTRUE;
 505   while(lFrom!=lTo) {
 506     if (lrFr) {
 507       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
 508       checkFirst = kFALSE;
 509     }
 510     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 511     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 512     //
 513     if (!GoToExitFromLayer(seed,lrTo,dir)) return kFALSE; // go the entrance of the layer, assuming no materials in between
 514     lrFr = lrTo;
 515   }
 516   return kTRUE;
 517   //
 518 }
 519
 520 //_________________________________________________________________________
 521 Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 522 {
 523   // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
 524   // If check is requested, do this only provided the track has not exited the layer already
 525   double xToGo = lr->GetR(dir);
 526   if (check) { // do we need to track till the surface of the current layer ?
 527     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 528     //    AliInfo(Form(" dir:%d Cur: %e Tgt: %e",dir,Sqrt(curR2),xToGo));
 529     if      (dir>0) { if (curR2-xToGo*xToGo>fgkToler) return kTRUE; } // on the surface or outside of the layer
 530     else if (dir<0) { if (xToGo*xToGo-curR2>fgkToler) return kTRUE; } // on the surface or outside of the layer
 531   }
 532   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 533   // go via layer to its boundary, applying material correction.
 534   if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
 535   //
 536   return kTRUE;
 537   //
 538 }
 539
 540 //_________________________________________________________________________
 541 Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 542 {
 543   // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
 544   // If check is requested, do this only provided the track has not exited the layer already
 545   double xToGo = lr->GetR(dir);
 546   if (check) { // do we need to track till the surface of the current layer ?
 547     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 548     if      (dir>0) { if (curR2-xToGo*xToGo>fgkToler) return kTRUE; } // on the surface or outside of the layer
 549     else if (dir<0) { if (xToGo*xToGo-curR2>fgkToler) return kTRUE; } // on the surface or outside of the layer
 550   }
 551   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 552   // go via layer to its boundary, applying material correction.
 553   if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
 554   //
 555   return kTRUE;
 556   //
 557 }
 558
 559
 560 //_________________________________________________________________________
 561 Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 562 {
 563   // go to the entrance of lr in direction dir, w/o applying material corrections.
 564   // If check is requested, do this only provided the track did not reach the layer already
 565   double xToGo = lr->GetR(-dir);
 566   if (check) { // do we need to track till the surface of the current layer ?
 567     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 568     if      (dir>0) { if (curR2-xToGo*xToGo>fgkToler) return kTRUE; } // already passed
 569     else if (dir<0) { if (xToGo*xToGo-curR2>fgkToler) return kTRUE; } // already passed
 570   }
 571   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 572   // go via layer to its boundary, applying material correction.
 573   if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
 574   return kTRUE;
 575   //
 576 }
 577
 578 //_________________________________________________________________________
 579 Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 580 {
 581   // go to the entrance of lr in direction dir, w/o applying material corrections.
 582   // If check is requested, do this only provided the track did not reach the layer already
 583   double xToGo = lr->GetR(-dir);
 584   if (check) { // do we need to track till the surface of the current layer ?
 585     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 586     if      (dir>0) { if (curR2-xToGo*xToGo>fgkToler) return kTRUE; } // already passed
 587     else if (dir<0) { if (xToGo*xToGo-curR2>fgkToler) return kTRUE; } // already passed
 588   }
 589   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 590   // go via layer to its boundary, applying material correction.
 591   if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
 592   return kTRUE;
 593   //
 594 }
 595
 596 //_________________________________________________________________________
 597 Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
 598 {
 599   // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
 600   // as well as some aux info
 601   double bz = GetBz();
 602   AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
 603   seed->GetXYZ(&fTrImpData[kTrXIn]);    // lab position at the entrance from above
 604   static AliExternalTrackParam sc;   // seed copy for manipalitions
 605   sc = *seed;
 606   //
 607   fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
 608   if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
 609   double dr  = lrA->GetDR();                              // approximate X dist at the inner radius
 610   if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
 611     // special case: track does not reach inner radius, might be tangential
 612     double r = sc.GetD(0,0,bz);
 613     double x;
 614     if (!sc.GetXatLabR(r,x,bz,-1)) {
 615       sc.Print();
 616       AliFatal(Form("This should not happen: r=%f",r));
 617     }
 618     dr = Abs(sc.GetX() - x);
 619     if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) {
 620       sc.Print();
 621       AliFatal(Form("This should not happen: x=%f",x));
 622     }
 623   }
 624   //
 625   fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
 626   double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
 627   double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
 628   sgy = Sqrt(sgy)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY();
 629   sgz = Sqrt(sgz)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ();
 630   fTrImpData[kTrPhi0] = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
 631   fTrImpData[kTrZ0]   = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
 632   fTrImpData[kTrDPhi] = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]) + sgy/lrA->GetR();
 633   fTrImpData[kTrDZ]   = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn])   + sgz;
 634   //
 635   return kTRUE;
 636 }
 637
 638 //________________________________________
 639 void AliITSUTrackerGlo::ResetSeedsPool()
 640 {
 641   // mark all seeds in the pool as unused
 642   AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
 643   fNFreeSeeds = 0;
 644   fSeedsPool.Clear(); // seeds don't allocate memory
 645 }
 646
 647
 648 //________________________________________
 649 void AliITSUTrackerGlo::MarkSeedFree(AliITSUSeed *sd)
 650 {
 651   // account that this seed is "deleted"
 652   int id = sd->GetPoolID();
 653   if (id<0) {
 654     AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
 655     return;
 656   }
 657   //  AliInfo(Form("%d %p",id, seed));
 658   fSeedsPool.RemoveAt(id);
 659   if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
 660   fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
 661 }
 662
 663 //_________________________________________________________________________
 664 Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
 665 {
 666   // Check if the cluster (in tracking frame!) is matching to track.
 667   // The track must be already propagated to sensor tracking frame.
 668   // Returns:  kStopSearchOnSensor if the search on given sensor should be stopped,
 669   //           kClusterMatching    if the cluster is matching
 670   //           kClusterMatching    otherwise
 671   //
 672   const double kTolerX = 5e-4;
 673   AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
 674   //
 675   Bool_t goodCl = kFALSE;
 676   int currLabel = Abs(fCurrESDtrack->GetTPCLabel());
 677   //
 678   if (cl->GetLabel(0)>=0) {
 679     for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}
 680   }
 681   else goodCl = kTRUE;
 682   //
 683   if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
 684     if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
 685       if (goodCl) {printf("Loose good cl: Failed propagation. |"); cl->Print();}
 686       return kStopSearchOnSensor; // propagation failed, seedT is intact
 687     }
 688   }
 689   double dy = cl->GetY()-track->GetY();
 690   double dz = cl->GetZ()-track->GetZ();
 691   //
 692   double dy2 = dy*dy;
 693   double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
 694     AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY(); // RS TOOPTIMIZE
 695   if (dy2>tol2) {                          // the clusters are sorted in Z(col) then in Y(row).
 696     if (goodCl) {printf("Loose good cl: dy2=%e > tol2=%e |",dy2,tol2); cl->Print();}
 697     if (dy>0) return kStopSearchOnSensor;  // No chance that other cluster of this sensor will match (all Y's will be even larger)
 698     else      return kClusterNotMatching;   // Other clusters may match
 699   }
 700   double dz2 = dz*dz;
 701   tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
 702     AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ(); // RS TOOPTIMIZE
 703   if (dz2>tol2) {
 704     if (goodCl) {printf("Loose good cl: dz2=%e > tol2=%e |",dz2,tol2); cl->Print();}
 705     return kClusterNotMatching; // Other clusters may match
 706   }
 707   //
 708   // check chi2
 709   Double_t p[2]={cl->GetY(), cl->GetZ()};
 710   Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
 711   double chi2 = track->GetPredictedChi2(p,cov);
 712   if (chi2>AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr)) {
 713     if (goodCl) {
 714       printf("Loose good cl: Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e\n",
 715              chi2,AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr),dy,dz);
 716       track->Print("etp");
 717       cl->Print("");
 718     }
 719     return kClusterNotMatching;
 720   }
 721   //
 722   track = NewSeedFromPool(track);  // input track will be reused, use its clone for updates
 723   if (!track->Update()) {
 724     if (goodCl) {printf("Loose good cl: Failed update |"); cl->Print();}
 725     MarkSeedFree(track);
 726     return kClusterNotMatching;
 727   }
 728   track->SetChi2Cl(chi2);
 729   track->SetLrClusterID(lr,clID);
 730   cl->IncreaseClusterUsage();
 731   //
 732   track->SetFake(!goodCl);
 733   //
 734   AliInfo(Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)",
 735                goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
 736   //
 737   AddProlongationHypothesis(track,lr);
 738   //
 739   return kClusterMatching;
 740 }
 741
 742 //_________________________________________________________________________
 743 Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
 744 {
 745   // check if the seed should not be discarded
 746   const UShort_t kMask = 0xffff;
 747   if (flag==kMissingCluster) {
 748     int lastChecked = seed->GetLayerID();
 749     UShort_t patt = seed->GetHitsPattern();
 750     if (lastChecked) patt |= ~(kMask<<lastChecked); // not all layers were checked, complete unchecked once by potential hits
 751     Bool_t seedOK = fCurrTrackCond->CheckPattern(patt);
 752     return !seedOK;
 753   }
 754   return kTRUE;
 755 }
 756
 757 //______________________________________________________________________________
 758 Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
 759 {
 760   // propagate seed to given x applying material correction if requested
 761   const Double_t kEpsilon = 1e-5;
 762   Double_t xpos     = seed->GetX();
 763   Int_t dir         = (xpos<xToGo) ? 1:-1;
 764   Double_t xyz0[3],xyz1[3],param[7];
 765   //
 766   Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
 767   if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
 768   while ( (xToGo-xpos)*dir > kEpsilon){
 769     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 770     Double_t x    = xpos+step;
 771     Double_t bz=GetBz();   // getting the local Bz
 772     if (!seed->PropagateToX(x,bz))  return kFALSE;
 773     double ds = 0;
 774     if (matCorr || updTime) {
 775       xyz0[0]=xyz1[0]; // global pos at the beginning of step
 776       xyz0[1]=xyz1[1];
 777       xyz0[2]=xyz1[2];
 778       seed->GetXYZ(xyz1);    //  // global pos at the end of step
 779       if (matCorr) {
 780         MeanMaterialBudget(xyz0,xyz1,param);
 781         Double_t xrho=param[0]*param[4], xx0=param[1];
 782         if (dir>0) xrho = -xrho; // outward should be negative
 783         if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) return kFALSE;
 784         ds = param[4];
 785       }
 786        else { // matCorr is not requested but time integral is
 787         double d0 = xyz1[0]-xyz0[0];
 788         double d1 = xyz1[1]-xyz0[1];
 789         double d2 = xyz1[2]-xyz0[2];
 790         ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
 791       }
 792     }
 793     if (updTime) seed->AddTimeStep(ds);
 794     xpos = seed->GetX();
 795   }
 796   return kTRUE;
 797 }
 798
 799 //______________________________________________________________________________
 800 Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
 801 {
 802   // propagate seed to given x applying material correction if requested
 803   const Double_t kEpsilon = 1e-5;
 804   Double_t xpos     = seed->GetX();
 805   Int_t dir         = (xpos<xToGo) ? 1:-1;
 806   Double_t xyz0[3],xyz1[3],param[7];
 807   //
 808   Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
 809   if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
 810   while ( (xToGo-xpos)*dir > kEpsilon){
 811     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 812     Double_t x    = xpos+step;
 813     Double_t bz=GetBz();   // getting the local Bz
 814     if (!seed->PropagateTo(x,bz))  return kFALSE;
 815     double ds = 0;
 816     if (matCorr || updTime) {
 817       xyz0[0]=xyz1[0]; // global pos at the beginning of step
 818       xyz0[1]=xyz1[1];
 819       xyz0[2]=xyz1[2];
 820       seed->GetXYZ(xyz1);    //  // global pos at the end of step
 821       //
 822       if (matCorr) {
 823         MeanMaterialBudget(xyz0,xyz1,param);
 824         Double_t xrho=param[0]*param[4], xx0=param[1];
 825         if (dir>0) xrho = -xrho; // outward should be negative
 826         if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 827         ds = param[4];
 828       }
 829       else { // matCorr is not requested but time integral is
 830         double d0 = xyz1[0]-xyz0[0];
 831         double d1 = xyz1[1]-xyz0[1];
 832         double d2 = xyz1[2]-xyz0[2];
 833         ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
 834       }
 835     }
 836     if (updTime) seed->AddTimeStep(ds);
 837     //
 838     xpos = seed->GetX();
 839   }
 840   return kTRUE;
 841 }
 842
 843 //______________________________________________________________________________
 844 void AliITSUTrackerGlo::SaveCurrentTrackHypotheses()
 845 {
 846   // RS: shall we clean up killed seeds?
 847   fCurrHyp = 0;
 848   // TODO
 849
 850 }
 851
 852 //______________________________________________________________________________
 853 void AliITSUTrackerGlo::FinalizeHypotheses()
 854 {
 855   // select winner for each hypothesis, remove cl. sharing conflicts
 856   AliInfo("TODO");
 857   //
 858   int nh = fHypStore.GetEntriesFast();
 859   for (int ih=0;ih<nh;ih++) {
 860     AliITSUTrackHyp* hyp = (AliITSUTrackHyp*) fHypStore.UncheckedAt(ih);
 861     if (!hyp || !hyp->DefineWinner()) continue; // TODO
 862     CookMCLabel(hyp);
 863     UpdateESDTrack(hyp,AliESDtrack::kITSin);
 864   }
 865
 866 }
 867
 868 //______________________________________________________________________________
 869 void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp,Int_t flag)
 870 {
 871   // update ESD track with current best hypothesis
 872   AliESDtrack* esdTr = hyp->GetESDTrack();
 873   if (!esdTr) return;
 874   AliITSUSeed* win = hyp->GetWinner();
 875   if (!win) return;
 876   switch (flag) {
 877   case AliESDtrack::kITSin:
 878     esdTr->UpdateTrackParams(hyp,flag); // update kinematics
 879     // TODO: set cluster info
 880     break;
 881     //
 882   case AliESDtrack::kITSout:
 883     esdTr->UpdateTrackParams(hyp,flag); // update kinematics
 884     // TODO: avoid setting friend
 885     break;
 886     //
 887   case AliESDtrack::kITSrefit:
 888     esdTr->UpdateTrackParams(hyp,flag); // update kinematics
 889     // TODO: avoid setting cluster info
 890     break;
 891   default:
 892     AliFatal(Form("Unknown flag %d",flag));
 893   }
 894   //
 895   // transfer module indices
 896   // TODO
 897 }
 898
 899 //______________________________________________________________________________
 900 Bool_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest, Bool_t stopAtLastCl)
 901 {
 902   // refit track till radius rDest
 903   AliITSUTrackHyp tmpTr;
 904   //
 905   double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
 906   int dir,lrStart,lrStop;
 907   //
 908   dir = rCurr<rDest ? 1 : -1;
 909   lrStart = fITS->FindFirstLayerID(rCurr,dir);
 910   lrStop  = fITS->FindLastLayerID(rDest,dir);
 911   printf("Refit %d: Lr:%d (%f) -> Lr:%d (%f)\n",dir,lrStart,rCurr, lrStop,rDest);
 912   printf("Before refit: "); trc->AliExternalTrackParam::Print();
 913   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));
 914   //
 915   int ncl = trc->FetchClusterInfo(fClInfo);
 916   fCurrMass = trc->GetMass();
 917   tmpTr.AliKalmanTrack::operator=(*trc);
 918   tmpTr.SetChi2(0);
 919   int iclLr[2],nclLr,clCount=0;
 920   //
 921   for (int ilr=lrStart;ilr!=lrStop;ilr+=dir) {
 922     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
 923     if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
 924     int ilrA2,ilrA = lr->GetActiveID();
 925     // passive layer or active w/o hits will be traversed on the way to next cluster
 926     if (!lr->IsActive() || fClInfo[ilrA2=(ilrA<<1)]<0) continue;
 927     //
 928     if (ilr!=lrStart && !TransportToLayer(&tmpTr,lrStart,ilr)) {
 929       AliInfo(Form("Failed to transport %d -> %d\n",lrStart,ilr));
 930       return kFALSE; // go to the entrance to the layer
 931     }
 932     lrStart = ilr;
 933     //
 934     // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
 935     nclLr=0;
 936     if (dir>0) { // clusters are stored in increasing radius order
 937       iclLr[nclLr++]=fClInfo[ilrA2++];
 938       if (fClInfo[ilrA2]>=0) iclLr[nclLr++]=fClInfo[ilrA2];
 939     }
 940     else {
 941       if ( fClInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=fClInfo[ilrA2+1];
 942       iclLr[nclLr++]=fClInfo[ilrA2];
 943     }
 944     //
 945     for (int icl=0;icl<nclLr;icl++) {
 946       AliITSUClusterPix* clus =  (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
 947       AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
 948       if (!tmpTr.Rotate(sens->GetPhiTF())) {
 949         AliInfo(Form("Failed on rotate to %f",sens->GetPhiTF()));
 950         return kFALSE;
 951       }
 952       //printf("Refit cl:%d on lr %d Need to go %.4f -> %.4f\n",icl,ilrA,tmpTr.GetX(),sens->GetXTF()+clus->GetX());
 953       if (!PropagateSeed(&tmpTr,sens->GetXTF()+clus->GetX(),fCurrMass)) {
 954         AliInfo(Form("Failed on propagate to %f",sens->GetXTF()+clus->GetX()));
 955         return kFALSE;
 956       }
 957       if (!tmpTr.Update(clus)) {
 958         AliInfo(Form("Failed on Update"));
 959         return kFALSE;
 960       }
 961       //printf("AfterRefit: "); tmpTr.AliExternalTrackParam::Print();
 962       if (stopAtLastCl && ++clCount==ncl) return kTRUE; // it was requested to not propagate after last update
 963     }
 964     //
 965   }
 966   // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
 967   // Still, try to go as close as possible to rDest.
 968   //
 969   if (lrStart!=lrStop) {
 970     //printf("Going to last layer %d -> %d\n",lrStart,lrStop);
 971     if (!TransportToLayer(&tmpTr,lrStart,lrStop)) {
 972       AliInfo(Form("Failed on TransportToLayer %d->%d",lrStart,lrStop));
 973       return kTRUE;
 974     }
 975     if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) {
 976       AliFatal(Form("Failed on GoToExitFromLayer %d",lrStop));
 977       return kTRUE; // go till the exit from layer
 978     }
 979     //
 980     //printf("On exit from last layer\n");
 981     tmpTr.AliExternalTrackParam::Print();
 982     // go to the destination radius
 983     if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),dir)) return kTRUE;
 984     if (!PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) return kTRUE;
 985   }
 986   trc->AliKalmanTrack::operator=(tmpTr);
 987   printf("After refit (now at lr %d): ",lrStart); trc->AliExternalTrackParam::Print();
 988   return kTRUE;
 989 }
 990
 991 //__________________________________________________________________
 992 void AliITSUTrackerGlo::CookMCLabel(AliITSUTrackHyp* hyp)
 993 {
 994   // build MC label
 995   //
 996   const int kMaxLbPerCl = 3;
 997   int lbID[kMaxLayers*6],lbStat[kMaxLayers*6];
 998   Int_t lr,nLab=0,nCl=0;
 999   AliITSUSeed *seed = hyp->GetWinner();
1000   while(seed) {
1001     int clID = seed->GetLrCluster(lr);
1002     if (clID>=0) {
1003       AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1004       nCl++;
1005       for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
1006         int trLb = cl->GetLabel(imc);
1007         if (trLb<0) break;
1008         // search this mc track in already accounted ones
1009         int iLab;
1010         for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
1011         if (iLab<nLab) lbStat[iLab]++;
1012         else {
1013           lbID[nLab] = trLb;
1014           lbStat[nLab++] = 1;
1015         }
1016       } // loop over given cluster's labels
1017     }
1018     seed = (AliITSUSeed*)seed->GetParent();
1019   } // loop over clusters
1020   //
1021   if (nCl && nLab) {
1022     int maxLab=0,nTPCok=0;
1023     AliESDtrack* esdTr = hyp->GetESDTrack();
1024     int tpcLab = esdTr ? Abs(esdTr->GetTPCLabel()) : -kDummyLabel;
1025     for (int ilb=nLab;ilb--;) {
1026       int st = lbStat[ilb];
1027       if (lbStat[maxLab]<st) maxLab=ilb;
1028       if (tpcLab==lbID[ilb]) nTPCok += st;
1029     }
1030     hyp->SetFakeRatio(1.-float(nTPCok)/nCl);
1031     hyp->SetLabel( nTPCok==nCl ? tpcLab : -tpcLab);
1032     hyp->SetITSLabel( lbStat[maxLab]==nCl ? lbStat[maxLab] : -lbStat[maxLab]); // winner label
1033     return;
1034   }
1035   //
1036   hyp->SetFakeRatio(-1.);
1037   hyp->SetLabel( kDummyLabel );
1038   hyp->SetITSLabel( kDummyLabel );
1039 }