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