ITS/AliITStrackerMI.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 /* $Id$ */
  16
  17 //-------------------------------------------------------------------------
  18 //               Implementation of the ITS tracker class
  19 //    It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
  20 //                   and fills with them the ESD
  21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  22 //          Current support and development:
  23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
  24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
  26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
  27 //-------------------------------------------------------------------------
  28
  29 #include <TMatrixD.h>
  30 #include <TTree.h>
  31 #include <TDatabasePDG.h>
  32 #include <TString.h>
  33 #include <TRandom.h>
  34 #include <TTreeStream.h>
  35
  36
  37 #include "AliLog.h"
  38 #include "AliITSPlaneEff.h"
  39 #include "AliITSCalibrationSPD.h"
  40 #include "AliITSCalibrationSDD.h"
  41 #include "AliITSCalibrationSSD.h"
  42 #include "AliCDBEntry.h"
  43 #include "AliCDBManager.h"
  44 #include "AliAlignObj.h"
  45 #include "AliTrackPointArray.h"
  46 #include "AliESDVertex.h"
  47 #include "AliESDEvent.h"
  48 #include "AliESDtrack.h"
  49 #include "AliV0.h"
  50 #include "AliHelix.h"
  51 #include "AliITSChannelStatus.h"
  52 #include "AliITSDetTypeRec.h"
  53 #include "AliITSRecPoint.h"
  54 #include "AliITSgeomTGeo.h"
  55 #include "AliITSReconstructor.h"
  56 #include "AliITSClusterParam.h"
  57 #include "AliITSsegmentation.h"
  58 #include "AliITSCalibration.h"
  59 #include "AliITSPlaneEffSPD.h"
  60 #include "AliITSPlaneEffSDD.h"
  61 #include "AliITSPlaneEffSSD.h"
  62 #include "AliITSV0Finder.h"
  63 #include "AliITStrackerMI.h"
  64
  65 ClassImp(AliITStrackerMI)
  66
  67 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
  68
  69 AliITStrackerMI::AliITStrackerMI():AliTracker(),
  70 fI(0),
  71 fBestTrack(),
  72 fTrackToFollow(),
  73 fTrackHypothesys(),
  74 fBestHypothesys(),
  75 fOriginal(),
  76 fCurrentEsdTrack(),
  77 fPass(0),
  78 fAfterV0(kFALSE),
  79 fLastLayerToTrackTo(0),
  80 fCoefficients(0),
  81 fEsd(0),
  82 fTrackingPhase("Default"),
  83 fUseTGeo(3),
  84 fNtracks(0),
  85 fxOverX0Pipe(-1.),
  86 fxTimesRhoPipe(-1.),
  87 fxOverX0PipeTrks(0),
  88 fxTimesRhoPipeTrks(0),
  89 fxOverX0ShieldTrks(0),
  90 fxTimesRhoShieldTrks(0),
  91 fxOverX0LayerTrks(0),
  92 fxTimesRhoLayerTrks(0),
  93 fDebugStreamer(0),
  94 fITSChannelStatus(0),
  95 fkDetTypeRec(0),
  96 fPlaneEff(0) {
  97   //Default constructor
  98   Int_t i;
  99   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
 100   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 101   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 102 }
 103 //------------------------------------------------------------------------
 104 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
 105 fI(AliITSgeomTGeo::GetNLayers()),
 106 fBestTrack(),
 107 fTrackToFollow(),
 108 fTrackHypothesys(),
 109 fBestHypothesys(),
 110 fOriginal(),
 111 fCurrentEsdTrack(),
 112 fPass(0),
 113 fAfterV0(kFALSE),
 114 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
 115 fCoefficients(0),
 116 fEsd(0),
 117 fTrackingPhase("Default"),
 118 fUseTGeo(3),
 119 fNtracks(0),
 120 fxOverX0Pipe(-1.),
 121 fxTimesRhoPipe(-1.),
 122 fxOverX0PipeTrks(0),
 123 fxTimesRhoPipeTrks(0),
 124 fxOverX0ShieldTrks(0),
 125 fxTimesRhoShieldTrks(0),
 126 fxOverX0LayerTrks(0),
 127 fxTimesRhoLayerTrks(0),
 128 fDebugStreamer(0),
 129 fITSChannelStatus(0),
 130 fkDetTypeRec(0),
 131 fPlaneEff(0) {
 132   //--------------------------------------------------------------------
 133   //This is the AliITStrackerMI constructor
 134   //--------------------------------------------------------------------
 135   if (geom) {
 136     AliWarning("\"geom\" is actually a dummy argument !");
 137   }
 138
 139   fCoefficients = 0;
 140   fAfterV0     = kFALSE;
 141
 142   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 143     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 144     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 145
 146     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 147     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 148     Double_t poff=TMath::ATan2(y,x);
 149     Double_t zoff=z;
 150     Double_t r=TMath::Sqrt(x*x + y*y);
 151
 152     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 153     r += TMath::Sqrt(x*x + y*y);
 154     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 155     r += TMath::Sqrt(x*x + y*y);
 156     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 157     r += TMath::Sqrt(x*x + y*y);
 158     r*=0.25;
 159
 160     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
 161
 162     for (Int_t j=1; j<nlad+1; j++) {
 163       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 164         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 165         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 166         m.Multiply(tm);
 167         Double_t txyz[3]={0.};
 168         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 169         m.LocalToMaster(txyz,xyz);
 170         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 171         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 172
 173         if (phi<0) phi+=TMath::TwoPi();
 174         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 175
 176         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 177         new(&det) AliITSdetector(r,phi);
 178         // compute the real radius (with misalignment)
 179         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 180         mmisal.Multiply(tm);
 181         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 182         mmisal.LocalToMaster(txyz,xyz);
 183         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 184         det.SetRmisal(rmisal);
 185
 186       } // end loop on detectors
 187     } // end loop on ladders
 188   } // end loop on layers
 189
 190
 191   fI=AliITSgeomTGeo::GetNLayers();
 192
 193   fPass=0;
 194   fConstraint[0]=1; fConstraint[1]=0;
 195
 196   Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
 197                      AliITSReconstructor::GetRecoParam()->GetYVdef(),
 198                      AliITSReconstructor::GetRecoParam()->GetZVdef()};
 199   Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
 200                      AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
 201                      AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
 202   SetVertex(xyzVtx,ersVtx);
 203
 204   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=AliITSRecoParam::GetLayersNotToSkip(i);
 205   fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
 206   for (Int_t i=0;i<100000;i++){
 207     fBestTrackIndex[i]=0;
 208   }
 209
 210   // store positions of centre of SPD modules (in z)
 211   Double_t tr[3];
 212   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 213   fSPDdetzcentre[0] = tr[2];
 214   AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 215   fSPDdetzcentre[1] = tr[2];
 216   AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 217   fSPDdetzcentre[2] = tr[2];
 218   AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 219   fSPDdetzcentre[3] = tr[2];
 220
 221   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
 222   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 223     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 224     fUseTGeo = 3;
 225   }
 226
 227   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 228   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 229
 230   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 231
 232   // only for plane efficiency evaluation
 233   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 234       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
 235     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
 236     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane))
 237       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
 238     if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
 239     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
 240     else fPlaneEff = new AliITSPlaneEffSSD();
 241     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
 242        if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
 243     if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
 244   }
 245 }
 246 //------------------------------------------------------------------------
 247 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
 248 fI(tracker.fI),
 249 fBestTrack(tracker.fBestTrack),
 250 fTrackToFollow(tracker.fTrackToFollow),
 251 fTrackHypothesys(tracker.fTrackHypothesys),
 252 fBestHypothesys(tracker.fBestHypothesys),
 253 fOriginal(tracker.fOriginal),
 254 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
 255 fPass(tracker.fPass),
 256 fAfterV0(tracker.fAfterV0),
 257 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
 258 fCoefficients(tracker.fCoefficients),
 259 fEsd(tracker.fEsd),
 260 fTrackingPhase(tracker.fTrackingPhase),
 261 fUseTGeo(tracker.fUseTGeo),
 262 fNtracks(tracker.fNtracks),
 263 fxOverX0Pipe(tracker.fxOverX0Pipe),
 264 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 265 fxOverX0PipeTrks(0),
 266 fxTimesRhoPipeTrks(0),
 267 fxOverX0ShieldTrks(0),
 268 fxTimesRhoShieldTrks(0),
 269 fxOverX0LayerTrks(0),
 270 fxTimesRhoLayerTrks(0),
 271 fDebugStreamer(tracker.fDebugStreamer),
 272 fITSChannelStatus(tracker.fITSChannelStatus),
 273 fkDetTypeRec(tracker.fkDetTypeRec),
 274 fPlaneEff(tracker.fPlaneEff) {
 275   //Copy constructor
 276   Int_t i;
 277   for(i=0;i<4;i++) {
 278     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 279   }
 280   for(i=0;i<6;i++) {
 281     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 282     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 283   }
 284   for(i=0;i<2;i++) {
 285     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 286     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 287   }
 288 }
 289 //------------------------------------------------------------------------
 290 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
 291   //Assignment operator
 292   this->~AliITStrackerMI();
 293   new(this) AliITStrackerMI(tracker);
 294   return *this;
 295 }
 296 //------------------------------------------------------------------------
 297 AliITStrackerMI::~AliITStrackerMI()
 298 {
 299   //
 300   //destructor
 301   //
 302   if (fCoefficients) delete [] fCoefficients;
 303   DeleteTrksMaterialLUT();
 304   if (fDebugStreamer) {
 305     //fDebugStreamer->Close();
 306     delete fDebugStreamer;
 307   }
 308   if(fITSChannelStatus) delete fITSChannelStatus;
 309   if(fPlaneEff) delete fPlaneEff;
 310 }
 311 //------------------------------------------------------------------------
 312 void AliITStrackerMI::SetLayersNotToSkip(const Int_t *l) {
 313   //--------------------------------------------------------------------
 314   //This function set masks of the layers which must be not skipped
 315   //--------------------------------------------------------------------
 316   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=l[i];
 317 }
 318 //------------------------------------------------------------------------
 319 void AliITStrackerMI::ReadBadFromDetTypeRec() {
 320   //--------------------------------------------------------------------
 321   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
 322   //i.e. from OCDB
 323   //--------------------------------------------------------------------
 324
 325   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
 326
 327   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
 328
 329   if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
 330
 331   // ITS channels map
 332   if(fITSChannelStatus) delete fITSChannelStatus;
 333   fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
 334
 335   // ITS detectors and chips
 336   Int_t i=0,j=0,k=0,ndet=0;
 337   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 338     Int_t nBadDetsPerLayer=0;
 339     ndet=AliITSgeomTGeo::GetNDetectors(i);
 340     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
 341       for (k=1; k<ndet+1; k++) {
 342         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 343         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
 344         if(det.IsBad()) {nBadDetsPerLayer++;}
 345       } // end loop on detectors
 346     } // end loop on ladders
 347     Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
 348   } // end loop on layers
 349
 350   return;
 351 }
 352 //------------------------------------------------------------------------
 353 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 354   //--------------------------------------------------------------------
 355   //This function loads ITS clusters
 356   //--------------------------------------------------------------------
 357   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 358   if (!branch) {
 359     Error("LoadClusters"," can't get the branch !\n");
 360     return 1;
 361   }
 362
 363   static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 364   branch->SetAddress(&clusters);
 365
 366   Int_t i=0,j=0,ndet=0;
 367   Int_t detector=0;
 368   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 369     ndet=fgLayers[i].GetNdetectors();
 370     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 371     for (; j<jmax; j++) {
 372       if (!cTree->GetEvent(j)) continue;
 373       Int_t ncl=clusters->GetEntriesFast();
 374       SignDeltas(clusters,GetZ());
 375
 376       while (ncl--) {
 377         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 378         detector=c->GetDetectorIndex();
 379
 380         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 381
 382         fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 383       }
 384       clusters->Delete();
 385       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 386       // zwindow cm from the dead zone
 387       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
 388         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
 389           Int_t lab[4]   = {0,0,0,detector};
 390           Int_t info[3]  = {0,0,i};
 391           Float_t q      = 0.; // this identifies virtual clusters
 392           Float_t hit[5] = {xdead,
 393                             0.,
 394                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
 395                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
 396                             q};
 397           Bool_t local   = kTRUE;
 398           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
 399           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 400           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 401             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 402           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 403           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 404             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 405           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 406           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 407             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 408           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 409           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 410             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 411           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 412           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 413             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 414           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 415           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 416             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 417         }
 418       } // "virtual" clusters in SPD
 419
 420     }
 421     //
 422     fgLayers[i].ResetRoad(); //road defined by the cluster density
 423     fgLayers[i].SortClusters();
 424   }
 425
 426   dummy.Clear();
 427
 428   return 0;
 429 }
 430 //------------------------------------------------------------------------
 431 void AliITStrackerMI::UnloadClusters() {
 432   //--------------------------------------------------------------------
 433   //This function unloads ITS clusters
 434   //--------------------------------------------------------------------
 435   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
 436 }
 437 //------------------------------------------------------------------------
 438 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
 439   //--------------------------------------------------------------------
 440   // Publishes all pointers to clusters known to the tracker into the
 441   // passed object array.
 442   // The ownership is not transfered - the caller is not expected to delete
 443   // the clusters.
 444   //--------------------------------------------------------------------
 445
 446   for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 447     for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
 448       AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
 449       array->AddLast(cl);
 450     }
 451   }
 452
 453   return;
 454 }
 455 //------------------------------------------------------------------------
 456 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
 457   //--------------------------------------------------------------------
 458   // Correction for the material between the TPC and the ITS
 459   //--------------------------------------------------------------------
 460   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 461       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 462       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 463       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 464   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 465       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 466       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 467       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 468   } else {
 469     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
 470     return 0;
 471   }
 472
 473   return 1;
 474 }
 475 //------------------------------------------------------------------------
 476 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
 477   //--------------------------------------------------------------------
 478   // This functions reconstructs ITS tracks
 479   // The clusters must be already loaded !
 480   //--------------------------------------------------------------------
 481
 482
 483   fTrackingPhase="Clusters2Tracks";
 484
 485   TObjArray itsTracks(15000);
 486   fOriginal.Clear();
 487   fEsd = event;         // store pointer to the esd
 488
 489   // temporary (for cosmics)
 490   if(event->GetVertex()) {
 491     TString title = event->GetVertex()->GetTitle();
 492     if(title.Contains("cosmics")) {
 493       Double_t xyz[3]={GetX(),GetY(),GetZ()};
 494       Double_t exyz[3]={0.1,0.1,0.1};
 495       SetVertex(xyz,exyz);
 496     }
 497   }
 498   // temporary
 499
 500   {/* Read ESD tracks */
 501     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 502     Int_t nentr=event->GetNumberOfTracks();
 503     Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 504     while (nentr--) {
 505       AliESDtrack *esd=event->GetTrack(nentr);
 506       //  ---- for debugging:
 507       //if(TMath::Abs(esd->GetX()-83.65)<0.1) { FILE *f=fopen("tpc.dat","a"); fprintf(f,"%f %f %f %f %f %f\n",(Float_t)event->GetEventNumberInFile(),(Float_t)TMath::Abs(esd->GetLabel()),(Float_t)esd->GetX(),(Float_t)esd->GetY(),(Float_t)esd->GetZ(),(Float_t)esd->Pt()); fclose(f); }
 508
 509       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 510       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 511       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 512       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 513       AliITStrackMI *t=0;
 514       try {
 515         t=new AliITStrackMI(*esd);
 516       } catch (const Char_t *msg) {
 517         //Warning("Clusters2Tracks",msg);
 518         delete t;
 519         continue;
 520       }
 521       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
 522       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 523
 524
 525       // look at the ESD mass hypothesys !
 526       if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
 527       // write expected q
 528       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 529
 530       if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
 531         //track - can be  V0 according to TPC
 532       } else {
 533         if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
 534           delete t;
 535           continue;
 536         }
 537         if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
 538           delete t;
 539           continue;
 540         }
 541         if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
 542           delete t;
 543           continue;
 544         }
 545         if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 546           delete t;
 547           continue;
 548         }
 549       }
 550       t->SetReconstructed(kFALSE);
 551       itsTracks.AddLast(t);
 552       fOriginal.AddLast(t);
 553     }
 554   } /* End Read ESD tracks */
 555
 556   itsTracks.Sort();
 557   fOriginal.Sort();
 558   Int_t nentr=itsTracks.GetEntriesFast();
 559   fTrackHypothesys.Expand(nentr);
 560   fBestHypothesys.Expand(nentr);
 561   MakeCoefficients(nentr);
 562   if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
 563   Int_t ntrk=0;
 564   // THE TWO TRACKING PASSES
 565   for (fPass=0; fPass<2; fPass++) {
 566      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 567      for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
 568        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
 569        if (t==0) continue;              //this track has been already tracked
 570        //cout<<"========== "<<fPass<<"    "<<fCurrentEsdTrack<<" =========\n";
 571        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 572        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
 573        if (fConstraint[fPass]) {
 574          if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
 575              TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
 576        }
 577
 578        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 579        AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
 580        fI = 6;
 581        ResetTrackToFollow(*t);
 582        ResetBestTrack();
 583
 584        FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
 585
 586
 587        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 588        //
 589        AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 590        if (!besttrack) continue;
 591        besttrack->SetLabel(tpcLabel);
 592        //       besttrack->CookdEdx();
 593        CookdEdx(besttrack);
 594        besttrack->SetFakeRatio(1.);
 595        CookLabel(besttrack,0.); //For comparison only
 596        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 597
 598        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 599
 600        t->SetReconstructed(kTRUE);
 601        ntrk++;
 602        AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
 603      }
 604      GetBestHypothesysMIP(itsTracks);
 605   } // end loop on the two tracking passes
 606
 607   if(event->GetNumberOfV0s()>0) UpdateTPCV0(event);
 608   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) FindV02(event);
 609   //if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
 610   //if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
 611   fAfterV0 = kTRUE;
 612   //
 613   itsTracks.Delete();
 614   //
 615   Int_t entries = fTrackHypothesys.GetEntriesFast();
 616   for (Int_t ientry=0; ientry<entries; ientry++) {
 617     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
 618     if (array) array->Delete();
 619     delete fTrackHypothesys.RemoveAt(ientry);
 620   }
 621
 622   fTrackHypothesys.Delete();
 623   fBestHypothesys.Delete();
 624   fOriginal.Clear();
 625   delete [] fCoefficients;
 626   fCoefficients=0;
 627   DeleteTrksMaterialLUT();
 628
 629   Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
 630
 631   fTrackingPhase="Default";
 632
 633   return 0;
 634 }
 635 //------------------------------------------------------------------------
 636 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
 637   //--------------------------------------------------------------------
 638   // This functions propagates reconstructed ITS tracks back
 639   // The clusters must be loaded !
 640   //--------------------------------------------------------------------
 641   fTrackingPhase="PropagateBack";
 642   Int_t nentr=event->GetNumberOfTracks();
 643   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 644
 645   Int_t ntrk=0;
 646   for (Int_t i=0; i<nentr; i++) {
 647      AliESDtrack *esd=event->GetTrack(i);
 648
 649      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 650      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 651
 652      AliITStrackMI *t=0;
 653      try {
 654         t=new AliITStrackMI(*esd);
 655      } catch (const Char_t *msg) {
 656        //Warning("PropagateBack",msg);
 657         delete t;
 658         continue;
 659      }
 660      t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 661
 662      ResetTrackToFollow(*t);
 663
 664      // propagate to vertex [SR, GSI 17.02.2003]
 665      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 666      if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 667        if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
 668          fTrackToFollow.StartTimeIntegral();
 669        // from vertex to outside pipe
 670        CorrectForPipeMaterial(&fTrackToFollow,"outward");
 671      }
 672
 673      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
 674      if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
 675         if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
 676           delete t;
 677           continue;
 678         }
 679         fTrackToFollow.SetLabel(t->GetLabel());
 680         //fTrackToFollow.CookdEdx();
 681         CookLabel(&fTrackToFollow,0.); //For comparison only
 682         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 683         //UseClusters(&fTrackToFollow);
 684         ntrk++;
 685      }
 686      delete t;
 687   }
 688
 689   Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
 690
 691   fTrackingPhase="Default";
 692
 693   return 0;
 694 }
 695 //------------------------------------------------------------------------
 696 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
 697   //--------------------------------------------------------------------
 698   // This functions refits ITS tracks using the
 699   // "inward propagated" TPC tracks
 700   // The clusters must be loaded !
 701   //--------------------------------------------------------------------
 702   fTrackingPhase="RefitInward";
 703
 704   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) RefitV02(event);
 705   //if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
 706
 707   Int_t nentr=event->GetNumberOfTracks();
 708   Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 709
 710   Int_t ntrk=0;
 711   for (Int_t i=0; i<nentr; i++) {
 712     AliESDtrack *esd=event->GetTrack(i);
 713
 714     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 715     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 716     if (esd->GetStatus()&AliESDtrack::kTPCout)
 717       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 718
 719     AliITStrackMI *t=0;
 720     try {
 721         t=new AliITStrackMI(*esd);
 722     } catch (const Char_t *msg) {
 723       //Warning("RefitInward",msg);
 724         delete t;
 725         continue;
 726     }
 727     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 728     if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 729        delete t;
 730        continue;
 731     }
 732
 733     ResetTrackToFollow(*t);
 734     fTrackToFollow.ResetClusters();
 735
 736     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
 737       fTrackToFollow.ResetCovariance(10.);
 738
 739     //Refitting...
 740     Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 741                AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
 742     AliDebug(2,Form("Refit LABEL %d  %d",t->GetLabel(),t->GetNumberOfClusters()));
 743     if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
 744        AliDebug(2,"  refit OK");
 745        fTrackToFollow.SetLabel(t->GetLabel());
 746        //       fTrackToFollow.CookdEdx();
 747        CookdEdx(&fTrackToFollow);
 748
 749        CookLabel(&fTrackToFollow,0.0); //For comparison only
 750
 751        //The beam pipe
 752        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 753          fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
 754          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 755          //printf("                                       %d\n",esdTrack->GetITSModuleIndex(0));
 756          //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
 757          Double_t r[3]={0.,0.,0.};
 758          Double_t maxD=3.;
 759          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 760          ntrk++;
 761        }
 762     }
 763     delete t;
 764   }
 765
 766   Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
 767
 768   fTrackingPhase="Default";
 769
 770   return 0;
 771 }
 772 //------------------------------------------------------------------------
 773 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 774   //--------------------------------------------------------------------
 775   //       Return pointer to a given cluster
 776   //--------------------------------------------------------------------
 777   Int_t l=(index & 0xf0000000) >> 28;
 778   Int_t c=(index & 0x0fffffff) >> 00;
 779   return fgLayers[l].GetCluster(c);
 780 }
 781 //------------------------------------------------------------------------
 782 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 783   //--------------------------------------------------------------------
 784   // Get track space point with index i
 785   //--------------------------------------------------------------------
 786
 787   Int_t l=(index & 0xf0000000) >> 28;
 788   Int_t c=(index & 0x0fffffff) >> 00;
 789   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 790   Int_t idet = cl->GetDetectorIndex();
 791
 792   Float_t xyz[3];
 793   Float_t cov[6];
 794   cl->GetGlobalXYZ(xyz);
 795   cl->GetGlobalCov(cov);
 796   p.SetXYZ(xyz, cov);
 797   p.SetCharge(cl->GetQ());
 798   p.SetDriftTime(cl->GetDriftTime());
 799   p.SetChargeRatio(cl->GetChargeRatio());
 800   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 801   switch (l) {
 802   case 0:
 803     iLayer = AliGeomManager::kSPD1;
 804     break;
 805   case 1:
 806     iLayer = AliGeomManager::kSPD2;
 807     break;
 808   case 2:
 809     iLayer = AliGeomManager::kSDD1;
 810     break;
 811   case 3:
 812     iLayer = AliGeomManager::kSDD2;
 813     break;
 814   case 4:
 815     iLayer = AliGeomManager::kSSD1;
 816     break;
 817   case 5:
 818     iLayer = AliGeomManager::kSSD2;
 819     break;
 820   default:
 821     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 822     break;
 823   };
 824   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 825   p.SetVolumeID((UShort_t)volid);
 826   return kTRUE;
 827 }
 828 //------------------------------------------------------------------------
 829 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
 830                         AliTrackPoint& p, const AliESDtrack *t) {
 831   //--------------------------------------------------------------------
 832   // Get track space point with index i
 833   // (assign error estimated during the tracking)
 834   //--------------------------------------------------------------------
 835
 836   Int_t l=(index & 0xf0000000) >> 28;
 837   Int_t c=(index & 0x0fffffff) >> 00;
 838   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 839   Int_t idet = cl->GetDetectorIndex();
 840
 841   const AliITSdetector &det=fgLayers[l].GetDetector(idet);
 842
 843   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 844   Float_t detxy[2];
 845   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 846   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 847   Double_t alpha = t->GetAlpha();
 848   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 849   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
 850   phi += alpha-det.GetPhi();
 851   Float_t tgphi = TMath::Tan(phi);
 852
 853   Float_t tgl = t->GetTgl(); // tgl about const along track
 854   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 855
 856   Float_t errlocalx,errlocalz;
 857   Bool_t addMisalErr=kFALSE;
 858   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
 859
 860   Float_t xyz[3];
 861   Float_t cov[6];
 862   cl->GetGlobalXYZ(xyz);
 863   //  cl->GetGlobalCov(cov);
 864   Float_t pos[3] = {0.,0.,0.};
 865   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
 866   tmpcl.GetGlobalCov(cov);
 867
 868   p.SetXYZ(xyz, cov);
 869   p.SetCharge(cl->GetQ());
 870   p.SetDriftTime(cl->GetDriftTime());
 871   p.SetChargeRatio(cl->GetChargeRatio());
 872
 873   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 874   switch (l) {
 875   case 0:
 876     iLayer = AliGeomManager::kSPD1;
 877     break;
 878   case 1:
 879     iLayer = AliGeomManager::kSPD2;
 880     break;
 881   case 2:
 882     iLayer = AliGeomManager::kSDD1;
 883     break;
 884   case 3:
 885     iLayer = AliGeomManager::kSDD2;
 886     break;
 887   case 4:
 888     iLayer = AliGeomManager::kSSD1;
 889     break;
 890   case 5:
 891     iLayer = AliGeomManager::kSSD2;
 892     break;
 893   default:
 894     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 895     break;
 896   };
 897   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 898
 899   p.SetVolumeID((UShort_t)volid);
 900   return kTRUE;
 901 }
 902 //------------------------------------------------------------------------
 903 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
 904 {
 905   //--------------------------------------------------------------------
 906   // Follow prolongation tree
 907   //--------------------------------------------------------------------
 908   //
 909   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
 910   Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 911
 912
 913   AliESDtrack * esd = otrack->GetESDtrack();
 914   if (esd->GetV0Index(0)>0) {
 915     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
 916     //                      mapping of ESD track is different as ITS track in Containers
 917     //                      Need something more stable
 918     //                      Indexes are set back again to the ESD track indexes in UpdateTPCV0
 919     for (Int_t i=0;i<3;i++){
 920       Int_t  index = esd->GetV0Index(i);
 921       if (index==0) break;
 922       AliESDv0 * vertex = fEsd->GetV0(index);
 923       if (vertex->GetStatus()<0) continue;     // rejected V0
 924       //
 925       if (esd->GetSign()>0) {
 926         vertex->SetIndex(0,esdindex);
 927       } else {
 928         vertex->SetIndex(1,esdindex);
 929       }
 930     }
 931   }
 932   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
 933   if (!bestarray){
 934     bestarray = new TObjArray(5);
 935     fBestHypothesys.AddAt(bestarray,esdindex);
 936   }
 937
 938   //
 939   //setup tree of the prolongations
 940   //
 941   static AliITStrackMI tracks[7][100];
 942   AliITStrackMI *currenttrack;
 943   static AliITStrackMI currenttrack1;
 944   static AliITStrackMI currenttrack2;
 945   static AliITStrackMI backuptrack;
 946   Int_t ntracks[7];
 947   Int_t nindexes[7][100];
 948   Float_t normalizedchi2[100];
 949   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
 950   otrack->SetNSkipped(0);
 951   new (&(tracks[6][0])) AliITStrackMI(*otrack);
 952   ntracks[6]=1;
 953   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
 954   Int_t modstatus = 1; // found
 955   Float_t xloc,zloc;
 956   //
 957   //
 958   // follow prolongations
 959   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
 960     AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
 961     fI = ilayer;
 962     //
 963     AliITSlayer &layer=fgLayers[ilayer];
 964     Double_t r = layer.GetR();
 965     ntracks[ilayer]=0;
 966     //
 967     //
 968     Int_t nskipped=0;
 969     Float_t nused =0;
 970     for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
 971       //set current track
 972       if (ntracks[ilayer]>=100) break;
 973       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
 974       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
 975       if (ntracks[ilayer]>15+ilayer){
 976         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
 977         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
 978       }
 979
 980       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
 981
 982       // material between SSD and SDD, SDD and SPD
 983       if (ilayer==3)
 984         if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
 985       if (ilayer==1)
 986         if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
 987
 988       // detector number
 989       Double_t phi,z;
 990       if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
 991       Int_t idet=layer.FindDetectorIndex(phi,z);
 992
 993       Double_t trackGlobXYZ1[3];
 994       if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
 995
 996       // Get the budget to the primary vertex for the current track being prolonged
 997       Double_t budgetToPrimVertex = GetEffectiveThickness();
 998
 999       // check if we allow a prolongation without point
1000       Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
1001       if (skip) {
1002         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1003         // propagate to the layer radius
1004         Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1005         if(!vtrack->Propagate(xToGo)) continue;
1006         // apply correction for material of the current layer
1007         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1008         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1009         vtrack->SetClIndex(ilayer,-1);
1010         modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1011         if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1012           vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1013         }
1014         if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1015         ntracks[ilayer]++;
1016         continue;
1017       }
1018
1019       // track outside layer acceptance in z
1020       if (idet<0) continue;
1021
1022       //propagate to the intersection with the detector plane
1023       const AliITSdetector &det=layer.GetDetector(idet);
1024       new(&currenttrack2)  AliITStrackMI(currenttrack1);
1025       if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1026       if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1027       currenttrack1.SetDetectorIndex(idet);
1028       currenttrack2.SetDetectorIndex(idet);
1029       if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
1030
1031       //***************
1032       // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1033       //
1034       // road in global (rphi,z) [i.e. in tracking ref. system]
1035       Double_t zmin,zmax,ymin,ymax;
1036       if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1037
1038       // select clusters in road
1039       layer.SelectClusters(zmin,zmax,ymin,ymax);
1040       //********************
1041
1042       // Define criteria for track-cluster association
1043       Double_t msz = currenttrack1.GetSigmaZ2() +
1044         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1045         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1046         AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1047       Double_t msy = currenttrack1.GetSigmaY2() +
1048         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1049         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1050         AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1051       if (constrain) {
1052         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1053         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1054       }  else {
1055         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1056         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1057       }
1058       msz = 1./msz; // 1/RoadZ^2
1059       msy = 1./msy; // 1/RoadY^2
1060
1061       //
1062       //
1063       // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1064       //
1065       const AliITSRecPoint *cl=0;
1066       Int_t clidx=-1;
1067       Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1068       Bool_t deadzoneSPD=kFALSE;
1069       currenttrack = &currenttrack1;
1070
1071       // check if the road contains a dead zone
1072       Bool_t noClusters = kFALSE;
1073       if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1074       if (noClusters) AliDebug(2,"no clusters in road");
1075       Double_t dz=0.5*(zmax-zmin);
1076       Double_t dy=0.5*(ymax-ymin);
1077       Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters);
1078       if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1079       // create a prolongation without clusters (check also if there are no clusters in the road)
1080       if (dead ||
1081           (noClusters &&
1082            AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1083         AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1084         updatetrack->SetClIndex(ilayer,-1);
1085         if (dead==0) {
1086           modstatus = 5; // no cls in road
1087         } else if (dead==1) {
1088           modstatus = 7; // holes in z in SPD
1089         } else if (dead==2 || dead==3) {
1090           modstatus = 2; // dead from OCDB
1091         }
1092         updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1093         // apply correction for material of the current layer
1094         CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1095         if (constrain) { // apply vertex constrain
1096           updatetrack->SetConstrain(constrain);
1097           Bool_t isPrim = kTRUE;
1098           if (ilayer<4) { // check that it's close to the vertex
1099             updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1100             if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1101                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1102                 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1103                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1104           }
1105           if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1106         }
1107         if (dead) {
1108           updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1109           if (dead==1) { // dead zone at z=0,+-7cm in SPD
1110             updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1111             deadzoneSPD=kTRUE;
1112           }
1113         }
1114         ntracks[ilayer]++;
1115       }
1116
1117       clidx=-1;
1118       // loop over clusters in the road
1119       while ((cl=layer.GetNextCluster(clidx))!=0) {
1120         if (ntracks[ilayer]>95) break; //space for skipped clusters
1121         Bool_t changedet =kFALSE;
1122         if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1123         Int_t idetc=cl->GetDetectorIndex();
1124
1125         if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1126           // take into account misalignment (bring track to real detector plane)
1127           Double_t xTrOrig = currenttrack->GetX();
1128           if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1129           // a first cut on track-cluster distance
1130           if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1131                (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1132             {  // cluster not associated to track
1133               AliDebug(2,"not associated");
1134               continue;
1135             }
1136           // bring track back to ideal detector plane
1137           if (!currenttrack->Propagate(xTrOrig)) continue;
1138         } else {                                      // have to move track to cluster's detector
1139           const AliITSdetector &detc=layer.GetDetector(idetc);
1140           // a first cut on track-cluster distance
1141           Double_t y;
1142           if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1143           if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1144                (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1145             continue; // cluster not associated to track
1146           //
1147           new (&backuptrack) AliITStrackMI(currenttrack2);
1148           changedet = kTRUE;
1149           currenttrack =&currenttrack2;
1150           if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1151             new (currenttrack) AliITStrackMI(backuptrack);
1152             changedet = kFALSE;
1153             continue;
1154           }
1155           currenttrack->SetDetectorIndex(idetc);
1156           // Get again the budget to the primary vertex
1157           // for the current track being prolonged, if had to change detector
1158           //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1159         }
1160
1161         // calculate track-clusters chi2
1162         chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1163         // chi2 cut
1164         AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1165         if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1166           if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1167           if (ntracks[ilayer]>=100) continue;
1168           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1169           updatetrack->SetClIndex(ilayer,-1);
1170           if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
1171
1172           if (cl->GetQ()!=0) { // real cluster
1173             if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1174               AliDebug(2,"update failed");
1175               continue;
1176             }
1177             updatetrack->SetSampledEdx(cl->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
1178             modstatus = 1; // found
1179           } else {             // virtual cluster in dead zone
1180             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1181             updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1182             modstatus = 7; // holes in z in SPD
1183           }
1184
1185           if (changedet) {
1186             Float_t xlocnewdet,zlocnewdet;
1187             if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1188               updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1189             }
1190           } else {
1191             updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1192           }
1193           if (cl->IsUsed()) updatetrack->IncrementNUsed();
1194
1195           // apply correction for material of the current layer
1196           CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1197
1198           if (constrain) { // apply vertex constrain
1199             updatetrack->SetConstrain(constrain);
1200             Bool_t isPrim = kTRUE;
1201             if (ilayer<4) { // check that it's close to the vertex
1202               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1203               if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1204                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1205                   TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1206                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1207             }
1208             if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1209           } //apply vertex constrain
1210           ntracks[ilayer]++;
1211         }  // create new hypothesis
1212         else {
1213           AliDebug(2,"chi2 too large");
1214         }
1215
1216       } // loop over possible prolongations
1217
1218       // allow one prolongation without clusters
1219       if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1220         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1221         // apply correction for material of the current layer
1222         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1223         vtrack->SetClIndex(ilayer,-1);
1224         modstatus = 3; // skipped
1225         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1226         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1227         vtrack->IncrementNSkipped();
1228         ntracks[ilayer]++;
1229       }
1230
1231       // allow one prolongation without clusters for tracks with |tgl|>1.1
1232       if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) {  //big theta - for low flux
1233         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1234         // apply correction for material of the current layer
1235         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1236         vtrack->SetClIndex(ilayer,-1);
1237         modstatus = 3; // skipped
1238         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1239         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1240         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1241         ntracks[ilayer]++;
1242       }
1243
1244
1245     } // loop over tracks in layer ilayer+1
1246
1247     //loop over track candidates for the current layer
1248     //
1249     //
1250     Int_t accepted=0;
1251
1252     Int_t golden=0;
1253     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1254       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1255       if (normalizedchi2[itrack] <
1256           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1257       if (ilayer>4) {
1258         accepted++;
1259       } else {
1260         if (constrain) { // constrain
1261           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1262             accepted++;
1263         } else {        // !constrain
1264           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1265             accepted++;
1266         }
1267       }
1268     }
1269     // sort tracks by increasing normalized chi2
1270     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1271     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1272     if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1273     if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1274   } // end loop over layers
1275
1276
1277   //
1278   // Now select tracks to be kept
1279   //
1280   Int_t max = constrain ? 20 : 5;
1281
1282   // tracks that reach layer 0 (SPD inner)
1283   for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1284     AliITStrackMI & track= tracks[0][nindexes[0][i]];
1285     if (track.GetNumberOfClusters()<2) continue;
1286     if (!constrain && track.GetNormChi2(0) >
1287         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1288       continue;
1289     }
1290     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1291   }
1292
1293   // tracks that reach layer 1 (SPD outer)
1294   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1295     AliITStrackMI & track= tracks[1][nindexes[1][i]];
1296     if (track.GetNumberOfClusters()<4) continue;
1297     if (!constrain && track.GetNormChi2(1) >
1298         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1299     if (constrain) track.IncrementNSkipped();
1300     if (!constrain) {
1301       track.SetD(0,track.GetD(GetX(),GetY()));
1302       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1303       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1304         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1305       }
1306     }
1307     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1308   }
1309
1310   // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1311   if (!constrain){
1312     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1313       AliITStrackMI & track= tracks[2][nindexes[2][i]];
1314       if (track.GetNumberOfClusters()<3) continue;
1315       if (!constrain && track.GetNormChi2(2) >
1316           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1317       if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1318       if (!constrain){
1319         track.SetD(0,track.GetD(GetX(),GetY()));
1320         track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1321         if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1322           track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1323         }
1324       }
1325       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1326     }
1327   }
1328
1329   if (!constrain) {
1330     //
1331     // register best track of each layer - important for V0 finder
1332     //
1333     for (Int_t ilayer=0;ilayer<5;ilayer++){
1334       if (ntracks[ilayer]==0) continue;
1335       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1336       if (track.GetNumberOfClusters()<1) continue;
1337       CookLabel(&track,0);
1338       bestarray->AddAt(new AliITStrackMI(track),ilayer);
1339     }
1340   }
1341   //
1342   // update TPC V0 information
1343   //
1344   if (otrack->GetESDtrack()->GetV0Index(0)>0){
1345     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1346     for (Int_t i=0;i<3;i++){
1347       Int_t  index = otrack->GetESDtrack()->GetV0Index(i);
1348       if (index==0) break;
1349       AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1350       if (vertex->GetStatus()<0) continue;     // rejected V0
1351       //
1352       if (otrack->GetSign()>0) {
1353         vertex->SetIndex(0,esdindex);
1354       }
1355       else{
1356         vertex->SetIndex(1,esdindex);
1357       }
1358       //find nearest layer with track info
1359       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1360       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1361       Int_t nearest     = nearestold;
1362       for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1363         if (ntracks[nearest]==0){
1364           nearest = ilayer;
1365         }
1366       }
1367       //
1368       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1369       if (nearestold<5&&nearest<5){
1370         Bool_t accept = track.GetNormChi2(nearest)<10;
1371         if (accept){
1372           if (track.GetSign()>0) {
1373             vertex->SetParamP(track);
1374             vertex->Update(fprimvertex);
1375             //vertex->SetIndex(0,track.fESDtrack->GetID());
1376             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1377           }else{
1378             vertex->SetParamN(track);
1379             vertex->Update(fprimvertex);
1380             //vertex->SetIndex(1,track.fESDtrack->GetID());
1381             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1382           }
1383           vertex->SetStatus(vertex->GetStatus()+1);
1384         }else{
1385           //vertex->SetStatus(-2);  // reject V0  - not enough clusters
1386         }
1387       }
1388     }
1389   }
1390
1391 }
1392 //------------------------------------------------------------------------
1393 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1394 {
1395   //--------------------------------------------------------------------
1396   //
1397   //
1398   return fgLayers[layer];
1399 }
1400 //------------------------------------------------------------------------
1401 AliITStrackerMI::AliITSlayer::AliITSlayer():
1402 fR(0),
1403 fPhiOffset(0),
1404 fNladders(0),
1405 fZOffset(0),
1406 fNdetectors(0),
1407 fDetectors(0),
1408 fN(0),
1409 fDy5(0),
1410 fDy10(0),
1411 fDy20(0),
1412 fClustersCs(0),
1413 fClusterIndexCs(0),
1414 fYcs(0),
1415 fZcs(0),
1416 fNcs(0),
1417 fCurrentSlice(-1),
1418 fZmax(0),
1419 fYmin(0),
1420 fYmax(0),
1421 fI(0),
1422 fImax(0),
1423 fSkip(0),
1424 fAccepted(0),
1425 fRoad(0){
1426   //--------------------------------------------------------------------
1427   //default AliITSlayer constructor
1428   //--------------------------------------------------------------------
1429   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1430     fClusterWeight[i]=0;
1431     fClusterTracks[0][i]=-1;
1432     fClusterTracks[1][i]=-1;
1433     fClusterTracks[2][i]=-1;
1434     fClusterTracks[3][i]=-1;
1435   }
1436 }
1437 //------------------------------------------------------------------------
1438 AliITStrackerMI::AliITSlayer::
1439 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1440 fR(r),
1441 fPhiOffset(p),
1442 fNladders(nl),
1443 fZOffset(z),
1444 fNdetectors(nd),
1445 fDetectors(0),
1446 fN(0),
1447 fDy5(0),
1448 fDy10(0),
1449 fDy20(0),
1450 fClustersCs(0),
1451 fClusterIndexCs(0),
1452 fYcs(0),
1453 fZcs(0),
1454 fNcs(0),
1455 fCurrentSlice(-1),
1456 fZmax(0),
1457 fYmin(0),
1458 fYmax(0),
1459 fI(0),
1460 fImax(0),
1461 fSkip(0),
1462 fAccepted(0),
1463 fRoad(0) {
1464   //--------------------------------------------------------------------
1465   //main AliITSlayer constructor
1466   //--------------------------------------------------------------------
1467   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1468   fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1469 }
1470 //------------------------------------------------------------------------
1471 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1472 fR(layer.fR),
1473 fPhiOffset(layer.fPhiOffset),
1474 fNladders(layer.fNladders),
1475 fZOffset(layer.fZOffset),
1476 fNdetectors(layer.fNdetectors),
1477 fDetectors(layer.fDetectors),
1478 fN(layer.fN),
1479 fDy5(layer.fDy5),
1480 fDy10(layer.fDy10),
1481 fDy20(layer.fDy20),
1482 fClustersCs(layer.fClustersCs),
1483 fClusterIndexCs(layer.fClusterIndexCs),
1484 fYcs(layer.fYcs),
1485 fZcs(layer.fZcs),
1486 fNcs(layer.fNcs),
1487 fCurrentSlice(layer.fCurrentSlice),
1488 fZmax(layer.fZmax),
1489 fYmin(layer.fYmin),
1490 fYmax(layer.fYmax),
1491 fI(layer.fI),
1492 fImax(layer.fImax),
1493 fSkip(layer.fSkip),
1494 fAccepted(layer.fAccepted),
1495 fRoad(layer.fRoad){
1496   //Copy constructor
1497 }
1498 //------------------------------------------------------------------------
1499 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1500   //--------------------------------------------------------------------
1501   // AliITSlayer destructor
1502   //--------------------------------------------------------------------
1503   delete [] fDetectors;
1504   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1505   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1506     fClusterWeight[i]=0;
1507     fClusterTracks[0][i]=-1;
1508     fClusterTracks[1][i]=-1;
1509     fClusterTracks[2][i]=-1;
1510     fClusterTracks[3][i]=-1;
1511   }
1512 }
1513 //------------------------------------------------------------------------
1514 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1515   //--------------------------------------------------------------------
1516   // This function removes loaded clusters
1517   //--------------------------------------------------------------------
1518   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1519   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1520     fClusterWeight[i]=0;
1521     fClusterTracks[0][i]=-1;
1522     fClusterTracks[1][i]=-1;
1523     fClusterTracks[2][i]=-1;
1524     fClusterTracks[3][i]=-1;
1525   }
1526
1527   fN=0;
1528   fI=0;
1529 }
1530 //------------------------------------------------------------------------
1531 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1532   //--------------------------------------------------------------------
1533   // This function reset weights of the clusters
1534   //--------------------------------------------------------------------
1535   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1536     fClusterWeight[i]=0;
1537     fClusterTracks[0][i]=-1;
1538     fClusterTracks[1][i]=-1;
1539     fClusterTracks[2][i]=-1;
1540     fClusterTracks[3][i]=-1;
1541   }
1542   for (Int_t i=0; i<fN;i++) {
1543     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1544     if (cl&&cl->IsUsed()) cl->Use();
1545   }
1546
1547 }
1548 //------------------------------------------------------------------------
1549 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1550   //--------------------------------------------------------------------
1551   // This function calculates the road defined by the cluster density
1552   //--------------------------------------------------------------------
1553   Int_t n=0;
1554   for (Int_t i=0; i<fN; i++) {
1555      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1556   }
1557   if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1558 }
1559 //------------------------------------------------------------------------
1560 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1561   //--------------------------------------------------------------------
1562   //This function adds a cluster to this layer
1563   //--------------------------------------------------------------------
1564   if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1565     ::Error("InsertCluster","Too many clusters !\n");
1566     return 1;
1567   }
1568   fCurrentSlice=-1;
1569   fClusters[fN]=cl;
1570   fN++;
1571   AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1572   if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1573   if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1574   if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1575   if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1576
1577   return 0;
1578 }
1579 //------------------------------------------------------------------------
1580 void  AliITStrackerMI::AliITSlayer::SortClusters()
1581 {
1582   //
1583   //sort clusters
1584   //
1585   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1586   Float_t *z                = new Float_t[fN];
1587   Int_t   * index           = new Int_t[fN];
1588   //
1589   for (Int_t i=0;i<fN;i++){
1590     z[i] = fClusters[i]->GetZ();
1591   }
1592   TMath::Sort(fN,z,index,kFALSE);
1593   for (Int_t i=0;i<fN;i++){
1594     clusters[i] = fClusters[index[i]];
1595   }
1596   //
1597   for (Int_t i=0;i<fN;i++){
1598     fClusters[i] = clusters[i];
1599     fZ[i]        = fClusters[i]->GetZ();
1600     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1601     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1602     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1603     fY[i] = y;
1604   }
1605   delete[] index;
1606   delete[] z;
1607   delete[] clusters;
1608   //
1609
1610   fYB[0]=10000000;
1611   fYB[1]=-10000000;
1612   for (Int_t i=0;i<fN;i++){
1613     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1614     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1615     fClusterIndex[i] = i;
1616   }
1617   //
1618   // fill slices
1619   fDy5 = (fYB[1]-fYB[0])/5.;
1620   fDy10 = (fYB[1]-fYB[0])/10.;
1621   fDy20 = (fYB[1]-fYB[0])/20.;
1622   for (Int_t i=0;i<6;i++)  fN5[i] =0;
1623   for (Int_t i=0;i<11;i++) fN10[i]=0;
1624   for (Int_t i=0;i<21;i++) fN20[i]=0;
1625   //
1626   for (Int_t i=0;i<6;i++) {fBy5[i][0] =  fYB[0]+(i-0.75)*fDy5; fBy5[i][1] =  fYB[0]+(i+0.75)*fDy5;}
1627   for (Int_t i=0;i<11;i++) {fBy10[i][0] =  fYB[0]+(i-0.75)*fDy10; fBy10[i][1] =  fYB[0]+(i+0.75)*fDy10;}
1628   for (Int_t i=0;i<21;i++) {fBy20[i][0] =  fYB[0]+(i-0.75)*fDy20; fBy20[i][1] =  fYB[0]+(i+0.75)*fDy20;}
1629   //
1630   //
1631   for (Int_t i=0;i<fN;i++)
1632     for (Int_t irot=-1;irot<=1;irot++){
1633       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1634       // slice 5
1635       for (Int_t slice=0; slice<6;slice++){
1636         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1637           fClusters5[slice][fN5[slice]] = fClusters[i];
1638           fY5[slice][fN5[slice]] = curY;
1639           fZ5[slice][fN5[slice]] = fZ[i];
1640           fClusterIndex5[slice][fN5[slice]]=i;
1641           fN5[slice]++;
1642         }
1643       }
1644       // slice 10
1645       for (Int_t slice=0; slice<11;slice++){
1646         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1647           fClusters10[slice][fN10[slice]] = fClusters[i];
1648           fY10[slice][fN10[slice]] = curY;
1649           fZ10[slice][fN10[slice]] = fZ[i];
1650           fClusterIndex10[slice][fN10[slice]]=i;
1651           fN10[slice]++;
1652         }
1653       }
1654       // slice 20
1655       for (Int_t slice=0; slice<21;slice++){
1656         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1657           fClusters20[slice][fN20[slice]] = fClusters[i];
1658           fY20[slice][fN20[slice]] = curY;
1659           fZ20[slice][fN20[slice]] = fZ[i];
1660           fClusterIndex20[slice][fN20[slice]]=i;
1661           fN20[slice]++;
1662         }
1663       }
1664     }
1665
1666   //
1667   // consistency check
1668   //
1669   for (Int_t i=0;i<fN-1;i++){
1670     if (fZ[i]>fZ[i+1]){
1671       printf("Bug\n");
1672     }
1673   }
1674   //
1675   for (Int_t slice=0;slice<21;slice++)
1676   for (Int_t i=0;i<fN20[slice]-1;i++){
1677     if (fZ20[slice][i]>fZ20[slice][i+1]){
1678       printf("Bug\n");
1679     }
1680   }
1681
1682
1683 }
1684 //------------------------------------------------------------------------
1685 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1686   //--------------------------------------------------------------------
1687   // This function returns the index of the nearest cluster
1688   //--------------------------------------------------------------------
1689   Int_t ncl=0;
1690   const Float_t *zcl;
1691   if (fCurrentSlice<0) {
1692     ncl = fN;
1693     zcl   = fZ;
1694   }
1695   else{
1696     ncl   = fNcs;
1697     zcl   = fZcs;;
1698   }
1699
1700   if (ncl==0) return 0;
1701   Int_t b=0, e=ncl-1, m=(b+e)/2;
1702   for (; b<e; m=(b+e)/2) {
1703     //    if (z > fClusters[m]->GetZ()) b=m+1;
1704     if (z > zcl[m]) b=m+1;
1705     else e=m;
1706   }
1707   return m;
1708 }
1709 //------------------------------------------------------------------------
1710 Bool_t AliITStrackerMI::ComputeRoad(AliITStrackMI* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1711   //--------------------------------------------------------------------
1712   // This function computes the rectangular road for this track
1713   //--------------------------------------------------------------------
1714
1715
1716   AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1717   // take into account the misalignment: propagate track to misaligned detector plane
1718   if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1719
1720   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1721                     TMath::Sqrt(track->GetSigmaZ2() +
1722                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1723                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1724                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1725   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1726                     TMath::Sqrt(track->GetSigmaY2() +
1727                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1728                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1729                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1730
1731   // track at boundary between detectors, enlarge road
1732   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1733   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1734        (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1735        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1736        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1737     Float_t tgl = TMath::Abs(track->GetTgl());
1738     if (tgl > 1.) tgl=1.;
1739     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1740     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1741     Float_t snp = TMath::Abs(track->GetSnp());
1742     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1743     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1744   } // boundary
1745
1746   // add to the road a term (up to 2-3 mm) to deal with misalignments
1747   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1748   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1749
1750   Double_t r = fgLayers[ilayer].GetR();
1751   zmin = track->GetZ() - dz;
1752   zmax = track->GetZ() + dz;
1753   ymin = track->GetY() + r*det.GetPhi() - dy;
1754   ymax = track->GetY() + r*det.GetPhi() + dy;
1755
1756   // bring track back to idead detector plane
1757   if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1758
1759   return kTRUE;
1760 }
1761 //------------------------------------------------------------------------
1762 void AliITStrackerMI::AliITSlayer::
1763 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1764   //--------------------------------------------------------------------
1765   // This function sets the "window"
1766   //--------------------------------------------------------------------
1767
1768   Double_t circle=2*TMath::Pi()*fR;
1769   fYmin = ymin; fYmax =ymax;
1770   Float_t ymiddle = (fYmin+fYmax)*0.5;
1771   if (ymiddle<fYB[0]) {
1772     fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1773   } else if (ymiddle>fYB[1]) {
1774     fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1775   }
1776
1777   //
1778   fCurrentSlice =-1;
1779   // defualt take all
1780   fClustersCs = fClusters;
1781   fClusterIndexCs = fClusterIndex;
1782   fYcs  = fY;
1783   fZcs  = fZ;
1784   fNcs  = fN;
1785   //
1786   //is in 20 slice?
1787   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1788     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1789     if (slice<0) slice=0;
1790     if (slice>20) slice=20;
1791     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1792     if (isOK) {
1793       fCurrentSlice=slice;
1794       fClustersCs = fClusters20[fCurrentSlice];
1795       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1796       fYcs  = fY20[fCurrentSlice];
1797       fZcs  = fZ20[fCurrentSlice];
1798       fNcs  = fN20[fCurrentSlice];
1799     }
1800   }
1801   //
1802   //is in 10 slice?
1803   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1804     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1805     if (slice<0) slice=0;
1806     if (slice>10) slice=10;
1807     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1808     if (isOK) {
1809       fCurrentSlice=slice;
1810       fClustersCs = fClusters10[fCurrentSlice];
1811       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1812       fYcs  = fY10[fCurrentSlice];
1813       fZcs  = fZ10[fCurrentSlice];
1814       fNcs  = fN10[fCurrentSlice];
1815     }
1816   }
1817   //
1818   //is in 5 slice?
1819   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1820     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1821     if (slice<0) slice=0;
1822     if (slice>5) slice=5;
1823     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1824     if (isOK) {
1825       fCurrentSlice=slice;
1826       fClustersCs = fClusters5[fCurrentSlice];
1827       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1828       fYcs  = fY5[fCurrentSlice];
1829       fZcs  = fZ5[fCurrentSlice];
1830       fNcs  = fN5[fCurrentSlice];
1831     }
1832   }
1833   //
1834   fI=FindClusterIndex(zmin); fZmax=zmax;
1835   fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1836   fSkip = 0;
1837   fAccepted =0;
1838
1839   return;
1840 }
1841 //------------------------------------------------------------------------
1842 Int_t AliITStrackerMI::AliITSlayer::
1843 FindDetectorIndex(Double_t phi, Double_t z) const {
1844   //--------------------------------------------------------------------
1845   //This function finds the detector crossed by the track
1846   //--------------------------------------------------------------------
1847   Double_t dphi;
1848   if (fZOffset<0)            // old geometry
1849     dphi = -(phi-fPhiOffset);
1850   else                       // new geometry
1851     dphi = phi-fPhiOffset;
1852
1853
1854   if      (dphi <  0) dphi += 2*TMath::Pi();
1855   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1856   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1857   if (np>=fNladders) np-=fNladders;
1858   if (np<0)          np+=fNladders;
1859
1860
1861   Double_t dz=fZOffset-z;
1862   Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1863   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1864   if (nz>=fNdetectors) return -1;
1865   if (nz<0)            return -1;
1866
1867   // ad hoc correction for 3rd ladder of SDD inner layer,
1868   // which is reversed (rotated by pi around local y)
1869   // this correction is OK only from AliITSv11Hybrid onwards
1870   if (GetR()>12. && GetR()<20.) { // SDD inner
1871     if(np==2) { // 3rd ladder
1872       nz = (fNdetectors-1) - nz;
1873     }
1874   }
1875   //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
1876
1877
1878   return np*fNdetectors + nz;
1879 }
1880 //------------------------------------------------------------------------
1881 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1882 {
1883   //--------------------------------------------------------------------
1884   // This function returns clusters within the "window"
1885   //--------------------------------------------------------------------
1886
1887   if (fCurrentSlice<0) {
1888     Double_t rpi2 = 2.*fR*TMath::Pi();
1889     for (Int_t i=fI; i<fImax; i++) {
1890       Double_t y = fY[i];
1891       if (fYmax<y) y -= rpi2;
1892       if (fYmin>y) y += rpi2;
1893       if (y<fYmin) continue;
1894       if (y>fYmax) continue;
1895       if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1896       ci=i;
1897       if (!test) fI=i+1;
1898       return fClusters[i];
1899     }
1900   } else {
1901     for (Int_t i=fI; i<fImax; i++) {
1902       if (fYcs[i]<fYmin) continue;
1903       if (fYcs[i]>fYmax) continue;
1904       if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1905       ci=fClusterIndexCs[i];
1906       if (!test) fI=i+1;
1907       return fClustersCs[i];
1908     }
1909   }
1910   return 0;
1911 }
1912 //------------------------------------------------------------------------
1913 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1914 const {
1915   //--------------------------------------------------------------------
1916   // This function returns the layer thickness at this point (units X0)
1917   //--------------------------------------------------------------------
1918   Double_t d=0.0085;
1919   x0=AliITSRecoParam::GetX0Air();
1920   if (43<fR&&fR<45) { //SSD2
1921      Double_t dd=0.0034;
1922      d=dd;
1923      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1924      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1925      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1926      for (Int_t i=0; i<12; i++) {
1927        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1928           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1929           d+=0.0034;
1930           break;
1931        }
1932        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1933           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1934           d+=0.0034;
1935           break;
1936        }
1937        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1938        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1939      }
1940   } else
1941   if (37<fR&&fR<41) { //SSD1
1942      Double_t dd=0.0034;
1943      d=dd;
1944      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1945      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1946      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1947      for (Int_t i=0; i<11; i++) {
1948        if (TMath::Abs(z-3.9*i)<0.15) {
1949           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1950           d+=dd;
1951           break;
1952        }
1953        if (TMath::Abs(z+3.9*i)<0.15) {
1954           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1955           d+=dd;
1956           break;
1957        }
1958        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1959        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1960      }
1961   } else
1962   if (13<fR&&fR<26) { //SDD
1963      Double_t dd=0.0033;
1964      d=dd;
1965      if (TMath::Abs(y-0.00)>3.30) d+=dd;
1966
1967      if (TMath::Abs(y-1.80)<0.55) {
1968         d+=0.016;
1969         for (Int_t j=0; j<20; j++) {
1970           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1971           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1972         }
1973      }
1974      if (TMath::Abs(y+1.80)<0.55) {
1975         d+=0.016;
1976         for (Int_t j=0; j<20; j++) {
1977           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1978           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1979         }
1980      }
1981
1982      for (Int_t i=0; i<4; i++) {
1983        if (TMath::Abs(z-7.3*i)<0.60) {
1984           d+=dd;
1985           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1986           break;
1987        }
1988        if (TMath::Abs(z+7.3*i)<0.60) {
1989           d+=dd;
1990           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1991           break;
1992        }
1993      }
1994   } else
1995   if (6<fR&&fR<8) {   //SPD2
1996      Double_t dd=0.0063; x0=21.5;
1997      d=dd;
1998      if (TMath::Abs(y-3.08)>0.5) d+=dd;
1999      if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2000   } else
2001   if (3<fR&&fR<5) {   //SPD1
2002      Double_t dd=0.0063; x0=21.5;
2003      d=dd;
2004      if (TMath::Abs(y+0.21)>0.6) d+=dd;
2005      if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2006   }
2007
2008   return d;
2009 }
2010 //------------------------------------------------------------------------
2011 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2012 fR(det.fR),
2013 fRmisal(det.fRmisal),
2014 fPhi(det.fPhi),
2015 fSinPhi(det.fSinPhi),
2016 fCosPhi(det.fCosPhi),
2017 fYmin(det.fYmin),
2018 fYmax(det.fYmax),
2019 fZmin(det.fZmin),
2020 fZmax(det.fZmax),
2021 fIsBad(det.fIsBad),
2022 fNChips(det.fNChips),
2023 fChipIsBad(det.fChipIsBad)
2024 {
2025   //Copy constructor
2026 }
2027 //------------------------------------------------------------------------
2028 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2029                                                const AliITSDetTypeRec *detTypeRec)
2030 {
2031   //--------------------------------------------------------------------
2032   // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2033   //--------------------------------------------------------------------
2034
2035   // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2036   // while in the tracker they start from 0 for each layer
2037   for(Int_t il=0; il<ilayer; il++)
2038     idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2039
2040   Int_t detType;
2041   if (ilayer==0 || ilayer==1) {        // ----------  SPD
2042     detType = 0;
2043   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
2044     detType = 1;
2045   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
2046     detType = 2;
2047   } else {
2048     printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2049     return;
2050   }
2051
2052   // Get calibration from AliITSDetTypeRec
2053   AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2054   calib->SetModuleIndex(idet);
2055   AliITSCalibration *calibSPDdead = 0;
2056   if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2057   if (calib->IsBad() ||
2058       (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2059     {
2060       SetBad();
2061       //      printf("lay %d bad %d\n",ilayer,idet);
2062     }
2063
2064   // Get segmentation from AliITSDetTypeRec
2065   AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2066
2067   // Read info about bad chips
2068   fNChips = segm->GetMaximumChipIndex()+1;
2069   //printf("ilayer %d  detType %d idet %d fNChips %d %d  GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2070   if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2071   fChipIsBad = new Bool_t[fNChips];
2072   for (Int_t iCh=0;iCh<fNChips;iCh++) {
2073     fChipIsBad[iCh] = calib->IsChipBad(iCh);
2074     if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2075     //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2076   }
2077
2078   return;
2079 }
2080 //------------------------------------------------------------------------
2081 Double_t AliITStrackerMI::GetEffectiveThickness()
2082 {
2083   //--------------------------------------------------------------------
2084   // Returns the thickness between the current layer and the vertex (units X0)
2085   //--------------------------------------------------------------------
2086
2087   if(fUseTGeo!=0) {
2088     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2089     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2090     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2091   }
2092
2093   // beam pipe
2094   Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2095   Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2096
2097   // layers
2098   Double_t x0=0;
2099   Double_t xn=fgLayers[fI].GetR();
2100   for (Int_t i=0; i<fI; i++) {
2101     Double_t xi=fgLayers[i].GetR();
2102     Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2103     d+=dLayer*xi*xi;
2104   }
2105
2106   // shields
2107   if (fI>1) {
2108     Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2109     d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2110   }
2111   if (fI>3) {
2112     Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2113     d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2114   }
2115   return d/(xn*xn);
2116 }
2117 //------------------------------------------------------------------------
2118 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2119   //-------------------------------------------------------------------
2120   // This function returns number of clusters within the "window"
2121   //--------------------------------------------------------------------
2122   Int_t ncl=0;
2123   for (Int_t i=fI; i<fN; i++) {
2124     const AliITSRecPoint *c=fClusters[i];
2125     if (c->GetZ() > fZmax) break;
2126     if (c->IsUsed()) continue;
2127     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2128     Double_t y=fR*det.GetPhi() + c->GetY();
2129
2130     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2131     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2132
2133     if (y<fYmin) continue;
2134     if (y>fYmax) continue;
2135     ncl++;
2136   }
2137   return ncl;
2138 }
2139 //------------------------------------------------------------------------
2140 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2141                                 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2142 {
2143   //--------------------------------------------------------------------
2144   // This function refits the track "track" at the position "x" using
2145   // the clusters from "clusters"
2146   // If "extra"==kTRUE,
2147   //    the clusters from overlapped modules get attached to "track"
2148   // If "planeff"==kTRUE,
2149   //    special approach for plane efficiency evaluation is applyed
2150   //--------------------------------------------------------------------
2151
2152   Int_t index[AliITSgeomTGeo::kNLayers];
2153   Int_t k;
2154   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2155   Int_t nc=clusters->GetNumberOfClusters();
2156   for (k=0; k<nc; k++) {
2157     Int_t idx=clusters->GetClusterIndex(k);
2158     Int_t ilayer=(idx&0xf0000000)>>28;
2159     index[ilayer]=idx;
2160   }
2161
2162   return RefitAt(xx,track,index,extra,planeeff); // call the method below
2163 }
2164 //------------------------------------------------------------------------
2165 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2166                                 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2167 {
2168   //--------------------------------------------------------------------
2169   // This function refits the track "track" at the position "x" using
2170   // the clusters from array
2171   // If "extra"==kTRUE,
2172   //    the clusters from overlapped modules get attached to "track"
2173   // If "planeff"==kTRUE,
2174   //    special approach for plane efficiency evaluation is applyed
2175   //--------------------------------------------------------------------
2176   Int_t index[AliITSgeomTGeo::kNLayers];
2177   Int_t k;
2178   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2179   //
2180   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2181     index[k]=clusters[k];
2182   }
2183
2184   // special for cosmics: check which the innermost layer crossed
2185   // by the track
2186   Int_t innermostlayer=5;
2187   Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2188   for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2189     if(drphi < fgLayers[innermostlayer].GetR()) break;
2190   }
2191   //printf(" drphi  %f  innermost %d\n",drphi,innermostlayer);
2192
2193   Int_t modstatus=1; // found
2194   Float_t xloc,zloc;
2195   Int_t from, to, step;
2196   if (xx > track->GetX()) {
2197       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2198       step=+1;
2199   } else {
2200       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2201       step=-1;
2202   }
2203   TString dir = (step>0 ? "outward" : "inward");
2204
2205   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2206      AliITSlayer &layer=fgLayers[ilayer];
2207      Double_t r=layer.GetR();
2208      if (step<0 && xx>r) break;
2209
2210      // material between SSD and SDD, SDD and SPD
2211      Double_t hI=ilayer-0.5*step;
2212      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2213        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2214      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2215        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2216
2217
2218      Double_t oldGlobXYZ[3];
2219      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2220
2221      Double_t phi,z;
2222      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2223
2224      Int_t idet=layer.FindDetectorIndex(phi,z);
2225
2226      // check if we allow a prolongation without point for large-eta tracks
2227      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2228      if (skip==2) {
2229        modstatus = 4; // out in z
2230        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2231          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2232        }
2233        // cross layer
2234        // apply correction for material of the current layer
2235        // add time if going outward
2236        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2237        continue;
2238      }
2239
2240      if (idet<0) return kFALSE;
2241
2242      const AliITSdetector &det=layer.GetDetector(idet);
2243      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2244
2245      track->SetDetectorIndex(idet);
2246      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2247
2248      Double_t dz,zmin,zmax,dy,ymin,ymax;
2249
2250      const AliITSRecPoint *clAcc=0;
2251      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2252
2253      Int_t idx=index[ilayer];
2254      if (idx>=0) { // cluster in this layer
2255        modstatus = 6; // no refit
2256        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2257        if (cl) {
2258          if (idet != cl->GetDetectorIndex()) {
2259            idet=cl->GetDetectorIndex();
2260            const AliITSdetector &detc=layer.GetDetector(idet);
2261            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2262            track->SetDetectorIndex(idet);
2263            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2264          }
2265          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2266          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2267          if (chi2<maxchi2) {
2268            clAcc=cl;
2269            maxchi2=chi2;
2270            modstatus = 1; // found
2271          } else {
2272             return kFALSE; //
2273          }
2274        }
2275      } else { // no cluster in this layer
2276        if (skip==1) {
2277          modstatus = 3; // skipped
2278          // Plane Eff determination:
2279          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2280            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2281               UseTrackForPlaneEff(track,ilayer);
2282          }
2283        } else {
2284          modstatus = 5; // no cls in road
2285          // check dead
2286          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2287          dz = 0.5*(zmax-zmin);
2288          dy = 0.5*(ymax-ymin);
2289          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2290          if (dead==1) modstatus = 7; // holes in z in SPD
2291          if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2292        }
2293      }
2294
2295      if (clAcc) {
2296        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2297        track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2298      }
2299      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2300
2301
2302      if (extra) { // search for extra clusters in overlapped modules
2303        AliITStrackV2 tmp(*track);
2304        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2305        layer.SelectClusters(zmin,zmax,ymin,ymax);
2306
2307        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2308        Int_t idetExtra=-1;
2309        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2310        Double_t tolerance=0.1;
2311        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2312          // only clusters in another module! (overlaps)
2313          idetExtra = clExtra->GetDetectorIndex();
2314          if (idet == idetExtra) continue;
2315
2316          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2317
2318          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2319          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2320          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2321          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2322
2323          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2324          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2325        }
2326        if (cci>=0) {
2327          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2328          track->SetExtraModule(ilayer,idetExtra);
2329        }
2330      } // end search for extra clusters in overlapped modules
2331
2332      // Correct for material of the current layer
2333      // cross material
2334      // add time if going outward
2335      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2336
2337   } // end loop on layers
2338
2339   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2340
2341   return kTRUE;
2342 }
2343 //------------------------------------------------------------------------
2344 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2345 {
2346   //
2347   // calculate normalized chi2
2348   //  return NormalizedChi2(track,0);
2349   Float_t chi2 = 0;
2350   Float_t sum=0;
2351   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2352   //  track->fdEdxMismatch=0;
2353   Float_t dedxmismatch =0;
2354   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2355   if (mode<100){
2356     for (Int_t i = 0;i<6;i++){
2357       if (track->GetClIndex(i)>=0){
2358         Float_t cerry, cerrz;
2359         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2360         else
2361           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2362         cerry*=cerry;
2363         cerrz*=cerrz;
2364         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2365         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2366           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2367           if (ratio<0.5) {
2368             cchi2+=(0.5-ratio)*10.;
2369             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2370             dedxmismatch+=(0.5-ratio)*10.;
2371           }
2372         }
2373         if (i<2 ||i>3){
2374           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2375           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2376           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2377           if (i<2) chi2+=2*cl->GetDeltaProbability();
2378         }
2379         chi2+=cchi2;
2380         sum++;
2381       }
2382     }
2383     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2384       track->SetdEdxMismatch(dedxmismatch);
2385     }
2386   }
2387   else{
2388     for (Int_t i = 0;i<4;i++){
2389       if (track->GetClIndex(i)>=0){
2390         Float_t cerry, cerrz;
2391         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2392         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2393         cerry*=cerry;
2394         cerrz*=cerrz;
2395         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2396         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2397         sum++;
2398       }
2399     }
2400     for (Int_t i = 4;i<6;i++){
2401       if (track->GetClIndex(i)>=0){
2402         Float_t cerry, cerrz;
2403         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2404         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2405         cerry*=cerry;
2406         cerrz*=cerrz;
2407         Float_t cerryb, cerrzb;
2408         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2409         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2410         cerryb*=cerryb;
2411         cerrzb*=cerrzb;
2412         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2413         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2414         sum++;
2415       }
2416     }
2417   }
2418   if (track->GetESDtrack()->GetTPCsignal()>85){
2419     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2420     if (ratio<0.5) {
2421       chi2+=(0.5-ratio)*5.;
2422     }
2423     if (ratio>2){
2424       chi2+=(ratio-2.0)*3;
2425     }
2426   }
2427   //
2428   Double_t match = TMath::Sqrt(track->GetChi22());
2429   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2430   if (!track->GetConstrain()) {
2431     if (track->GetNumberOfClusters()>2) {
2432       match/=track->GetNumberOfClusters()-2.;
2433     } else {
2434       match=0;
2435     }
2436   }
2437   if (match<0) match=0;
2438   Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2439   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2440     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2441                                 1./(1.+track->GetNSkipped()));
2442
2443  return normchi2;
2444 }
2445 //------------------------------------------------------------------------
2446 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2447 {
2448   //
2449   // return matching chi2 between two tracks
2450   Double_t largeChi2=1000.;
2451
2452   AliITStrackMI track3(*track2);
2453   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2454   TMatrixD vec(5,1);
2455   vec(0,0)=track1->GetY()   - track3.GetY();
2456   vec(1,0)=track1->GetZ()   - track3.GetZ();
2457   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2458   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2459   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2460   //
2461   TMatrixD cov(5,5);
2462   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2463   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2464   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2465   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2466   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2467
2468   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2469   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2470   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2471   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2472   //
2473   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2474   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2475   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2476   //
2477   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2478   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2479   //
2480   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2481
2482   cov.Invert();
2483   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2484   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2485   return chi2(0,0);
2486 }
2487 //------------------------------------------------------------------------
2488 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2489 {
2490   //
2491   //  return probability that given point (characterized by z position and error)
2492   //  is in SPD dead zone
2493   //
2494   Double_t probability = 0.;
2495   Double_t absz = TMath::Abs(zpos);
2496   Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2497                                   0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2498   if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2499   Double_t zmin, zmax;
2500   if (zpos<-6.) { // dead zone at z = -7
2501     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2502     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2503   } else if (zpos>6.) { // dead zone at z = +7
2504     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2505     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2506   } else if (absz<2.) { // dead zone at z = 0
2507     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2508     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2509   } else {
2510     zmin = 0.;
2511     zmax = 0.;
2512   }
2513   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2514   // dead zone)
2515   probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2516                       TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2517   return probability;
2518 }
2519 //------------------------------------------------------------------------
2520 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2521 {
2522   //
2523   // calculate normalized chi2
2524   Float_t chi2[6];
2525   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2526   Float_t ncl = 0;
2527   for (Int_t i = 0;i<6;i++){
2528     if (TMath::Abs(track->GetDy(i))>0){
2529       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2530       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2531       ncl++;
2532     }
2533     else{chi2[i]=10000;}
2534   }
2535   Int_t index[6];
2536   TMath::Sort(6,chi2,index,kFALSE);
2537   Float_t max = float(ncl)*fac-1.;
2538   Float_t sumchi=0, sumweight=0;
2539   for (Int_t i=0;i<max+1;i++){
2540     Float_t weight = (i<max)?1.:(max+1.-i);
2541     sumchi+=weight*chi2[index[i]];
2542     sumweight+=weight;
2543   }
2544   Double_t normchi2 = sumchi/sumweight;
2545   return normchi2;
2546 }
2547 //------------------------------------------------------------------------
2548 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2549 {
2550   //
2551   // calculate normalized chi2
2552   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2553   Int_t npoints = 0;
2554   Double_t res =0;
2555   for (Int_t i=0;i<6;i++){
2556     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2557     Double_t sy1 = forwardtrack->GetSigmaY(i);
2558     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2559     Double_t sy2 = backtrack->GetSigmaY(i);
2560     Double_t sz2 = backtrack->GetSigmaZ(i);
2561     if (i<2){ sy2=1000.;sz2=1000;}
2562     //
2563     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2564     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2565     //
2566     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2567     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2568     //
2569     res+= nz0*nz0+ny0*ny0;
2570     npoints++;
2571   }
2572   if (npoints>1) return
2573                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2574                    //2*forwardtrack->fNUsed+
2575                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2576                                   1./(1.+forwardtrack->GetNSkipped()));
2577   return 1000;
2578 }
2579 //------------------------------------------------------------------------
2580 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2581   //--------------------------------------------------------------------
2582   //       Return pointer to a given cluster
2583   //--------------------------------------------------------------------
2584   Int_t l=(index & 0xf0000000) >> 28;
2585   Int_t c=(index & 0x0fffffff) >> 00;
2586   return fgLayers[l].GetWeight(c);
2587 }
2588 //------------------------------------------------------------------------
2589 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2590 {
2591   //---------------------------------------------
2592   // register track to the list
2593   //
2594   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2595   //
2596   //
2597   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2598     Int_t index = track->GetClusterIndex(icluster);
2599     Int_t l=(index & 0xf0000000) >> 28;
2600     Int_t c=(index & 0x0fffffff) >> 00;
2601     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2602     for (Int_t itrack=0;itrack<4;itrack++){
2603       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2604         fgLayers[l].SetClusterTracks(itrack,c,id);
2605         break;
2606       }
2607     }
2608   }
2609 }
2610 //------------------------------------------------------------------------
2611 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2612 {
2613   //---------------------------------------------
2614   // unregister track from the list
2615   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2616     Int_t index = track->GetClusterIndex(icluster);
2617     Int_t l=(index & 0xf0000000) >> 28;
2618     Int_t c=(index & 0x0fffffff) >> 00;
2619     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2620     for (Int_t itrack=0;itrack<4;itrack++){
2621       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2622         fgLayers[l].SetClusterTracks(itrack,c,-1);
2623       }
2624     }
2625   }
2626 }
2627 //------------------------------------------------------------------------
2628 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2629 {
2630   //-------------------------------------------------------------
2631   //get number of shared clusters
2632   //-------------------------------------------------------------
2633   Float_t shared=0;
2634   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2635   // mean  number of clusters
2636   Float_t *ny = GetNy(id), *nz = GetNz(id);
2637
2638
2639   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2640     Int_t index = track->GetClusterIndex(icluster);
2641     Int_t l=(index & 0xf0000000) >> 28;
2642     Int_t c=(index & 0x0fffffff) >> 00;
2643     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2644     if (ny[l]==0){
2645       printf("problem\n");
2646     }
2647     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2648     Float_t weight=1;
2649     //
2650     Float_t deltan = 0;
2651     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2652     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2653       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2654     if (l<2 || l>3){
2655       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2656     }
2657     else{
2658       deltan = (cl->GetNz()-nz[l]);
2659     }
2660     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2661     weight = 2./TMath::Max(3.+deltan,2.);
2662     //
2663     for (Int_t itrack=0;itrack<4;itrack++){
2664       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2665         list[l]=index;
2666         clist[l] = (AliITSRecPoint*)GetCluster(index);
2667         shared+=weight;
2668         break;
2669       }
2670     }
2671   }
2672   track->SetNUsed(shared);
2673   return shared;
2674 }
2675 //------------------------------------------------------------------------
2676 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2677 {
2678   //
2679   // find first shared track
2680   //
2681   // mean  number of clusters
2682   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2683   //
2684   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2685   Int_t sharedtrack=100000;
2686   Int_t tracks[24],trackindex=0;
2687   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2688   //
2689   for (Int_t icluster=0;icluster<6;icluster++){
2690     if (clusterlist[icluster]<0) continue;
2691     Int_t index = clusterlist[icluster];
2692     Int_t l=(index & 0xf0000000) >> 28;
2693     Int_t c=(index & 0x0fffffff) >> 00;
2694     if (ny[l]==0){
2695       printf("problem\n");
2696     }
2697     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2698     //if (l>3) continue;
2699     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2700     //
2701     Float_t deltan = 0;
2702     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2703     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2704       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2705     if (l<2 || l>3){
2706       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2707     }
2708     else{
2709       deltan = (cl->GetNz()-nz[l]);
2710     }
2711     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2712     //
2713     for (Int_t itrack=3;itrack>=0;itrack--){
2714       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2715       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2716        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2717        trackindex++;
2718       }
2719     }
2720   }
2721   if (trackindex==0) return -1;
2722   if (trackindex==1){
2723     sharedtrack = tracks[0];
2724   }else{
2725     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2726     else{
2727       //
2728       Int_t tracks2[24], cluster[24];
2729       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2730       Int_t index =0;
2731       //
2732       for (Int_t i=0;i<trackindex;i++){
2733         if (tracks[i]<0) continue;
2734         tracks2[index] = tracks[i];
2735         cluster[index]++;
2736         for (Int_t j=i+1;j<trackindex;j++){
2737           if (tracks[j]<0) continue;
2738           if (tracks[j]==tracks[i]){
2739             cluster[index]++;
2740             tracks[j]=-1;
2741           }
2742         }
2743         index++;
2744       }
2745       Int_t max=0;
2746       for (Int_t i=0;i<index;i++){
2747         if (cluster[index]>max) {
2748           sharedtrack=tracks2[index];
2749           max=cluster[index];
2750         }
2751       }
2752     }
2753   }
2754
2755   if (sharedtrack>=100000) return -1;
2756   //
2757   // make list of overlaps
2758   shared =0;
2759   for (Int_t icluster=0;icluster<6;icluster++){
2760     if (clusterlist[icluster]<0) continue;
2761     Int_t index = clusterlist[icluster];
2762     Int_t l=(index & 0xf0000000) >> 28;
2763     Int_t c=(index & 0x0fffffff) >> 00;
2764     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2765     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2766     if (l==0 || l==1){
2767       if (cl->GetNy()>2) continue;
2768       if (cl->GetNz()>2) continue;
2769     }
2770     if (l==4 || l==5){
2771       if (cl->GetNy()>3) continue;
2772       if (cl->GetNz()>3) continue;
2773     }
2774     //
2775     for (Int_t itrack=3;itrack>=0;itrack--){
2776       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2777       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2778         overlist[l]=index;
2779         shared++;
2780       }
2781     }
2782   }
2783   return sharedtrack;
2784 }
2785 //------------------------------------------------------------------------
2786 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2787   //
2788   // try to find track hypothesys without conflicts
2789   // with minimal chi2;
2790   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2791   Int_t entries1 = arr1->GetEntriesFast();
2792   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2793   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2794   Int_t entries2 = arr2->GetEntriesFast();
2795   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2796   //
2797   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2798   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2799   if (track10->Pt()>0.5+track20->Pt()) return track10;
2800
2801   for (Int_t itrack=0;itrack<entries1;itrack++){
2802     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2803     UnRegisterClusterTracks(track,trackID1);
2804   }
2805   //
2806   for (Int_t itrack=0;itrack<entries2;itrack++){
2807     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2808     UnRegisterClusterTracks(track,trackID2);
2809   }
2810   Int_t index1=0;
2811   Int_t index2=0;
2812   Float_t maxconflicts=6;
2813   Double_t maxchi2 =1000.;
2814   //
2815   // get weight of hypothesys - using best hypothesy
2816   Double_t w1,w2;
2817
2818   Int_t list1[6],list2[6];
2819   AliITSRecPoint *clist1[6], *clist2[6] ;
2820   RegisterClusterTracks(track10,trackID1);
2821   RegisterClusterTracks(track20,trackID2);
2822   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2823   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2824   UnRegisterClusterTracks(track10,trackID1);
2825   UnRegisterClusterTracks(track20,trackID2);
2826   //
2827   // normalized chi2
2828   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2829   Float_t nerry[6],nerrz[6];
2830   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2831   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2832   for (Int_t i=0;i<6;i++){
2833      if ( (erry1[i]>0) && (erry2[i]>0)) {
2834        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2835        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2836      }else{
2837        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2838        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2839      }
2840      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2841        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2842        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2843        ncl1++;
2844      }
2845      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2846        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2847        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2848        ncl2++;
2849      }
2850   }
2851   chi21/=ncl1;
2852   chi22/=ncl2;
2853   //
2854   //
2855   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2856   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2857   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2858   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2859   //
2860   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2861         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2862         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2863         );
2864   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2865         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2866         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2867         );
2868
2869   Double_t sumw = w1+w2;
2870   w1/=sumw;
2871   w2/=sumw;
2872   if (w1<w2*0.5) {
2873     w1 = (d2+0.5)/(d1+d2+1.);
2874     w2 = (d1+0.5)/(d1+d2+1.);
2875   }
2876   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2877   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2878   //
2879   // get pair of "best" hypothesys
2880   //
2881   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2882   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2883
2884   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2885     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2886     //if (track1->fFakeRatio>0) continue;
2887     RegisterClusterTracks(track1,trackID1);
2888     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2889       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2890
2891       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2892       //if (track2->fFakeRatio>0) continue;
2893       Float_t nskipped=0;
2894       RegisterClusterTracks(track2,trackID2);
2895       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2896       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2897       UnRegisterClusterTracks(track2,trackID2);
2898       //
2899       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2900       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2901       if (nskipped>0.5) continue;
2902       //
2903       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2904       if (conflict1+1<cconflict1) continue;
2905       if (conflict2+1<cconflict2) continue;
2906       Float_t conflict=0;
2907       Float_t sumchi2=0;
2908       Float_t sum=0;
2909       for (Int_t i=0;i<6;i++){
2910         //
2911         Float_t c1 =0.; // conflict coeficients
2912         Float_t c2 =0.;
2913         if (clist1[i]&&clist2[i]){
2914           Float_t deltan = 0;
2915           if (i<2 || i>3){
2916             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2917           }
2918           else{
2919             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2920           }
2921           c1  = 2./TMath::Max(3.+deltan,2.);
2922           c2  = 2./TMath::Max(3.+deltan,2.);
2923         }
2924         else{
2925           if (clist1[i]){
2926             Float_t deltan = 0;
2927             if (i<2 || i>3){
2928               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2929             }
2930             else{
2931               deltan = (clist1[i]->GetNz()-nz1[i]);
2932             }
2933             c1  = 2./TMath::Max(3.+deltan,2.);
2934             c2  = 0;
2935           }
2936
2937           if (clist2[i]){
2938             Float_t deltan = 0;
2939             if (i<2 || i>3){
2940               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2941             }
2942             else{
2943               deltan = (clist2[i]->GetNz()-nz2[i]);
2944             }
2945             c2  = 2./TMath::Max(3.+deltan,2.);
2946             c1  = 0;
2947           }
2948         }
2949         //
2950         chi21=0;chi22=0;
2951         if (TMath::Abs(track1->GetDy(i))>0.) {
2952           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2953             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2954           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2955           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2956         }else{
2957           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2958         }
2959         //
2960         if (TMath::Abs(track2->GetDy(i))>0.) {
2961           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2962             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2963           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2964           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2965         }
2966         else{
2967           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2968         }
2969         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2970         if (chi21>0) sum+=w1;
2971         if (chi22>0) sum+=w2;
2972         conflict+=(c1+c2);
2973       }
2974       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2975       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2976       Double_t normchi2 = 2*conflict+sumchi2/sum;
2977       if ( normchi2 <maxchi2 ){
2978         index1 = itrack1;
2979         index2 = itrack2;
2980         maxconflicts = conflict;
2981         maxchi2 = normchi2;
2982       }
2983     }
2984     UnRegisterClusterTracks(track1,trackID1);
2985   }
2986   //
2987   //  if (maxconflicts<4 && maxchi2<th0){
2988   if (maxchi2<th0*2.){
2989     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
2990     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2991     track1->SetChi2MIP(5,maxconflicts);
2992     track1->SetChi2MIP(6,maxchi2);
2993     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
2994     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
2995     track1->SetChi2MIP(8,index1);
2996     fBestTrackIndex[trackID1] =index1;
2997     UpdateESDtrack(track1, AliESDtrack::kITSin);
2998   }
2999   else if (track10->GetChi2MIP(0)<th1){
3000     track10->SetChi2MIP(5,maxconflicts);
3001     track10->SetChi2MIP(6,maxchi2);
3002     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3003     UpdateESDtrack(track10,AliESDtrack::kITSin);
3004   }
3005
3006   for (Int_t itrack=0;itrack<entries1;itrack++){
3007     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3008     UnRegisterClusterTracks(track,trackID1);
3009   }
3010   //
3011   for (Int_t itrack=0;itrack<entries2;itrack++){
3012     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3013     UnRegisterClusterTracks(track,trackID2);
3014   }
3015
3016   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3017       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3018     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3019   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3020     RegisterClusterTracks(track10,trackID1);
3021   }
3022   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3023       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3024     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3025     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3026     RegisterClusterTracks(track20,trackID2);
3027   }
3028   return track10;
3029
3030 }
3031 //------------------------------------------------------------------------
3032 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3033   //--------------------------------------------------------------------
3034   // This function marks clusters assigned to the track
3035   //--------------------------------------------------------------------
3036   AliTracker::UseClusters(t,from);
3037
3038   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3039   //if (c->GetQ()>2) c->Use();
3040   if (c->GetSigmaZ2()>0.1) c->Use();
3041   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3042   //if (c->GetQ()>2) c->Use();
3043   if (c->GetSigmaZ2()>0.1) c->Use();
3044
3045 }
3046 //------------------------------------------------------------------------
3047 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3048 {
3049   //------------------------------------------------------------------
3050   // add track to the list of hypothesys
3051   //------------------------------------------------------------------
3052
3053   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3054     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3055   //
3056   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3057   if (!array) {
3058     array = new TObjArray(10);
3059     fTrackHypothesys.AddAt(array,esdindex);
3060   }
3061   array->AddLast(track);
3062 }
3063 //------------------------------------------------------------------------
3064 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3065 {
3066   //-------------------------------------------------------------------
3067   // compress array of track hypothesys
3068   // keep only maxsize best hypothesys
3069   //-------------------------------------------------------------------
3070   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3071   if (! (fTrackHypothesys.At(esdindex)) ) return;
3072   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3073   Int_t entries = array->GetEntriesFast();
3074   //
3075   //- find preliminary besttrack as a reference
3076   Float_t minchi2=10000;
3077   Int_t maxn=0;
3078   AliITStrackMI * besttrack=0;
3079   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3080     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3081     if (!track) continue;
3082     Float_t chi2 = NormalizedChi2(track,0);
3083     //
3084     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3085     track->SetLabel(tpcLabel);
3086     CookdEdx(track);
3087     track->SetFakeRatio(1.);
3088     CookLabel(track,0.); //For comparison only
3089     //
3090     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3091     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3092       if (track->GetNumberOfClusters()<maxn) continue;
3093       maxn = track->GetNumberOfClusters();
3094       if (chi2<minchi2){
3095         minchi2=chi2;
3096         besttrack=track;
3097       }
3098     }
3099     else{
3100       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3101         delete array->RemoveAt(itrack);
3102       }
3103     }
3104   }
3105   if (!besttrack) return;
3106   //
3107   //
3108   //take errors of best track as a reference
3109   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3110   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3111   for (Int_t j=0;j<6;j++) {
3112     if (besttrack->GetClIndex(j)>=0){
3113       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3114       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3115       ny[j]   = besttrack->GetNy(j);
3116       nz[j]   = besttrack->GetNz(j);
3117     }
3118   }
3119   //
3120   // calculate normalized chi2
3121   //
3122   Float_t * chi2        = new Float_t[entries];
3123   Int_t * index         = new Int_t[entries];
3124   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3125   for (Int_t itrack=0;itrack<entries;itrack++){
3126     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3127     if (track){
3128       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3129       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3130         chi2[itrack] = track->GetChi2MIP(0);
3131       else{
3132         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3133           delete array->RemoveAt(itrack);
3134         }
3135       }
3136     }
3137   }
3138   //
3139   TMath::Sort(entries,chi2,index,kFALSE);
3140   besttrack = (AliITStrackMI*)array->At(index[0]);
3141   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3142     for (Int_t j=0;j<6;j++){
3143       if (besttrack->GetClIndex(j)>=0){
3144         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3145         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3146         ny[j]   = besttrack->GetNy(j);
3147         nz[j]   = besttrack->GetNz(j);
3148       }
3149     }
3150   }
3151   //
3152   // calculate one more time with updated normalized errors
3153   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3154   for (Int_t itrack=0;itrack<entries;itrack++){
3155     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3156     if (track){
3157       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3158       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3159         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3160       else
3161         {
3162           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3163             delete array->RemoveAt(itrack);
3164           }
3165         }
3166     }
3167   }
3168   entries = array->GetEntriesFast();
3169   //
3170   //
3171   if (entries>0){
3172     TObjArray * newarray = new TObjArray();
3173     TMath::Sort(entries,chi2,index,kFALSE);
3174     besttrack = (AliITStrackMI*)array->At(index[0]);
3175     if (besttrack){
3176       //
3177       for (Int_t j=0;j<6;j++){
3178         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3179           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3180           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3181           ny[j]   = besttrack->GetNy(j);
3182           nz[j]   = besttrack->GetNz(j);
3183         }
3184       }
3185       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3186       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3187       Float_t minn = besttrack->GetNumberOfClusters()-3;
3188       Int_t accepted=0;
3189       for (Int_t i=0;i<entries;i++){
3190         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3191         if (!track) continue;
3192         if (accepted>maxcut) break;
3193         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3194         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3195           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3196             delete array->RemoveAt(index[i]);
3197             continue;
3198           }
3199         }
3200         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3201         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3202           if (!shortbest) accepted++;
3203           //
3204           newarray->AddLast(array->RemoveAt(index[i]));
3205           for (Int_t j=0;j<6;j++){
3206             if (nz[j]==0){
3207               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3208               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3209               ny[j]   = track->GetNy(j);
3210               nz[j]   = track->GetNz(j);
3211             }
3212           }
3213         }
3214         else{
3215           delete array->RemoveAt(index[i]);
3216         }
3217       }
3218       array->Delete();
3219       delete fTrackHypothesys.RemoveAt(esdindex);
3220       fTrackHypothesys.AddAt(newarray,esdindex);
3221     }
3222     else{
3223       array->Delete();
3224       delete fTrackHypothesys.RemoveAt(esdindex);
3225     }
3226   }
3227   delete [] chi2;
3228   delete [] index;
3229 }
3230 //------------------------------------------------------------------------
3231 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3232 {
3233   //-------------------------------------------------------------
3234   // try to find best hypothesy
3235   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3236   //-------------------------------------------------------------
3237   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3238   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3239   if (!array) return 0;
3240   Int_t entries = array->GetEntriesFast();
3241   if (!entries) return 0;
3242   Float_t minchi2 = 100000;
3243   AliITStrackMI * besttrack=0;
3244   //
3245   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3246   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3247   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3248   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3249   //
3250   for (Int_t i=0;i<entries;i++){
3251     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3252     if (!track) continue;
3253     Float_t sigmarfi,sigmaz;
3254     GetDCASigma(track,sigmarfi,sigmaz);
3255     track->SetDnorm(0,sigmarfi);
3256     track->SetDnorm(1,sigmaz);
3257     //
3258     track->SetChi2MIP(1,1000000);
3259     track->SetChi2MIP(2,1000000);
3260     track->SetChi2MIP(3,1000000);
3261     //
3262     // backtrack
3263     backtrack = new(backtrack) AliITStrackMI(*track);
3264     if (track->GetConstrain()) {
3265       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3266       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;
3267       backtrack->ResetCovariance(10.);
3268     }else{
3269       backtrack->ResetCovariance(10.);
3270     }
3271     backtrack->ResetClusters();
3272
3273     Double_t x = original->GetX();
3274     if (!RefitAt(x,backtrack,track)) continue;
3275     //
3276     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3277     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3278     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3279     track->SetChi22(GetMatchingChi2(backtrack,original));
3280
3281     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3282     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3283     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3284
3285
3286     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3287     //
3288     //forward track - without constraint
3289     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3290     forwardtrack->ResetClusters();
3291     x = track->GetX();
3292     RefitAt(x,forwardtrack,track);
3293     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3294     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3295     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3296
3297     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3298     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3299     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3300     forwardtrack->SetD(0,track->GetD(0));
3301     forwardtrack->SetD(1,track->GetD(1));
3302     {
3303       Int_t list[6];
3304       AliITSRecPoint* clist[6];
3305       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3306       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3307     }
3308
3309     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3310     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3311     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3312       track->SetChi2MIP(3,1000);
3313       continue;
3314     }
3315     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3316     //
3317     for (Int_t ichi=0;ichi<5;ichi++){
3318       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3319     }
3320     if (chi2 < minchi2){
3321       //besttrack = new AliITStrackMI(*forwardtrack);
3322       besttrack = track;
3323       besttrack->SetLabel(track->GetLabel());
3324       besttrack->SetFakeRatio(track->GetFakeRatio());
3325       minchi2   = chi2;
3326       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3327       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3328       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3329     }
3330   }
3331   delete backtrack;
3332   delete forwardtrack;
3333   Int_t accepted=0;
3334   for (Int_t i=0;i<entries;i++){
3335     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3336
3337     if (!track) continue;
3338
3339     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3340         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3341         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3342       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3343         delete array->RemoveAt(i);
3344         continue;
3345       }
3346     }
3347     else{
3348       accepted++;
3349     }
3350   }
3351   //
3352   array->Compress();
3353   SortTrackHypothesys(esdindex,checkmax,1);
3354
3355   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3356   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3357   besttrack = (AliITStrackMI*)array->At(0);
3358   if (!besttrack)  return 0;
3359   besttrack->SetChi2MIP(8,0);
3360   fBestTrackIndex[esdindex]=0;
3361   entries = array->GetEntriesFast();
3362   AliITStrackMI *longtrack =0;
3363   minchi2 =1000;
3364   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3365   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3366     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3367     if (!track->GetConstrain()) continue;
3368     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3369     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3370     if (track->GetChi2MIP(0)>4.) continue;
3371     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3372     longtrack =track;
3373   }
3374   //if (longtrack) besttrack=longtrack;
3375
3376   Int_t list[6];
3377   AliITSRecPoint * clist[6];
3378   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3379   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3380       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3381     RegisterClusterTracks(besttrack,esdindex);
3382   }
3383   //
3384   //
3385   if (shared>0.0){
3386     Int_t nshared;
3387     Int_t overlist[6];
3388     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3389     if (sharedtrack>=0){
3390       //
3391       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3392       if (besttrack){
3393         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3394       }
3395       else return 0;
3396     }
3397   }
3398
3399   if (shared>2.5) return 0;
3400   if (shared>1.0) return besttrack;
3401   //
3402   // Don't sign clusters if not gold track
3403   //
3404   if (!besttrack->IsGoldPrimary()) return besttrack;
3405   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3406   //
3407   if (fConstraint[fPass]){
3408     //
3409     // sign clusters
3410     //
3411     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3412     for (Int_t i=0;i<6;i++){
3413       Int_t index = besttrack->GetClIndex(i);
3414       if (index<0) continue;
3415       Int_t ilayer =  (index & 0xf0000000) >> 28;
3416       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3417       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3418       if (!c) continue;
3419       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3420       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3421       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3422       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3423         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3424       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3425
3426       Bool_t cansign = kTRUE;
3427       for (Int_t itrack=0;itrack<entries; itrack++){
3428         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3429         if (!track) continue;
3430         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3431         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3432           cansign = kFALSE;
3433           break;
3434         }
3435       }
3436       if (cansign){
3437         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3438         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3439         if (!c->IsUsed()) c->Use();
3440       }
3441     }
3442   }
3443   return besttrack;
3444 }
3445 //------------------------------------------------------------------------
3446 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3447 {
3448   //
3449   // get "best" hypothesys
3450   //
3451
3452   Int_t nentries = itsTracks.GetEntriesFast();
3453   for (Int_t i=0;i<nentries;i++){
3454     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3455     if (!track) continue;
3456     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3457     if (!array) continue;
3458     if (array->GetEntriesFast()<=0) continue;
3459     //
3460     AliITStrackMI* longtrack=0;
3461     Float_t minn=0;
3462     Float_t maxchi2=1000;
3463     for (Int_t j=0;j<array->GetEntriesFast();j++){
3464       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3465       if (!trackHyp) continue;
3466       if (trackHyp->GetGoldV0()) {
3467         longtrack = trackHyp;   //gold V0 track taken
3468         break;
3469       }
3470       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3471       Float_t chi2 = trackHyp->GetChi2MIP(0);
3472       if (fAfterV0){
3473         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3474       }
3475       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3476       //
3477       if (chi2 > maxchi2) continue;
3478       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3479       maxchi2 = chi2;
3480       longtrack=trackHyp;
3481     }
3482     //
3483     //
3484     //
3485     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3486     if (!longtrack) {longtrack = besttrack;}
3487     else besttrack= longtrack;
3488     //
3489     if (besttrack) {
3490       Int_t list[6];
3491       AliITSRecPoint * clist[6];
3492       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3493       //
3494       track->SetNUsed(shared);
3495       track->SetNSkipped(besttrack->GetNSkipped());
3496       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3497       if (shared>0) {
3498         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3499         Int_t nshared;
3500         Int_t overlist[6];
3501         //
3502         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3503         //if (sharedtrack==-1) sharedtrack=0;
3504         if (sharedtrack>=0) {
3505           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3506         }
3507       }
3508       if (besttrack&&fAfterV0) {
3509         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3510       }
3511       if (besttrack&&fConstraint[fPass])
3512         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3513       if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3514         if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3515              TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3516       }
3517
3518     }
3519   }
3520 }
3521 //------------------------------------------------------------------------
3522 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3523   //--------------------------------------------------------------------
3524   //This function "cooks" a track label. If label<0, this track is fake.
3525   //--------------------------------------------------------------------
3526   Int_t tpcLabel=-1;
3527
3528   if ( track->GetESDtrack())   tpcLabel =  TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3529
3530    track->SetChi2MIP(9,0);
3531    Int_t nwrong=0;
3532    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3533      Int_t cindex = track->GetClusterIndex(i);
3534      Int_t l=(cindex & 0xf0000000) >> 28;
3535      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3536      Int_t isWrong=1;
3537      for (Int_t ind=0;ind<3;ind++){
3538        if (tpcLabel>0)
3539          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3540        AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3541      }
3542      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3543      nwrong+=isWrong;
3544    }
3545    Int_t nclusters = track->GetNumberOfClusters();
3546    if (nclusters > 0) //PH Some tracks don't have any cluster
3547      track->SetFakeRatio(double(nwrong)/double(nclusters));
3548    if (tpcLabel>0){
3549      if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3550      else
3551        track->SetLabel(tpcLabel);
3552    }
3553    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3554
3555 }
3556 //------------------------------------------------------------------------
3557 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3558 {
3559   //
3560   //
3561   //  Int_t list[6];
3562   //AliITSRecPoint * clist[6];
3563   //  Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3564   Float_t dedx[4];
3565   Int_t accepted=0;
3566   track->SetChi2MIP(9,0);
3567   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3568     Int_t cindex = track->GetClusterIndex(i);
3569     Int_t l=(cindex & 0xf0000000) >> 28;
3570     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3571     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3572     Int_t isWrong=1;
3573     for (Int_t ind=0;ind<3;ind++){
3574       if (cl->GetLabel(ind)==lab) isWrong=0;
3575     }
3576     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3577     if (l<2) continue;
3578     //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue;  //shared track
3579     //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3580     //if (l<4&& !(cl->GetType()==1)) continue;
3581     dedx[accepted]= track->GetSampledEdx(i);
3582     //dedx[accepted]= track->fNormQ[l];
3583     accepted++;
3584   }
3585   if (accepted<1) {
3586     track->SetdEdx(0);
3587     return;
3588   }
3589   Int_t indexes[4];
3590   TMath::Sort(accepted,dedx,indexes,kFALSE);
3591   Double_t low=0.;
3592   Double_t up=0.51;
3593   Double_t nl=low*accepted, nu =up*accepted;
3594   Float_t sumamp = 0;
3595   Float_t sumweight =0;
3596   for (Int_t i=0; i<accepted; i++) {
3597     Float_t weight =1;
3598     if (i<nl+0.1)     weight = TMath::Max(1.-(nl-i),0.);
3599     if (i>nu-1)     weight = TMath::Max(nu-i,0.);
3600     sumamp+= dedx[indexes[i]]*weight;
3601     sumweight+=weight;
3602   }
3603   track->SetdEdx(sumamp/sumweight);
3604 }
3605 //------------------------------------------------------------------------
3606 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3607   //
3608   // Create some arrays
3609   //
3610   if (fCoefficients) delete []fCoefficients;
3611   fCoefficients = new Float_t[ntracks*48];
3612   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3613 }
3614 //------------------------------------------------------------------------
3615 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3616 {
3617   //
3618   // Compute predicted chi2
3619   //
3620   Float_t erry,errz;
3621   Float_t theta = track->GetTgl();
3622   Float_t phi   = track->GetSnp();
3623   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3624   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3625   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3626   // Take into account the mis-alignment (bring track to cluster plane)
3627   Double_t xTrOrig=track->GetX();
3628   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3629   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3630   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3631   // Bring the track back to detector plane in ideal geometry
3632   // [mis-alignment will be accounted for in UpdateMI()]
3633   if (!track->Propagate(xTrOrig)) return 1000.;
3634   Float_t ny,nz;
3635   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3636   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3637   if (delta>1){
3638     chi2+=0.5*TMath::Min(delta/2,2.);
3639     chi2+=2.*cluster->GetDeltaProbability();
3640   }
3641   //
3642   track->SetNy(layer,ny);
3643   track->SetNz(layer,nz);
3644   track->SetSigmaY(layer,erry);
3645   track->SetSigmaZ(layer, errz);
3646   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3647   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3648   return chi2;
3649
3650 }
3651 //------------------------------------------------------------------------
3652 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3653 {
3654   //
3655   // Update ITS track
3656   //
3657   Int_t layer = (index & 0xf0000000) >> 28;
3658   track->SetClIndex(layer, index);
3659   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3660     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3661       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3662       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3663     }
3664   }
3665
3666   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
3667
3668
3669   // Take into account the mis-alignment (bring track to cluster plane)
3670   Double_t xTrOrig=track->GetX();
3671   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3672   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3673   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3674   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
3675
3676   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3677
3678
3679   AliCluster c(*cl);
3680   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3681   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3682
3683
3684   Int_t updated = track->UpdateMI(&c,chi2,index);
3685
3686   // Bring the track back to detector plane in ideal geometry
3687   if (!track->Propagate(xTrOrig)) return 0;
3688
3689   if(!updated) AliDebug(2,"update failed");
3690   return updated;
3691 }
3692
3693 //------------------------------------------------------------------------
3694 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3695 {
3696   //
3697   //DCA sigmas parameterization
3698   //to be paramterized using external parameters in future
3699   //
3700   //
3701   sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3702   sigmaz   = 0.0110+4.37*TMath::Abs(track->GetC());
3703 }
3704 //------------------------------------------------------------------------
3705 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3706 {
3707   //
3708   // Clusters from delta electrons?
3709   //
3710   Int_t entries = clusterArray->GetEntriesFast();
3711   if (entries<4) return;
3712   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3713   Int_t layer = cluster->GetLayer();
3714   if (layer>1) return;
3715   Int_t index[10000];
3716   Int_t ncandidates=0;
3717   Float_t r = (layer>0)? 7:4;
3718   //
3719   for (Int_t i=0;i<entries;i++){
3720     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3721     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3722     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3723     index[ncandidates] = i;  //candidate to belong to delta electron track
3724     ncandidates++;
3725     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3726       cl0->SetDeltaProbability(1);
3727     }
3728   }
3729   //
3730   //
3731   //
3732   for (Int_t i=0;i<ncandidates;i++){
3733     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3734     if (cl0->GetDeltaProbability()>0.8) continue;
3735     //
3736     Int_t ncl = 0;
3737     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3738     sumy=sumz=sumy2=sumyz=sumw=0.0;
3739     for (Int_t j=0;j<ncandidates;j++){
3740       if (i==j) continue;
3741       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3742       //
3743       Float_t dz = cl0->GetZ()-cl1->GetZ();
3744       Float_t dy = cl0->GetY()-cl1->GetY();
3745       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3746         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3747         y[ncl] = cl1->GetY();
3748         z[ncl] = cl1->GetZ();
3749         sumy+= y[ncl]*weight;
3750         sumz+= z[ncl]*weight;
3751         sumy2+=y[ncl]*y[ncl]*weight;
3752         sumyz+=y[ncl]*z[ncl]*weight;
3753         sumw+=weight;
3754         ncl++;
3755       }
3756     }
3757     if (ncl<4) continue;
3758     Float_t det = sumw*sumy2  - sumy*sumy;
3759     Float_t delta=1000;
3760     if (TMath::Abs(det)>0.01){
3761       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3762       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3763       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3764     }
3765     else{
3766       Float_t z0  = sumyz/sumy;
3767       delta = TMath::Abs(cl0->GetZ()-z0);
3768     }
3769     if ( delta<0.05) {
3770       cl0->SetDeltaProbability(1-20.*delta);
3771     }
3772   }
3773 }
3774 //------------------------------------------------------------------------
3775 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3776 {
3777   //
3778   // Update ESD track
3779   //
3780   track->UpdateESDtrack(flags);
3781   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3782   if (oldtrack) delete oldtrack;
3783   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3784   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3785     printf("Problem\n");
3786   }
3787 }
3788 //------------------------------------------------------------------------
3789 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3790   //
3791   // Get nearest upper layer close to the point xr.
3792   // rough approximation
3793   //
3794   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3795   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3796   Int_t res =6;
3797   for (Int_t i=0;i<6;i++){
3798     if (radius<kRadiuses[i]){
3799       res =i;
3800       break;
3801     }
3802   }
3803   return res;
3804 }
3805
3806 //------------------------------------------------------------------------
3807 void AliITStrackerMI::UpdateTPCV0(const AliESDEvent *event){
3808   //
3809   //try to update, or reject TPC  V0s
3810   //
3811   Int_t nv0s = event->GetNumberOfV0s();
3812   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3813
3814   for (Int_t i=0;i<nv0s;i++){
3815     AliESDv0 * vertex = event->GetV0(i);
3816     Int_t ip = vertex->GetIndex(0);
3817     Int_t im = vertex->GetIndex(1);
3818     //
3819     TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3820     TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3821     AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3822     AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3823     //
3824     //
3825     if (trackp){
3826       if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3827         if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3828         if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3829       }
3830     }
3831
3832     if (trackm){
3833       if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3834         if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3835         if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3836       }
3837     }
3838     if (vertex->GetStatus()==-100) continue;
3839     //
3840     Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
3841     Int_t clayer = GetNearestLayer(xrp);                    //I.B.
3842     vertex->SetNBefore(clayer);        //
3843     vertex->SetChi2Before(9*clayer);   //
3844     vertex->SetNAfter(6-clayer);       //
3845     vertex->SetChi2After(0);           //
3846     //
3847     if (clayer >1 ){ // calculate chi2 before vertex
3848       Float_t chi2p = 0, chi2m=0;
3849       //
3850       if (trackp){
3851         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3852           if (trackp->GetClIndex(ilayer)>=0){
3853             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3854               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3855           }
3856           else{
3857             chi2p+=9;
3858           }
3859         }
3860       }else{
3861         chi2p = 9*clayer;
3862       }
3863       //
3864       if (trackm){
3865         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3866           if (trackm->GetClIndex(ilayer)>=0){
3867             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3868               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3869           }
3870           else{
3871             chi2m+=9;
3872           }
3873         }
3874       }else{
3875         chi2m = 9*clayer;
3876       }
3877       vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3878       if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10);  // track exist before vertex
3879     }
3880
3881     if (clayer < 5 ){ // calculate chi2 after vertex
3882       Float_t chi2p = 0, chi2m=0;
3883       //
3884       if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3885         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3886           if (trackp->GetClIndex(ilayer)>=0){
3887             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3888               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3889           }
3890           else{
3891             chi2p+=9;
3892           }
3893         }
3894       }else{
3895         chi2p = 0;
3896       }
3897       //
3898       if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3899         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3900           if (trackm->GetClIndex(ilayer)>=0){
3901             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3902               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3903           }
3904           else{
3905             chi2m+=9;
3906           }
3907         }
3908       }else{
3909         chi2m = 0;
3910       }
3911       vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3912       if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20);  // track not found in ITS
3913     }
3914   }
3915   //
3916 }
3917 //------------------------------------------------------------------------
3918 void AliITStrackerMI::FindV02(AliESDEvent *event)
3919 {
3920   //
3921   // V0 finder
3922   //
3923   //  Cuts on DCA -  R dependent
3924   //          max distance DCA between 2 tracks cut
3925   //          maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3926   //
3927   const Float_t kMaxDist0      = 0.1;
3928   const Float_t kMaxDist1      = 0.1;
3929   const Float_t kMaxDist       = 1;
3930   const Float_t kMinPointAngle  = 0.85;
3931   const Float_t kMinPointAngle2 = 0.99;
3932   const Float_t kMinR           = 0.5;
3933   const Float_t kMaxR           = 220;
3934   //const Float_t kCausality0Cut   = 0.19;
3935   //const Float_t kLikelihood01Cut = 0.25;
3936   //const Float_t kPointAngleCut   = 0.9996;
3937   const Float_t kCausality0Cut   = 0.19;
3938   const Float_t kLikelihood01Cut = 0.45;
3939   const Float_t kLikelihood1Cut  = 0.5;
3940   const Float_t kCombinedCut     = 0.55;
3941
3942   //
3943   //
3944   TTreeSRedirector &cstream = *fDebugStreamer;
3945   Int_t ntracks    = event->GetNumberOfTracks();
3946   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3947   fOriginal.Expand(ntracks);
3948   fTrackHypothesys.Expand(ntracks);
3949   fBestHypothesys.Expand(ntracks);
3950   //
3951   AliHelix * helixes   = new AliHelix[ntracks+2];
3952   TObjArray trackarray(ntracks+2);     //array with tracks - with vertex constrain
3953   TObjArray trackarrayc(ntracks+2);    //array of "best    tracks" - without vertex constrain
3954   TObjArray trackarrayl(ntracks+2);    //array of "longest tracks" - without vertex constrain
3955   Bool_t * forbidden   = new Bool_t [ntracks+2];
3956   Int_t   *itsmap      = new Int_t  [ntracks+2];
3957   Float_t *dist        = new Float_t[ntracks+2];
3958   Float_t *normdist0   = new Float_t[ntracks+2];
3959   Float_t *normdist1   = new Float_t[ntracks+2];
3960   Float_t *normdist    = new Float_t[ntracks+2];
3961   Float_t *norm        = new Float_t[ntracks+2];
3962   Float_t *maxr        = new Float_t[ntracks+2];
3963   Float_t *minr        = new Float_t[ntracks+2];
3964   Float_t *minPointAngle= new Float_t[ntracks+2];
3965   //
3966   AliV0 *pvertex      = new AliV0;
3967   AliITStrackMI * dummy= new AliITStrackMI;
3968   dummy->SetLabel(0);
3969   AliITStrackMI  trackat0;    //temporary track for DCA calculation
3970   //
3971   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3972   //
3973   // make ITS -  ESD map
3974   //
3975   for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3976     itsmap[itrack]        = -1;
3977     forbidden[itrack]     = kFALSE;
3978     maxr[itrack]          = kMaxR;
3979     minr[itrack]          = kMinR;
3980     minPointAngle[itrack] = kMinPointAngle;
3981   }
3982   for (Int_t itrack=0;itrack<nitstracks;itrack++){
3983     AliITStrackMI * original =   (AliITStrackMI*)(fOriginal.At(itrack));
3984     Int_t           esdindex =   original->GetESDtrack()->GetID();
3985     itsmap[esdindex]         =   itrack;
3986   }
3987   //
3988   // create ITS tracks from ESD tracks if not done before
3989   //
3990   for (Int_t itrack=0;itrack<ntracks;itrack++){
3991     if (itsmap[itrack]>=0) continue;
3992     AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3993     //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3994     //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3995     tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP());   //I.B.
3996     if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3997       // tracks which can reach inner part of ITS
3998       // propagate track to outer its volume - with correction for material
3999       CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
4000     }
4001     itsmap[itrack] = nitstracks;
4002     fOriginal.AddAt(tpctrack,nitstracks);
4003     nitstracks++;
4004   }
4005   //
4006   // fill temporary arrays
4007   //
4008   for (Int_t itrack=0;itrack<ntracks;itrack++){
4009     AliESDtrack *  esdtrack = event->GetTrack(itrack);
4010     Int_t          itsindex = itsmap[itrack];
4011     AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
4012     if (!original) continue;
4013     AliITStrackMI *bestConst  = 0;
4014     AliITStrackMI *bestLong   = 0;
4015     AliITStrackMI *best       = 0;
4016     //
4017     //
4018     TObjArray * array    = (TObjArray*)  fTrackHypothesys.At(itsindex);
4019     Int_t       hentries = (array==0) ?  0 : array->GetEntriesFast();
4020     // Get best track with vertex constrain
4021     for (Int_t ih=0;ih<hentries;ih++){
4022       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4023       if (!trackh->GetConstrain()) continue;
4024       if (!bestConst) bestConst = trackh;
4025       if (trackh->GetNumberOfClusters()>5.0){
4026         bestConst  = trackh;                         // full track -  with minimal chi2
4027         break;
4028       }
4029       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone())  continue;
4030       bestConst = trackh;
4031       break;
4032     }
4033     // Get best long track without vertex constrain and best track without vertex constrain
4034     for (Int_t ih=0;ih<hentries;ih++){
4035       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4036       if (trackh->GetConstrain()) continue;
4037       if (!best)     best     = trackh;
4038       if (!bestLong) bestLong = trackh;
4039       if (trackh->GetNumberOfClusters()>5.0){
4040         bestLong  = trackh;                         // full track -  with minimal chi2
4041         break;
4042       }
4043       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone())  continue;
4044       bestLong = trackh;
4045     }
4046     if (!best) {
4047       best     = original;
4048       bestLong = original;
4049     }
4050     //I.B. trackat0 = *bestLong;
4051     new (&trackat0) AliITStrackMI(*bestLong);
4052     Double_t xx,yy,zz,alpha;
4053     if (!bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz)) continue;
4054     alpha = TMath::ATan2(yy,xx);
4055     trackat0.Propagate(alpha,0); //PH The check on the return value is temporarily disabled (bug 45751)
4056     // calculate normalized distances to the vertex
4057     //
4058     Float_t ptfac  = (1.+100.*TMath::Abs(trackat0.GetC()));
4059     if ( bestLong->GetNumberOfClusters()>3 ){
4060       dist[itsindex]      = trackat0.GetY();
4061       norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4062       normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4063       normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4064       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4065       if (!bestConst){
4066         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
4067         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
4068         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
4069       }else{
4070         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
4071         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
4072       }
4073     }
4074     else{
4075       if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
4076         dist[itsindex] = bestConst->GetD(0);
4077         norm[itsindex] = bestConst->GetDnorm(0);
4078         normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4079         normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4080         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4081       }else{
4082         dist[itsindex]      = trackat0.GetY();
4083         norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4084         normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4085         normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4086         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4087         if (TMath::Abs(trackat0.GetTgl())>1.05){
4088           if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
4089           if (normdist[itsindex]>3) {
4090             minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
4091           }
4092         }
4093       }
4094     }
4095     //
4096     //-----------------------------------------------------------
4097     //Forbid primary track candidates -
4098     //
4099     //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
4100     //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
4101     //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
4102     //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
4103     //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
4104     //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
4105     //-----------------------------------------------------------
4106     if (bestConst){
4107       if (bestLong->GetNumberOfClusters()<4       && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5)               forbidden[itsindex]=kTRUE;
4108       if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5)               forbidden[itsindex]=kTRUE;
4109       if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>=0 && bestConst->GetClIndex(1)>=0 ) forbidden[itsindex]=kTRUE;
4110       if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>=0)                              forbidden[itsindex]=kTRUE;
4111       if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2)                             forbidden[itsindex]=kTRUE;
4112       if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1)                             forbidden[itsindex]=kTRUE;
4113       if (bestConst->GetNormChi2(0)<2.5) {
4114         minPointAngle[itsindex]= 0.9999;
4115         maxr[itsindex]         = 10;
4116       }
4117     }
4118     //
4119     //forbid daughter kink candidates
4120     //
4121     if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
4122     Bool_t isElectron = kTRUE;
4123     Bool_t isProton   = kTRUE;
4124     Double_t pid[5];
4125     esdtrack->GetESDpid(pid);
4126     for (Int_t i=1;i<5;i++){
4127       if (pid[0]<pid[i]) isElectron= kFALSE;
4128       if (pid[4]<pid[i]) isProton= kFALSE;
4129     }
4130     if (isElectron){
4131       forbidden[itsindex]=kFALSE;
4132       normdist[itsindex]*=-1;
4133     }
4134     if (isProton){
4135       if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
4136       normdist[itsindex]*=-1;
4137     }
4138
4139     //
4140     // Causality cuts in TPC volume
4141     //
4142     if (esdtrack->GetTPCdensity(0,10) >0.6)  maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
4143     if (esdtrack->GetTPCdensity(10,30)>0.6)  maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
4144     if (esdtrack->GetTPCdensity(20,40)>0.6)  maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
4145     if (esdtrack->GetTPCdensity(30,50)>0.6)  maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
4146     //
4147     if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
4148     //
4149     //
4150     if (kFALSE){
4151       cstream<<"Track"<<
4152         "Tr0.="<<best<<
4153         "Tr1.="<<((bestConst)? bestConst:dummy)<<
4154         "Tr2.="<<bestLong<<
4155         "Tr3.="<<&trackat0<<
4156         "Esd.="<<esdtrack<<
4157         "Dist="<<dist[itsindex]<<
4158         "ND0="<<normdist0[itsindex]<<
4159         "ND1="<<normdist1[itsindex]<<
4160         "ND="<<normdist[itsindex]<<
4161         "Pz="<<primvertex[2]<<
4162         "Forbid="<<forbidden[itsindex]<<
4163         "\n";
4164       //
4165     }
4166     trackarray.AddAt(best,itsindex);
4167     trackarrayc.AddAt(bestConst,itsindex);
4168     trackarrayl.AddAt(bestLong,itsindex);
4169     new (&helixes[itsindex]) AliHelix(*best);
4170   }
4171   //
4172   //
4173   //
4174   // first iterration of V0 finder
4175   //
4176   for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
4177     Int_t itrack0 = itsmap[iesd0];
4178     if (forbidden[itrack0]) continue;
4179     AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
4180     if (!btrack0) continue;
4181     if (btrack0->GetSign()>0 && !AliITSReconstructor::GetRecoParam()->GetStoreLikeSignV0s()) continue;
4182     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
4183     //
4184     for (Int_t iesd1=iesd0+1;iesd1<ntracks;iesd1++){
4185       Int_t itrack1 = itsmap[iesd1];
4186       if (forbidden[itrack1]) continue;
4187
4188       AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
4189       if (!btrack1) continue;
4190       if (btrack1->GetSign()<0 && !AliITSReconstructor::GetRecoParam()->GetStoreLikeSignV0s()) continue;
4191       Bool_t isGold = kFALSE;
4192       if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
4193         isGold = kTRUE;
4194       }
4195       AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
4196       AliHelix &h1 = helixes[itrack0];
4197       AliHelix &h2 = helixes[itrack1];
4198       //
4199       // find linear distance
4200       Double_t rmin =0;
4201       //
4202       //
4203       //
4204       Double_t phase[2][2],radius[2];
4205       Int_t  points = h1.GetRPHIintersections(h2, phase, radius);
4206       if    (points==0)  continue;
4207       Double_t delta[2]={1000000,1000000};
4208       rmin = radius[0];
4209       h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
4210       if (points==2){
4211         if (radius[1]<rmin) rmin = radius[1];
4212         h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
4213       }
4214       rmin = TMath::Sqrt(rmin);
4215       Double_t distance = 0;
4216       Double_t radiusC  = 0;
4217       Int_t    iphase   = 0;
4218       if (points==1 || delta[0]<delta[1]){
4219         distance = TMath::Sqrt(delta[0]);
4220         radiusC  = TMath::Sqrt(radius[0]);
4221       }else{
4222         distance = TMath::Sqrt(delta[1]);
4223         radiusC  = TMath::Sqrt(radius[1]);
4224         iphase=1;
4225       }
4226       if (radiusC<TMath::Max(minr[itrack0],minr[itrack1]))    continue;
4227       if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1]))     continue;
4228       Float_t maxDist  = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
4229       if (distance>maxDist) continue;
4230       Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
4231       if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
4232       //
4233       //
4234       //      Double_t distance = TestV0(h1,h2,pvertex,rmin);
4235       //
4236       //       if (distance>maxDist)           continue;
4237       //       if (pvertex->GetRr()<kMinR)     continue;
4238       //       if (pvertex->GetRr()>kMaxR)     continue;
4239       AliITStrackMI * track0=btrack0;
4240       AliITStrackMI * track1=btrack1;
4241       //      if (pvertex->GetRr()<3.5){
4242       if (radiusC<3.5){
4243         //use longest tracks inside the pipe
4244         track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4245         track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4246       }
4247       //
4248       //
4249       pvertex->SetParamN(*track0);
4250       pvertex->SetParamP(*track1);
4251       pvertex->Update(primvertex);
4252       pvertex->SetClusters(track0->ClIndex(),track1->ClIndex());  // register clusters
4253
4254       if (pvertex->GetRr()<kMinR) continue;
4255       if (pvertex->GetRr()>kMaxR) continue;
4256       if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4257 //Bo:      if (pvertex->GetDist2()>maxDist) continue;
4258       if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4259 //Bo:        pvertex->SetLab(0,track0->GetLabel());
4260 //Bo:        pvertex->SetLab(1,track1->GetLabel());
4261       pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4262       pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4263       //
4264       AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4265       AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4266
4267       //
4268       //
4269       TObjArray * array0b     = (TObjArray*)fBestHypothesys.At(itrack0);
4270       if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4271         fCurrentEsdTrack = itrack0;
4272         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4273       }
4274       TObjArray * array1b    = (TObjArray*)fBestHypothesys.At(itrack1);
4275       if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4276         fCurrentEsdTrack = itrack1;
4277         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4278       }
4279       //
4280       AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4281       AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4282       AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4283       AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4284
4285       Float_t minchi2before0=16;
4286       Float_t minchi2before1=16;
4287       Float_t minchi2after0 =16;
4288       Float_t minchi2after1 =16;
4289       Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
4290       Int_t maxLayer = GetNearestLayer(xrp);                   //I.B.
4291
4292       if (array0b) for (Int_t i=0;i<5;i++){
4293         // best track after vertex
4294         AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4295         if (!btrack) continue;
4296         if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4297         //      if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4298         if (btrack->GetX()<pvertex->GetRr()-2.) {
4299           if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4300             Float_t sumchi2= 0;
4301             Float_t sumn   = 0;
4302             if (maxLayer<3){   // take prim vertex as additional measurement
4303               if (normdist[itrack0]>0 && htrackc0){
4304                 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4305               }else{
4306                 sumchi2 +=  TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4307               }
4308               sumn    +=  3-maxLayer;
4309             }
4310             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4311               sumn+=1.;
4312               if (btrack->GetClIndex(ilayer)<0){
4313                 sumchi2+=25;
4314                 continue;
4315               }else{
4316                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4317                 for (Int_t itrack=0;itrack<4;itrack++){
4318                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4319                     sumchi2+=18.;  //shared cluster
4320                     break;
4321                   }
4322                 }
4323                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4324                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4325               }
4326             }
4327             sumchi2/=sumn;
4328             if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4329           }
4330           continue;   //safety space - Geo manager will give exact layer
4331         }
4332         track0b       = btrack;
4333         minchi2after0 = btrack->GetNormChi2(i);
4334         break;
4335       }
4336       if (array1b) for (Int_t i=0;i<5;i++){
4337         // best track after vertex
4338         AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4339         if (!btrack) continue;
4340         if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4341         //      if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4342         if (btrack->GetX()<pvertex->GetRr()-2){
4343           if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4344             Float_t sumchi2= 0;
4345             Float_t sumn   = 0;
4346             if (maxLayer<3){   // take prim vertex as additional measurement
4347               if (normdist[itrack1]>0 && htrackc1){
4348                 sumchi2 +=  TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4349               }else{
4350                 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4351               }
4352               sumn    +=  3-maxLayer;
4353             }
4354             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4355               sumn+=1.;
4356               if (btrack->GetClIndex(ilayer)<0){
4357                 sumchi2+=30;
4358                 continue;
4359               }else{
4360                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4361                 for (Int_t itrack=0;itrack<4;itrack++){
4362                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4363                     sumchi2+=18.;  //shared cluster
4364                     break;
4365                   }
4366                 }
4367                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4368                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4369               }
4370             }
4371             sumchi2/=sumn;
4372             if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4373           }
4374           continue;   //safety space - Geo manager will give exact layer
4375         }
4376         track1b = btrack;
4377         minchi2after1 = btrack->GetNormChi2(i);
4378         break;
4379       }
4380       //
4381       // position resolution - used for DCA cut
4382       Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4383         (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4384         (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4385       sigmad =TMath::Sqrt(sigmad)+0.04;
4386       if (pvertex->GetRr()>50){
4387         Double_t cov0[15],cov1[15];
4388         track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4389         track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4390         sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4391           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4392           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4393         sigmad =TMath::Sqrt(sigmad)+0.05;
4394       }
4395       //
4396       AliV0 vertex2;
4397       vertex2.SetParamN(*track0b);
4398       vertex2.SetParamP(*track1b);
4399       vertex2.Update(primvertex);
4400       //Bo:      if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4401       if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4402         pvertex->SetParamN(*track0b);
4403         pvertex->SetParamP(*track1b);
4404         pvertex->Update(primvertex);
4405         pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex());  // register clusters
4406         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4407         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4408       }
4409       pvertex->SetDistSigma(sigmad);
4410       //Bo:      pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4411       pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4412       //
4413       // define likelihhod and causalities
4414       //
4415       Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4416       if (maxLayer<1){
4417         Float_t fnorm0 = normdist[itrack0];
4418         if (fnorm0<0) fnorm0*=-3;
4419         Float_t fnorm1 = normdist[itrack1];
4420         if (fnorm1<0) fnorm1*=-3;
4421         if ((pvertex->GetAnglep()[2]>0.1) || ( (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 ) || (pvertex->GetRr()<3)){
4422           pb0    =  TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4423           pb1    =  TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4424         }
4425         pvertex->SetChi2Before(normdist[itrack0]);
4426         pvertex->SetChi2After(normdist[itrack1]);
4427         pvertex->SetNAfter(0);
4428         pvertex->SetNBefore(0);
4429       }else{
4430         pvertex->SetChi2Before(minchi2before0);
4431         pvertex->SetChi2After(minchi2before1);
4432          if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4433            pb0    =  TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4434            pb1    =  TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4435          }
4436          pvertex->SetNAfter(maxLayer);
4437          pvertex->SetNBefore(maxLayer);
4438       }
4439       if (pvertex->GetRr()<90){
4440         pa0  *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4441         pa1  *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4442       }
4443       if (pvertex->GetRr()<20){
4444         pa0  *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4445         pa1  *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4446       }
4447       //
4448       pvertex->SetCausality(pb0,pb1,pa0,pa1);
4449       //
4450       //  Likelihood calculations  - apply cuts
4451       //
4452       Bool_t v0OK = kTRUE;
4453       Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4454       p12        += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4455       p12         = TMath::Sqrt(p12);                             // "mean" momenta
4456       Float_t    sigmap0   = 0.0001+0.001/(0.1+pvertex->GetRr());
4457       Float_t    sigmap    = 0.5*sigmap0*(0.6+0.4*p12);           // "resolution: of point angle - as a function of radius and momenta
4458
4459       Float_t causalityA  = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4460       Float_t causalityB  = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4461                                         TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4462       //
4463       //Bo:      Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4464       Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4465       Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4466
4467       Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4468         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4469         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4470         0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4471       //
4472       if (causalityA<kCausality0Cut)                                          v0OK = kFALSE;
4473       if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut)              v0OK = kFALSE;
4474       if (likelihood1<kLikelihood1Cut)                                        v0OK = kFALSE;
4475       if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4476
4477       //
4478       //
4479       if (kFALSE){
4480         Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4481         cstream<<"It0"<<
4482           "Tr0.="<<track0<<                       //best without constrain
4483           "Tr1.="<<track1<<                       //best without constrain
4484           "Tr0B.="<<track0b<<                     //best without constrain  after vertex
4485           "Tr1B.="<<track1b<<                     //best without constrain  after vertex
4486           "Tr0C.="<<htrackc0<<                    //best with constrain     if exist
4487           "Tr1C.="<<htrackc1<<                    //best with constrain     if exist
4488           "Tr0L.="<<track0l<<                     //longest best
4489           "Tr1L.="<<track1l<<                     //longest best
4490           "Esd0.="<<track0->GetESDtrack()<<           // esd track0 params
4491           "Esd1.="<<track1->GetESDtrack()<<           // esd track1 params
4492           "V0.="<<pvertex<<                       //vertex properties
4493           "V0b.="<<&vertex2<<                       //vertex properties at "best" track
4494           "ND0="<<normdist[itrack0]<<             //normalize distance for track0
4495           "ND1="<<normdist[itrack1]<<             //normalize distance for track1
4496           "Gold="<<gold<<                         //
4497           //      "RejectBase="<<rejectBase<<             //rejection in First itteration
4498           "OK="<<v0OK<<
4499           "rmin="<<rmin<<
4500           "sigmad="<<sigmad<<
4501           "\n";
4502       }
4503       //if (rejectBase) continue;
4504       //
4505       pvertex->SetStatus(0);
4506       //      if (rejectBase) {
4507       //        pvertex->SetStatus(-100);
4508       //}
4509       if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4510         //Bo:     pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4511         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4512         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4513         if (v0OK){
4514           //      AliV0vertex vertexjuri(*track0,*track1);
4515           //      vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4516           //      event->AddV0(&vertexjuri);
4517           pvertex->SetStatus(100);
4518         }
4519         pvertex->SetOnFlyStatus(kTRUE);
4520         pvertex->ChangeMassHypothesis(kK0Short);
4521         event->AddV0(pvertex);
4522       }
4523     }
4524   }
4525   //
4526   //
4527   // delete temporary arrays
4528   //
4529   delete[] forbidden;
4530   delete[] minPointAngle;
4531   delete[] maxr;
4532   delete[] minr;
4533   delete[] norm;
4534   delete[] normdist;
4535   delete[] normdist1;
4536   delete[] normdist0;
4537   delete[] dist;
4538   delete[] itsmap;
4539   delete[] helixes;
4540   delete   pvertex;
4541 }
4542 //------------------------------------------------------------------------
4543 void AliITStrackerMI::RefitV02(const AliESDEvent *event)
4544 {
4545   //
4546   //try to refit  V0s in the third path of the reconstruction
4547   //
4548   TTreeSRedirector &cstream = *fDebugStreamer;
4549   //
4550   Int_t  nv0s = event->GetNumberOfV0s();
4551   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4552   AliV0 v0temp;
4553   for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4554     AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4555     if (!v0mi) continue;
4556     Int_t     itrack0   = v0mi->GetIndex(0);
4557     Int_t     itrack1   = v0mi->GetIndex(1);
4558     AliESDtrack *esd0   = event->GetTrack(itrack0);
4559     AliESDtrack *esd1   = event->GetTrack(itrack1);
4560     if (!esd0||!esd1) continue;
4561     AliITStrackMI tpc0(*esd0);
4562     AliITStrackMI tpc1(*esd1);
4563     Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4564     Double_t alpha =TMath::ATan2(y,x);   //I.B.
4565     if (v0mi->GetRr()>85){
4566       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4567         v0temp.SetParamN(tpc0);
4568         v0temp.SetParamP(tpc1);
4569         v0temp.Update(primvertex);
4570         if (kFALSE) cstream<<"Refit"<<
4571           "V0.="<<v0mi<<
4572           "V0refit.="<<&v0temp<<
4573           "Tr0.="<<&tpc0<<
4574           "Tr1.="<<&tpc1<<
4575           "\n";
4576         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4577         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4578           v0mi->SetParamN(tpc0);
4579           v0mi->SetParamP(tpc1);
4580           v0mi->Update(primvertex);
4581         }
4582       }
4583       continue;
4584     }
4585     if (v0mi->GetRr()>35){
4586       CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4587       CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4588       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4589         v0temp.SetParamN(tpc0);
4590         v0temp.SetParamP(tpc1);
4591         v0temp.Update(primvertex);
4592         if (kFALSE) cstream<<"Refit"<<
4593           "V0.="<<v0mi<<
4594           "V0refit.="<<&v0temp<<
4595           "Tr0.="<<&tpc0<<
4596           "Tr1.="<<&tpc1<<
4597           "\n";
4598         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4599         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4600           v0mi->SetParamN(tpc0);
4601           v0mi->SetParamP(tpc1);
4602           v0mi->Update(primvertex);
4603         }
4604       }
4605       continue;
4606     }
4607     CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4608     CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4609     //    if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4610     if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4611       v0temp.SetParamN(tpc0);
4612       v0temp.SetParamP(tpc1);
4613       v0temp.Update(primvertex);
4614       if (kFALSE) cstream<<"Refit"<<
4615         "V0.="<<v0mi<<
4616         "V0refit.="<<&v0temp<<
4617         "Tr0.="<<&tpc0<<
4618         "Tr1.="<<&tpc1<<
4619         "\n";
4620       //Bo:      if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4621       if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4622         v0mi->SetParamN(tpc0);
4623         v0mi->SetParamP(tpc1);
4624         v0mi->Update(primvertex);
4625       }
4626     }
4627   }
4628 }
4629
4630 //------------------------------------------------------------------------
4631 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4632   //--------------------------------------------------------------------
4633   // Fill a look-up table with mean material
4634   //--------------------------------------------------------------------
4635
4636   Int_t n=1000;
4637   Double_t mparam[7];
4638   Double_t point1[3],point2[3];
4639   Double_t phi,cosphi,sinphi,z;
4640   // 0-5 layers, 6 pipe, 7-8 shields
4641   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4642   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4643
4644   Int_t ifirst=0,ilast=0;
4645   if(material.Contains("Pipe")) {
4646     ifirst=6; ilast=6;
4647   } else if(material.Contains("Shields")) {
4648     ifirst=7; ilast=8;
4649   } else if(material.Contains("Layers")) {
4650     ifirst=0; ilast=5;
4651   } else {
4652     Error("BuildMaterialLUT","Wrong layer name\n");
4653   }
4654
4655   for(Int_t imat=ifirst; imat<=ilast; imat++) {
4656     Double_t param[5]={0.,0.,0.,0.,0.};
4657     for (Int_t i=0; i<n; i++) {
4658       phi = 2.*TMath::Pi()*gRandom->Rndm();
4659       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4660       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4661       point1[0] = rmin[imat]*cosphi;
4662       point1[1] = rmin[imat]*sinphi;
4663       point1[2] = z;
4664       point2[0] = rmax[imat]*cosphi;
4665       point2[1] = rmax[imat]*sinphi;
4666       point2[2] = z;
4667       AliTracker::MeanMaterialBudget(point1,point2,mparam);
4668       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4669     }
4670     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4671     if(imat<=5) {
4672       fxOverX0Layer[imat] = param[1];
4673       fxTimesRhoLayer[imat] = param[0]*param[4];
4674     } else if(imat==6) {
4675       fxOverX0Pipe = param[1];
4676       fxTimesRhoPipe = param[0]*param[4];
4677     } else if(imat==7) {
4678       fxOverX0Shield[0] = param[1];
4679       fxTimesRhoShield[0] = param[0]*param[4];
4680     } else if(imat==8) {
4681       fxOverX0Shield[1] = param[1];
4682       fxTimesRhoShield[1] = param[0]*param[4];
4683     }
4684   }
4685   /*
4686   printf("%s\n",material.Data());
4687   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4688   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4689   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4690   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4691   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4692   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4693   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4694   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4695   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4696   */
4697   return;
4698 }
4699 //------------------------------------------------------------------------
4700 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4701                                               TString direction) {
4702   //-------------------------------------------------------------------
4703   // Propagate beyond beam pipe and correct for material
4704   // (material budget in different ways according to fUseTGeo value)
4705   // Add time if going outward (PropagateTo or PropagateToTGeo)
4706   //-------------------------------------------------------------------
4707
4708   // Define budget mode:
4709   // 0: material from AliITSRecoParam (hard coded)
4710   // 1: material from TGeo in one step (on the fly)
4711   // 2: material from lut
4712   // 3: material from TGeo in one step (same for all hypotheses)
4713   Int_t mode;
4714   switch(fUseTGeo) {
4715   case 0:
4716     mode=0;
4717     break;
4718   case 1:
4719     mode=1;
4720     break;
4721   case 2:
4722     mode=2;
4723     break;
4724   case 3:
4725     if(fTrackingPhase.Contains("Clusters2Tracks"))
4726       { mode=3; } else { mode=1; }
4727     break;
4728   case 4:
4729     if(fTrackingPhase.Contains("Clusters2Tracks"))
4730       { mode=3; } else { mode=2; }
4731     break;
4732   default:
4733     mode=0;
4734     break;
4735   }
4736   if(fTrackingPhase.Contains("Default")) mode=0;
4737
4738   Int_t index=fCurrentEsdTrack;
4739
4740   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4741   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4742   Double_t xToGo;
4743   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4744
4745   Double_t xOverX0,x0,lengthTimesMeanDensity;
4746
4747   switch(mode) {
4748   case 0:
4749     xOverX0 = AliITSRecoParam::GetdPipe();
4750     x0 = AliITSRecoParam::GetX0Be();
4751     lengthTimesMeanDensity = xOverX0*x0;
4752     lengthTimesMeanDensity *= dir;
4753     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4754     break;
4755   case 1:
4756     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4757     break;
4758   case 2:
4759     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4760     xOverX0 = fxOverX0Pipe;
4761     lengthTimesMeanDensity = fxTimesRhoPipe;
4762     lengthTimesMeanDensity *= dir;
4763     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4764     break;
4765   case 3:
4766     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4767     if(fxOverX0PipeTrks[index]<0) {
4768       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4769       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4770                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4771       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4772       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4773       return 1;
4774     }
4775     xOverX0 = fxOverX0PipeTrks[index];
4776     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4777     lengthTimesMeanDensity *= dir;
4778     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4779     break;
4780   }
4781
4782   return 1;
4783 }
4784 //------------------------------------------------------------------------
4785 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4786                                                 TString shield,
4787                                                 TString direction) {
4788   //-------------------------------------------------------------------
4789   // Propagate beyond SPD or SDD shield and correct for material
4790   // (material budget in different ways according to fUseTGeo value)
4791   // Add time if going outward (PropagateTo or PropagateToTGeo)
4792   //-------------------------------------------------------------------
4793
4794   // Define budget mode:
4795   // 0: material from AliITSRecoParam (hard coded)
4796   // 1: material from TGeo in steps of X cm (on the fly)
4797   //    X = AliITSRecoParam::GetStepSizeTGeo()
4798   // 2: material from lut
4799   // 3: material from TGeo in one step (same for all hypotheses)
4800   Int_t mode;
4801   switch(fUseTGeo) {
4802   case 0:
4803     mode=0;
4804     break;
4805   case 1:
4806     mode=1;
4807     break;
4808   case 2:
4809     mode=2;
4810     break;
4811   case 3:
4812     if(fTrackingPhase.Contains("Clusters2Tracks"))
4813       { mode=3; } else { mode=1; }
4814     break;
4815   case 4:
4816     if(fTrackingPhase.Contains("Clusters2Tracks"))
4817       { mode=3; } else { mode=2; }
4818     break;
4819   default:
4820     mode=0;
4821     break;
4822   }
4823   if(fTrackingPhase.Contains("Default")) mode=0;
4824
4825   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4826   Double_t rToGo;
4827   Int_t    shieldindex=0;
4828   if (shield.Contains("SDD")) { // SDDouter
4829     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4830     shieldindex=1;
4831   } else if (shield.Contains("SPD")) {        // SPDouter
4832     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4833     shieldindex=0;
4834   } else {
4835     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4836     return 0;
4837   }
4838   Double_t xToGo;
4839   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4840
4841   Int_t index=2*fCurrentEsdTrack+shieldindex;
4842
4843   Double_t xOverX0,x0,lengthTimesMeanDensity;
4844   Int_t nsteps=1;
4845
4846   switch(mode) {
4847   case 0:
4848     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4849     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4850     lengthTimesMeanDensity = xOverX0*x0;
4851     lengthTimesMeanDensity *= dir;
4852     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4853     break;
4854   case 1:
4855     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4856     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4857     break;
4858   case 2:
4859     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4860     xOverX0 = fxOverX0Shield[shieldindex];
4861     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4862     lengthTimesMeanDensity *= dir;
4863     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4864     break;
4865   case 3:
4866     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4867     if(fxOverX0ShieldTrks[index]<0) {
4868       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4869       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4870                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4871       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4872       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4873       return 1;
4874     }
4875     xOverX0 = fxOverX0ShieldTrks[index];
4876     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4877     lengthTimesMeanDensity *= dir;
4878     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4879     break;
4880   }
4881
4882   return 1;
4883 }
4884 //------------------------------------------------------------------------
4885 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4886                                                Int_t layerindex,
4887                                                Double_t oldGlobXYZ[3],
4888                                                TString direction) {
4889   //-------------------------------------------------------------------
4890   // Propagate beyond layer and correct for material
4891   // (material budget in different ways according to fUseTGeo value)
4892   // Add time if going outward (PropagateTo or PropagateToTGeo)
4893   //-------------------------------------------------------------------
4894
4895   // Define budget mode:
4896   // 0: material from AliITSRecoParam (hard coded)
4897   // 1: material from TGeo in stepsof X cm (on the fly)
4898   //    X = AliITSRecoParam::GetStepSizeTGeo()
4899   // 2: material from lut
4900   // 3: material from TGeo in one step (same for all hypotheses)
4901   Int_t mode;
4902   switch(fUseTGeo) {
4903   case 0:
4904     mode=0;
4905     break;
4906   case 1:
4907     mode=1;
4908     break;
4909   case 2:
4910     mode=2;
4911     break;
4912   case 3:
4913     if(fTrackingPhase.Contains("Clusters2Tracks"))
4914       { mode=3; } else { mode=1; }
4915     break;
4916   case 4:
4917     if(fTrackingPhase.Contains("Clusters2Tracks"))
4918       { mode=3; } else { mode=2; }
4919     break;
4920   default:
4921     mode=0;
4922     break;
4923   }
4924   if(fTrackingPhase.Contains("Default")) mode=0;
4925
4926   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4927
4928   Double_t r=fgLayers[layerindex].GetR();
4929   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4930
4931   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4932   Double_t xToGo;
4933   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4934
4935   Int_t index=6*fCurrentEsdTrack+layerindex;
4936
4937
4938   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4939   Int_t nsteps=1;
4940
4941   // back before material (no correction)
4942   Double_t rOld,xOld;
4943   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4944   if (!t->GetLocalXat(rOld,xOld)) return 0;
4945   if (!t->Propagate(xOld)) return 0;
4946
4947   switch(mode) {
4948   case 0:
4949     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4950     lengthTimesMeanDensity = xOverX0*x0;
4951     lengthTimesMeanDensity *= dir;
4952     // Bring the track beyond the material
4953     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4954     break;
4955   case 1:
4956     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4957     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4958     break;
4959   case 2:
4960     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4961     xOverX0 = fxOverX0Layer[layerindex];
4962     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4963     lengthTimesMeanDensity *= dir;
4964     // Bring the track beyond the material
4965     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4966     break;
4967   case 3:
4968     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4969     if(fxOverX0LayerTrks[index]<0) {
4970       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4971       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4972       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4973                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4974       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4975       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4976       return 1;
4977     }
4978     xOverX0 = fxOverX0LayerTrks[index];
4979     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4980     lengthTimesMeanDensity *= dir;
4981     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4982     break;
4983   }
4984
4985
4986   return 1;
4987 }
4988 //------------------------------------------------------------------------
4989 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4990   //-----------------------------------------------------------------
4991   // Initialize LUT for storing material for each prolonged track
4992   //-----------------------------------------------------------------
4993   fxOverX0PipeTrks = new Float_t[ntracks];
4994   fxTimesRhoPipeTrks = new Float_t[ntracks];
4995   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4996   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4997   fxOverX0LayerTrks = new Float_t[ntracks*6];
4998   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4999
5000   for(Int_t i=0; i<ntracks; i++) {
5001     fxOverX0PipeTrks[i] = -1.;
5002     fxTimesRhoPipeTrks[i] = -1.;
5003   }
5004   for(Int_t j=0; j<ntracks*2; j++) {
5005     fxOverX0ShieldTrks[j] = -1.;
5006     fxTimesRhoShieldTrks[j] = -1.;
5007   }
5008   for(Int_t k=0; k<ntracks*6; k++) {
5009     fxOverX0LayerTrks[k] = -1.;
5010     fxTimesRhoLayerTrks[k] = -1.;
5011   }
5012
5013   fNtracks = ntracks;
5014
5015   return;
5016 }
5017 //------------------------------------------------------------------------
5018 void AliITStrackerMI::DeleteTrksMaterialLUT() {
5019   //-----------------------------------------------------------------
5020   // Delete LUT for storing material for each prolonged track
5021   //-----------------------------------------------------------------
5022   if(fxOverX0PipeTrks) {
5023     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
5024   }
5025   if(fxOverX0ShieldTrks) {
5026     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
5027   }
5028
5029   if(fxOverX0LayerTrks) {
5030     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
5031   }
5032   if(fxTimesRhoPipeTrks) {
5033     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
5034   }
5035   if(fxTimesRhoShieldTrks) {
5036     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
5037   }
5038   if(fxTimesRhoLayerTrks) {
5039     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
5040   }
5041   return;
5042 }
5043 //------------------------------------------------------------------------
5044 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
5045                                       Int_t ilayer,Int_t idet) const {
5046   //-----------------------------------------------------------------
5047   // This method is used to decide whether to allow a prolongation
5048   // without clusters, because we want to skip the layer.
5049   // In this case the return value is > 0:
5050   // return 1: the user requested to skip a layer
5051   // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
5052   //-----------------------------------------------------------------
5053
5054   if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
5055
5056   if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
5057     // check if track will cross SPD outer layer
5058     Double_t phiAtSPD2,zAtSPD2;
5059     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
5060       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
5061     }
5062   }
5063
5064   return 0;
5065 }
5066 //------------------------------------------------------------------------
5067 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
5068                                      Int_t ilayer,Int_t idet,
5069                                      Double_t dz,Double_t dy,
5070                                      Bool_t noClusters) const {
5071   //-----------------------------------------------------------------
5072   // This method is used to decide whether to allow a prolongation
5073   // without clusters, because there is a dead zone in the road.
5074   // In this case the return value is > 0:
5075   // return 1: dead zone at z=0,+-7cm in SPD
5076   // return 2: all road is "bad" (dead or noisy) from the OCDB
5077   // return 3: something "bad" (dead or noisy) from the OCDB
5078   //-----------------------------------------------------------------
5079
5080   // check dead zones at z=0,+-7cm in the SPD
5081   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
5082     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5083                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5084                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
5085     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5086                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5087                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
5088     for (Int_t i=0; i<3; i++)
5089       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
5090         AliDebug(2,Form("crack SPD %d",ilayer));
5091         return 1;
5092       }
5093   }
5094
5095   // check bad zones from OCDB
5096   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
5097
5098   if (idet<0) return 0;
5099
5100   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
5101
5102   Int_t detType=-1;
5103   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
5104   if (ilayer==0 || ilayer==1) {        // ----------  SPD
5105     detType = 0;
5106   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
5107     detType = 1;
5108     detSizeFactorX *= 2.;
5109   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
5110     detType = 2;
5111   }
5112   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
5113   if (detType==2) segm->SetLayer(ilayer+1);
5114   Float_t detSizeX = detSizeFactorX*segm->Dx();
5115   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
5116
5117   // check if the road overlaps with bad chips
5118   Float_t xloc,zloc;
5119   LocalModuleCoord(ilayer,idet,track,xloc,zloc);
5120   Float_t zlocmin = zloc-dz;
5121   Float_t zlocmax = zloc+dz;
5122   Float_t xlocmin = xloc-dy;
5123   Float_t xlocmax = xloc+dy;
5124   Int_t chipsInRoad[100];
5125
5126   // check if road goes out of detector
5127   Bool_t touchNeighbourDet=kFALSE;
5128   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
5129   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
5130   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5131   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5132   AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f   %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
5133
5134   // check if this detector is bad
5135   if (det.IsBad()) {
5136     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
5137     if(!touchNeighbourDet) {
5138       return 2; // all detectors in road are bad
5139     } else {
5140       return 3; // at least one is bad
5141     }
5142   }
5143
5144   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
5145   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
5146   if (!nChipsInRoad) return 0;
5147
5148   Bool_t anyBad=kFALSE,anyGood=kFALSE;
5149   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
5150     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
5151     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
5152     if (det.IsChipBad(chipsInRoad[iCh])) {
5153       anyBad=kTRUE;
5154     } else {
5155       anyGood=kTRUE;
5156     }
5157   }
5158
5159   if (!anyGood) {
5160     if(!touchNeighbourDet) {
5161       AliDebug(2,"all bad in road");
5162       return 2;  // all chips in road are bad
5163     } else {
5164       return 3; // at least a bad chip in road
5165     }
5166   }
5167
5168   if (anyBad) {
5169     AliDebug(2,"at least a bad in road");
5170     return 3; // at least a bad chip in road
5171   }
5172
5173
5174   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
5175       || ilayer==4 || ilayer==5     // SSD
5176       || !noClusters) return 0;
5177
5178   // There are no clusters in road: check if there is at least
5179   // a bad SPD pixel or SDD anode
5180
5181   Int_t idetInITS=idet;
5182   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
5183
5184   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
5185     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
5186     return 3;
5187   }
5188   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
5189
5190   return 0;
5191 }
5192 //------------------------------------------------------------------------
5193 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
5194                                        const AliITStrackMI *track,
5195                                        Float_t &xloc,Float_t &zloc) const {
5196   //-----------------------------------------------------------------
5197   // Gives position of track in local module ref. frame
5198   //-----------------------------------------------------------------
5199
5200   xloc=0.;
5201   zloc=0.;
5202
5203   if(idet<0) return kFALSE;
5204
5205   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
5206
5207   Int_t lad = Int_t(idet/ndet) + 1;
5208
5209   Int_t det = idet - (lad-1)*ndet + 1;
5210
5211   Double_t xyzGlob[3],xyzLoc[3];
5212
5213   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
5214   // take into account the misalignment: xyz at real detector plane
5215   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
5216
5217   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
5218
5219   xloc = (Float_t)xyzLoc[0];
5220   zloc = (Float_t)xyzLoc[2];
5221
5222   return kTRUE;
5223 }
5224 //------------------------------------------------------------------------
5225 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
5226 //
5227 // Method to be optimized further:
5228 // Aim: decide whether a track can be used for PlaneEff evaluation
5229 //      the decision is taken based on the track quality at the layer under study
5230 //      no information on the clusters on this layer has to be used
5231 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
5232 //      the cut is done on number of sigmas from the boundaries
5233 //
5234 //  Input: Actual track, layer [0,5] under study
5235 //  Output: none
5236 //  Return: kTRUE if this is a good track
5237 //
5238 // it will apply a pre-selection to obtain good quality tracks.
5239 // Here also  you will have the possibility to put a control on the
5240 // impact point of the track on the basic block, in order to exclude border regions
5241 // this will be done by calling a proper method of the AliITSPlaneEff class.
5242 //
5243 // input: AliITStrackMI* track, ilayer= layer number [0,5]
5244 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
5245 //
5246   Int_t index[AliITSgeomTGeo::kNLayers];
5247   Int_t k;
5248   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
5249   //
5250   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
5251     index[k]=clusters[k];
5252   }
5253
5254   if(!fPlaneEff)
5255     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
5256   AliITSlayer &layer=fgLayers[ilayer];
5257   Double_t r=layer.GetR();
5258   AliITStrackMI tmp(*track);
5259
5260 // require a minimal number of cluster in other layers and eventually clusters in closest layers
5261   Int_t ncl=0;
5262   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
5263     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
5264                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
5265     if (tmp.GetClIndex(lay)>=0) ncl++;
5266   }
5267   Bool_t nextout = kFALSE;
5268   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
5269   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
5270   Bool_t nextin = kFALSE;
5271   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
5272   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
5273   if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
5274      return kFALSE;
5275   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
5276   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
5277   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
5278  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
5279
5280 // detector number
5281   Double_t phi,z;
5282   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
5283   Int_t idet=layer.FindDetectorIndex(phi,z);
5284   if(idet<0) { AliInfo(Form("cannot find detector"));
5285     return kFALSE;}
5286
5287   // here check if it has good Chi Square.
5288
5289   //propagate to the intersection with the detector plane
5290   const AliITSdetector &det=layer.GetDetector(idet);
5291   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
5292
5293   Float_t locx; //
5294   Float_t locz; //
5295   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
5296   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5297   if(key>fPlaneEff->Nblock()) return kFALSE;
5298   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
5299   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
5300   //***************
5301   // DEFINITION OF SEARCH ROAD FOR accepting a track
5302   //
5303   //For the time being they are hard-wired, later on from AliITSRecoParam
5304   // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
5305   // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
5306   Double_t nsigz=4;
5307   Double_t nsigx=4;
5308   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
5309   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
5310                                                 // done for RecPoints
5311
5312   // exclude tracks at boundary between detectors
5313   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
5314   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
5315   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
5316   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
5317   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5318
5319   if ( (locx-dx < blockXmn+boundaryWidth) ||
5320        (locx+dx > blockXmx-boundaryWidth) ||
5321        (locz-dz < blockZmn+boundaryWidth) ||
5322        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5323   return kTRUE;
5324 }
5325 //------------------------------------------------------------------------
5326 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
5327 //
5328 // This Method has to be optimized! For the time-being it uses the same criteria
5329 // as those used in the search of extra clusters for overlapping modules.
5330 //
5331 // Method Purpose: estabilish whether a track has produced a recpoint or not
5332 //                 in the layer under study (For Plane efficiency)
5333 //
5334 // inputs: AliITStrackMI* track  (pointer to a usable track)
5335 // outputs: none
5336 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5337 //               traversed by this very track. In details:
5338 //               - if a cluster can be associated to the track then call
5339 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5340 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5341 //
5342   if(!fPlaneEff)
5343     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5344   AliITSlayer &layer=fgLayers[ilayer];
5345   Double_t r=layer.GetR();
5346   AliITStrackMI tmp(*track);
5347
5348 // detector number
5349   Double_t phi,z;
5350   if (!tmp.GetPhiZat(r,phi,z)) return;
5351   Int_t idet=layer.FindDetectorIndex(phi,z);
5352
5353   if(idet<0) { AliInfo(Form("cannot find detector"));
5354     return;}
5355
5356
5357 //propagate to the intersection with the detector plane
5358   const AliITSdetector &det=layer.GetDetector(idet);
5359   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5360
5361
5362 //***************
5363 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5364 //
5365   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5366                     TMath::Sqrt(tmp.GetSigmaZ2() +
5367                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5368                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5369                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5370   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5371                     TMath::Sqrt(tmp.GetSigmaY2() +
5372                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5373                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5374                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5375
5376 // road in global (rphi,z) [i.e. in tracking ref. system]
5377   Double_t zmin = tmp.GetZ() - dz;
5378   Double_t zmax = tmp.GetZ() + dz;
5379   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5380   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5381
5382 // select clusters in road
5383   layer.SelectClusters(zmin,zmax,ymin,ymax);
5384
5385 // Define criteria for track-cluster association
5386   Double_t msz = tmp.GetSigmaZ2() +
5387   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5388   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5389   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5390   Double_t msy = tmp.GetSigmaY2() +
5391   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5392   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5393   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5394   if (tmp.GetConstrain()) {
5395     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5396     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5397   }  else {
5398     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5399     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5400   }
5401   msz = 1./msz; // 1/RoadZ^2
5402   msy = 1./msy; // 1/RoadY^2
5403 //
5404
5405   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5406   Int_t idetc=-1;
5407   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5408   //Double_t  tolerance=0.2;
5409   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5410     idetc = cl->GetDetectorIndex();
5411     if(idet!=idetc) continue;
5412     //Int_t ilay = cl->GetLayer();
5413
5414     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5415     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5416
5417     Double_t chi2=tmp.GetPredictedChi2(cl);
5418     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5419   }*/
5420   Float_t locx; //
5421   Float_t locz; //
5422   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5423 //
5424   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5425   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5426   if(key>fPlaneEff->Nblock()) return;
5427   Bool_t found=kFALSE;
5428   //if (ci>=0) {
5429   Double_t chi2;
5430   while ((cl=layer.GetNextCluster(clidx))!=0) {
5431     idetc = cl->GetDetectorIndex();
5432     if(idet!=idetc) continue;
5433     // here real control to see whether the cluster can be associated to the track.
5434     // cluster not associated to track
5435     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5436          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
5437     // calculate track-clusters chi2
5438     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5439                                                // in particular, the error associated to the cluster
5440     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5441     // chi2 cut
5442     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5443     found=kTRUE;
5444     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5445    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5446    // track->SetExtraModule(ilayer,idetExtra);
5447   }
5448   if(!fPlaneEff->UpDatePlaneEff(found,key))
5449        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5450   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5451     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
5452     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5453     Int_t cltype[2]={-999,-999};
5454
5455     tr[0]=locx;
5456     tr[1]=locz;
5457     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
5458     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
5459
5460     if (found){
5461       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5462       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5463       cltype[0]=layer.GetCluster(ci)->GetNy();
5464       cltype[1]=layer.GetCluster(ci)->GetNz();
5465
5466      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5467      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
5468      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5469      // It is computed properly by calling the method
5470      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5471      // T
5472      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5473       //if (tmp.PropagateTo(x,0.,0.)) {
5474         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5475         AliCluster c(*layer.GetCluster(ci));
5476         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5477         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5478         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
5479         clu[2]=TMath::Sqrt(c.GetSigmaY2());
5480         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5481       //}
5482     }
5483     fPlaneEff->FillHistos(key,found,tr,clu,cltype);
5484   }
5485 return;
5486 }