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 <TTreeStream.h>
  32 #include <TDatabasePDG.h>
  33 #include <TString.h>
  34 #include <TRandom.h>
  35
  36
  37 #include "AliLog.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDtrack.h"
  40 #include "AliESDVertex.h"
  41 #include "AliV0.h"
  42 #include "AliHelix.h"
  43 #include "AliITSRecPoint.h"
  44 #include "AliITSgeomTGeo.h"
  45 #include "AliITSReconstructor.h"
  46 #include "AliTrackPointArray.h"
  47 #include "AliAlignObj.h"
  48 #include "AliITSClusterParam.h"
  49 #include "AliCDBManager.h"
  50 #include "AliCDBEntry.h"
  51 #include "AliITSsegmentation.h"
  52 #include "AliITSCalibration.h"
  53 #include "AliITSCalibrationSPD.h"
  54 #include "AliITSCalibrationSDD.h"
  55 #include "AliITSCalibrationSSD.h"
  56 #include "AliITSPlaneEff.h"
  57 #include "AliITSPlaneEffSPD.h"
  58 #include "AliITSPlaneEffSDD.h"
  59 #include "AliITSPlaneEffSSD.h"
  60 #include "AliITStrackerMI.h"
  61
  62 ClassImp(AliITStrackerMI)
  63
  64 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
  65
  66 AliITStrackerMI::AliITStrackerMI():AliTracker(),
  67 fI(0),
  68 fBestTrack(),
  69 fTrackToFollow(),
  70 fTrackHypothesys(),
  71 fBestHypothesys(),
  72 fOriginal(),
  73 fCurrentEsdTrack(),
  74 fPass(0),
  75 fAfterV0(kFALSE),
  76 fLastLayerToTrackTo(0),
  77 fCoefficients(0),
  78 fEsd(0),
  79 fTrackingPhase("Default"),
  80 fUseTGeo(3),
  81 fNtracks(0),
  82 fxOverX0Pipe(-1.),
  83 fxTimesRhoPipe(-1.),
  84 fxOverX0PipeTrks(0),
  85 fxTimesRhoPipeTrks(0),
  86 fxOverX0ShieldTrks(0),
  87 fxTimesRhoShieldTrks(0),
  88 fxOverX0LayerTrks(0),
  89 fxTimesRhoLayerTrks(0),
  90 fDebugStreamer(0),
  91 fITSChannelStatus(0),
  92 fDetTypeRec(0),
  93 fPlaneEff(0) {
  94   //Default constructor
  95   Int_t i;
  96   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
  97   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
  98   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
  99 }
 100 //------------------------------------------------------------------------
 101 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
 102 fI(AliITSgeomTGeo::GetNLayers()),
 103 fBestTrack(),
 104 fTrackToFollow(),
 105 fTrackHypothesys(),
 106 fBestHypothesys(),
 107 fOriginal(),
 108 fCurrentEsdTrack(),
 109 fPass(0),
 110 fAfterV0(kFALSE),
 111 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
 112 fCoefficients(0),
 113 fEsd(0),
 114 fTrackingPhase("Default"),
 115 fUseTGeo(3),
 116 fNtracks(0),
 117 fxOverX0Pipe(-1.),
 118 fxTimesRhoPipe(-1.),
 119 fxOverX0PipeTrks(0),
 120 fxTimesRhoPipeTrks(0),
 121 fxOverX0ShieldTrks(0),
 122 fxTimesRhoShieldTrks(0),
 123 fxOverX0LayerTrks(0),
 124 fxTimesRhoLayerTrks(0),
 125 fDebugStreamer(0),
 126 fITSChannelStatus(0),
 127 fDetTypeRec(0),
 128 fPlaneEff(0) {
 129   //--------------------------------------------------------------------
 130   //This is the AliITStrackerMI constructor
 131   //--------------------------------------------------------------------
 132   if (geom) {
 133     AliWarning("\"geom\" is actually a dummy argument !");
 134   }
 135
 136   fCoefficients = 0;
 137   fAfterV0     = kFALSE;
 138
 139   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 140     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 141     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 142
 143     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 144     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 145     Double_t poff=TMath::ATan2(y,x);
 146     Double_t zoff=z;
 147     Double_t r=TMath::Sqrt(x*x + y*y);
 148
 149     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 150     r += TMath::Sqrt(x*x + y*y);
 151     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 152     r += TMath::Sqrt(x*x + y*y);
 153     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 154     r += TMath::Sqrt(x*x + y*y);
 155     r*=0.25;
 156
 157     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
 158
 159     for (Int_t j=1; j<nlad+1; j++) {
 160       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 161         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 162         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 163         m.Multiply(tm);
 164         Double_t txyz[3]={0.};
 165         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 166         m.LocalToMaster(txyz,xyz);
 167         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 168         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 169
 170         if (phi<0) phi+=TMath::TwoPi();
 171         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 172
 173         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 174         new(&det) AliITSdetector(r,phi);
 175         // compute the real radius (with misalignment)
 176         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 177         mmisal.Multiply(tm);
 178         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 179         mmisal.LocalToMaster(txyz,xyz);
 180         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 181         det.SetRmisal(rmisal);
 182
 183       } // end loop on detectors
 184     } // end loop on ladders
 185   } // end loop on layers
 186
 187
 188   fI=AliITSgeomTGeo::GetNLayers();
 189
 190   fPass=0;
 191   fConstraint[0]=1; fConstraint[1]=0;
 192
 193   Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
 194                      AliITSReconstructor::GetRecoParam()->GetYVdef(),
 195                      AliITSReconstructor::GetRecoParam()->GetZVdef()};
 196   Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
 197                      AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
 198                      AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
 199   SetVertex(xyzVtx,ersVtx);
 200
 201   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=AliITSRecoParam::GetLayersNotToSkip(i);
 202   fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
 203   for (Int_t i=0;i<100000;i++){
 204     fBestTrackIndex[i]=0;
 205   }
 206
 207   // store positions of centre of SPD modules (in z)
 208   Double_t tr[3];
 209   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 210   fSPDdetzcentre[0] = tr[2];
 211   AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 212   fSPDdetzcentre[1] = tr[2];
 213   AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 214   fSPDdetzcentre[2] = tr[2];
 215   AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 216   fSPDdetzcentre[3] = tr[2];
 217
 218   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
 219   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 220     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 221     fUseTGeo = 3;
 222   }
 223
 224   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 225   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 226
 227   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 228
 229   // only for plane efficiency evaluation
 230   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff()) {
 231     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
 232     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane))
 233       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
 234     if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
 235     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
 236     else fPlaneEff = new AliITSPlaneEffSSD();
 237     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
 238        if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
 239     if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
 240   }
 241 }
 242 //------------------------------------------------------------------------
 243 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
 244 fI(tracker.fI),
 245 fBestTrack(tracker.fBestTrack),
 246 fTrackToFollow(tracker.fTrackToFollow),
 247 fTrackHypothesys(tracker.fTrackHypothesys),
 248 fBestHypothesys(tracker.fBestHypothesys),
 249 fOriginal(tracker.fOriginal),
 250 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
 251 fPass(tracker.fPass),
 252 fAfterV0(tracker.fAfterV0),
 253 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
 254 fCoefficients(tracker.fCoefficients),
 255 fEsd(tracker.fEsd),
 256 fTrackingPhase(tracker.fTrackingPhase),
 257 fUseTGeo(tracker.fUseTGeo),
 258 fNtracks(tracker.fNtracks),
 259 fxOverX0Pipe(tracker.fxOverX0Pipe),
 260 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 261 fxOverX0PipeTrks(0),
 262 fxTimesRhoPipeTrks(0),
 263 fxOverX0ShieldTrks(0),
 264 fxTimesRhoShieldTrks(0),
 265 fxOverX0LayerTrks(0),
 266 fxTimesRhoLayerTrks(0),
 267 fDebugStreamer(tracker.fDebugStreamer),
 268 fITSChannelStatus(tracker.fITSChannelStatus),
 269 fDetTypeRec(tracker.fDetTypeRec),
 270 fPlaneEff(tracker.fPlaneEff) {
 271   //Copy constructor
 272   Int_t i;
 273   for(i=0;i<4;i++) {
 274     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 275   }
 276   for(i=0;i<6;i++) {
 277     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 278     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 279   }
 280   for(i=0;i<2;i++) {
 281     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 282     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 283   }
 284 }
 285 //------------------------------------------------------------------------
 286 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
 287   //Assignment operator
 288   this->~AliITStrackerMI();
 289   new(this) AliITStrackerMI(tracker);
 290   return *this;
 291 }
 292 //------------------------------------------------------------------------
 293 AliITStrackerMI::~AliITStrackerMI()
 294 {
 295   //
 296   //destructor
 297   //
 298   if (fCoefficients) delete [] fCoefficients;
 299   DeleteTrksMaterialLUT();
 300   if (fDebugStreamer) {
 301     //fDebugStreamer->Close();
 302     delete fDebugStreamer;
 303   }
 304   if(fITSChannelStatus) delete fITSChannelStatus;
 305   if(fPlaneEff) delete fPlaneEff;
 306 }
 307 //------------------------------------------------------------------------
 308 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
 309   //--------------------------------------------------------------------
 310   //This function set masks of the layers which must be not skipped
 311   //--------------------------------------------------------------------
 312   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=l[i];
 313 }
 314 //------------------------------------------------------------------------
 315 void AliITStrackerMI::ReadBadFromDetTypeRec() {
 316   //--------------------------------------------------------------------
 317   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
 318   //i.e. from OCDB
 319   //--------------------------------------------------------------------
 320
 321   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
 322
 323   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
 324
 325   if(!fDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
 326
 327   // ITS channels map
 328   if(fITSChannelStatus) delete fITSChannelStatus;
 329   fITSChannelStatus = new AliITSChannelStatus(fDetTypeRec);
 330
 331   // ITS detectors and chips
 332   Int_t i=0,j=0,k=0,ndet=0;
 333   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 334     Int_t nBadDetsPerLayer=0;
 335     ndet=AliITSgeomTGeo::GetNDetectors(i);
 336     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
 337       for (k=1; k<ndet+1; k++) {
 338         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 339         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fDetTypeRec);
 340         if(det.IsBad()) {nBadDetsPerLayer++;}
 341       } // end loop on detectors
 342     } // end loop on ladders
 343     Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
 344   } // end loop on layers
 345
 346   return;
 347 }
 348 //------------------------------------------------------------------------
 349 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 350   //--------------------------------------------------------------------
 351   //This function loads ITS clusters
 352   //--------------------------------------------------------------------
 353   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 354   if (!branch) {
 355     Error("LoadClusters"," can't get the branch !\n");
 356     return 1;
 357   }
 358
 359   static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 360   branch->SetAddress(&clusters);
 361
 362   Int_t i=0,j=0,ndet=0;
 363   Int_t detector=0;
 364   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 365     ndet=fgLayers[i].GetNdetectors();
 366     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 367     for (; j<jmax; j++) {
 368       if (!cTree->GetEvent(j)) continue;
 369       Int_t ncl=clusters->GetEntriesFast();
 370       SignDeltas(clusters,GetZ());
 371
 372       while (ncl--) {
 373         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 374         detector=c->GetDetectorIndex();
 375
 376         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 377
 378         fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 379       }
 380       clusters->Delete();
 381       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 382       // zwindow cm from the dead zone
 383       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
 384         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
 385           Int_t lab[4]   = {0,0,0,detector};
 386           Int_t info[3]  = {0,0,i};
 387           Float_t q      = 0.; // this identifies virtual clusters
 388           Float_t hit[5] = {xdead,
 389                             0.,
 390                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
 391                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
 392                             q};
 393           Bool_t local   = kTRUE;
 394           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
 395           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 396           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 397             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 398           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 399           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 400             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 401           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 402           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 403             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 404           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 405           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 406             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 407           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 408           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 409             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 410           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 411           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 412             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 413         }
 414       } // "virtual" clusters in SPD
 415
 416     }
 417     //
 418     fgLayers[i].ResetRoad(); //road defined by the cluster density
 419     fgLayers[i].SortClusters();
 420   }
 421
 422   dummy.Clear();
 423
 424   return 0;
 425 }
 426 //------------------------------------------------------------------------
 427 void AliITStrackerMI::UnloadClusters() {
 428   //--------------------------------------------------------------------
 429   //This function unloads ITS clusters
 430   //--------------------------------------------------------------------
 431   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
 432 }
 433 //------------------------------------------------------------------------
 434 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
 435   //--------------------------------------------------------------------
 436   // Publishes all pointers to clusters known to the tracker into the
 437   // passed object array.
 438   // The ownership is not transfered - the caller is not expected to delete
 439   // the clusters.
 440   //--------------------------------------------------------------------
 441
 442   for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 443     for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
 444       AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
 445       array->AddLast(cl);
 446     }
 447   }
 448
 449   return;
 450 }
 451 //------------------------------------------------------------------------
 452 static Int_t CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
 453   //--------------------------------------------------------------------
 454   // Correction for the material between the TPC and the ITS
 455   //--------------------------------------------------------------------
 456   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 457       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 458       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 459       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 460   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 461       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 462       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 463       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 464   } else {
 465     Error("CorrectForTPCtoITSDeadZoneMaterial","Track is already in the dead zone !");
 466     return 0;
 467   }
 468
 469   return 1;
 470 }
 471 //------------------------------------------------------------------------
 472 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
 473   //--------------------------------------------------------------------
 474   // This functions reconstructs ITS tracks
 475   // The clusters must be already loaded !
 476   //--------------------------------------------------------------------
 477
 478
 479   fTrackingPhase="Clusters2Tracks";
 480
 481   TObjArray itsTracks(15000);
 482   fOriginal.Clear();
 483   fEsd = event;         // store pointer to the esd
 484
 485   // temporary (for cosmics)
 486   if(event->GetVertex()) {
 487     TString title = event->GetVertex()->GetTitle();
 488     if(title.Contains("cosmics")) {
 489       Double_t xyz[3]={GetX(),GetY(),GetZ()};
 490       Double_t exyz[3]={0.1,0.1,0.1};
 491       SetVertex(xyz,exyz);
 492     }
 493   }
 494   // temporary
 495
 496   {/* Read ESD tracks */
 497     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 498     Int_t nentr=event->GetNumberOfTracks();
 499     Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 500     while (nentr--) {
 501       AliESDtrack *esd=event->GetTrack(nentr);
 502       //  ---- for debugging:
 503       //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); }
 504
 505       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 506       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 507       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 508       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 509       AliITStrackMI *t=0;
 510       try {
 511         t=new AliITStrackMI(*esd);
 512       } catch (const Char_t *msg) {
 513         //Warning("Clusters2Tracks",msg);
 514         delete t;
 515         continue;
 516       }
 517       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
 518       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 519
 520
 521       // look at the ESD mass hypothesys !
 522       if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
 523       // write expected q
 524       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 525
 526       if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
 527         //track - can be  V0 according to TPC
 528       } else {
 529         if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
 530           delete t;
 531           continue;
 532         }
 533         if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
 534           delete t;
 535           continue;
 536         }
 537         if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
 538           delete t;
 539           continue;
 540         }
 541         if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 542           delete t;
 543           continue;
 544         }
 545       }
 546       t->SetReconstructed(kFALSE);
 547       itsTracks.AddLast(t);
 548       fOriginal.AddLast(t);
 549     }
 550   } /* End Read ESD tracks */
 551
 552   itsTracks.Sort();
 553   fOriginal.Sort();
 554   Int_t nentr=itsTracks.GetEntriesFast();
 555   fTrackHypothesys.Expand(nentr);
 556   fBestHypothesys.Expand(nentr);
 557   MakeCoefficients(nentr);
 558   if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
 559   Int_t ntrk=0;
 560   // THE TWO TRACKING PASSES
 561   for (fPass=0; fPass<2; fPass++) {
 562      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 563      for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
 564        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
 565        if (t==0) continue;              //this track has been already tracked
 566        //cout<<"========== "<<fPass<<"    "<<fCurrentEsdTrack<<" =========\n";
 567        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 568        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
 569        if (fConstraint[fPass]) {
 570          if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
 571              TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
 572        }
 573
 574        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 575        AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
 576        fI = 6;
 577        ResetTrackToFollow(*t);
 578        ResetBestTrack();
 579
 580        FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
 581
 582
 583        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 584        //
 585        AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 586        if (!besttrack) continue;
 587        besttrack->SetLabel(tpcLabel);
 588        //       besttrack->CookdEdx();
 589        CookdEdx(besttrack);
 590        besttrack->SetFakeRatio(1.);
 591        CookLabel(besttrack,0.); //For comparison only
 592        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 593
 594        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 595
 596        t->SetReconstructed(kTRUE);
 597        ntrk++;
 598        AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
 599      }
 600      GetBestHypothesysMIP(itsTracks);
 601   } // end loop on the two tracking passes
 602
 603   if(event->GetNumberOfV0s()>0) UpdateTPCV0(event);
 604   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) FindV02(event);
 605   fAfterV0 = kTRUE;
 606   //
 607   itsTracks.Delete();
 608   //
 609   Int_t entries = fTrackHypothesys.GetEntriesFast();
 610   for (Int_t ientry=0; ientry<entries; ientry++) {
 611     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
 612     if (array) array->Delete();
 613     delete fTrackHypothesys.RemoveAt(ientry);
 614   }
 615
 616   fTrackHypothesys.Delete();
 617   fBestHypothesys.Delete();
 618   fOriginal.Clear();
 619   delete [] fCoefficients;
 620   fCoefficients=0;
 621   DeleteTrksMaterialLUT();
 622
 623   Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
 624
 625   fTrackingPhase="Default";
 626
 627   return 0;
 628 }
 629 //------------------------------------------------------------------------
 630 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
 631   //--------------------------------------------------------------------
 632   // This functions propagates reconstructed ITS tracks back
 633   // The clusters must be loaded !
 634   //--------------------------------------------------------------------
 635   fTrackingPhase="PropagateBack";
 636   Int_t nentr=event->GetNumberOfTracks();
 637   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 638
 639   Int_t ntrk=0;
 640   for (Int_t i=0; i<nentr; i++) {
 641      AliESDtrack *esd=event->GetTrack(i);
 642
 643      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 644      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 645
 646      AliITStrackMI *t=0;
 647      try {
 648         t=new AliITStrackMI(*esd);
 649      } catch (const Char_t *msg) {
 650        //Warning("PropagateBack",msg);
 651         delete t;
 652         continue;
 653      }
 654      t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 655
 656      ResetTrackToFollow(*t);
 657
 658      // propagate to vertex [SR, GSI 17.02.2003]
 659      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 660      if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 661        if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
 662          fTrackToFollow.StartTimeIntegral();
 663        // from vertex to outside pipe
 664        CorrectForPipeMaterial(&fTrackToFollow,"outward");
 665      }
 666
 667      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
 668      if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
 669         if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
 670           delete t;
 671           continue;
 672         }
 673         fTrackToFollow.SetLabel(t->GetLabel());
 674         //fTrackToFollow.CookdEdx();
 675         CookLabel(&fTrackToFollow,0.); //For comparison only
 676         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 677         //UseClusters(&fTrackToFollow);
 678         ntrk++;
 679      }
 680      delete t;
 681   }
 682
 683   Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
 684
 685   fTrackingPhase="Default";
 686
 687   return 0;
 688 }
 689 //------------------------------------------------------------------------
 690 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
 691   //--------------------------------------------------------------------
 692   // This functions refits ITS tracks using the
 693   // "inward propagated" TPC tracks
 694   // The clusters must be loaded !
 695   //--------------------------------------------------------------------
 696   fTrackingPhase="RefitInward";
 697   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) RefitV02(event);
 698   Int_t nentr=event->GetNumberOfTracks();
 699   Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 700
 701   Int_t ntrk=0;
 702   for (Int_t i=0; i<nentr; i++) {
 703     AliESDtrack *esd=event->GetTrack(i);
 704
 705     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 706     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 707     if (esd->GetStatus()&AliESDtrack::kTPCout)
 708       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 709
 710     AliITStrackMI *t=0;
 711     try {
 712         t=new AliITStrackMI(*esd);
 713     } catch (const Char_t *msg) {
 714       //Warning("RefitInward",msg);
 715         delete t;
 716         continue;
 717     }
 718     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 719     if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 720        delete t;
 721        continue;
 722     }
 723
 724     ResetTrackToFollow(*t);
 725     fTrackToFollow.ResetClusters();
 726
 727     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
 728       fTrackToFollow.ResetCovariance(10.);
 729
 730     //Refitting...
 731     Bool_t pe=AliITSReconstructor::GetRecoParam()->GetComputePlaneEff();
 732     AliDebug(2,Form("Refit LABEL %d  %d",t->GetLabel(),t->GetNumberOfClusters()));
 733     if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
 734        AliDebug(2,"  refit OK");
 735        fTrackToFollow.SetLabel(t->GetLabel());
 736        //       fTrackToFollow.CookdEdx();
 737        CookdEdx(&fTrackToFollow);
 738
 739        CookLabel(&fTrackToFollow,0.0); //For comparison only
 740
 741        //The beam pipe
 742        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 743          fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
 744          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 745          //printf("                                       %d\n",esdTrack->GetITSModuleIndex(0));
 746          //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
 747          Float_t r[3]={0.,0.,0.};
 748          Double_t maxD=3.;
 749          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 750          ntrk++;
 751        }
 752     }
 753     delete t;
 754   }
 755
 756   Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
 757
 758   fTrackingPhase="Default";
 759
 760   return 0;
 761 }
 762 //------------------------------------------------------------------------
 763 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 764   //--------------------------------------------------------------------
 765   //       Return pointer to a given cluster
 766   //--------------------------------------------------------------------
 767   Int_t l=(index & 0xf0000000) >> 28;
 768   Int_t c=(index & 0x0fffffff) >> 00;
 769   return fgLayers[l].GetCluster(c);
 770 }
 771 //------------------------------------------------------------------------
 772 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 773   //--------------------------------------------------------------------
 774   // Get track space point with index i
 775   //--------------------------------------------------------------------
 776
 777   Int_t l=(index & 0xf0000000) >> 28;
 778   Int_t c=(index & 0x0fffffff) >> 00;
 779   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 780   Int_t idet = cl->GetDetectorIndex();
 781
 782   Float_t xyz[3];
 783   Float_t cov[6];
 784   cl->GetGlobalXYZ(xyz);
 785   cl->GetGlobalCov(cov);
 786   p.SetXYZ(xyz, cov);
 787   p.SetCharge(cl->GetQ());
 788   p.SetDriftTime(cl->GetDriftTime());
 789   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 790   switch (l) {
 791   case 0:
 792     iLayer = AliGeomManager::kSPD1;
 793     break;
 794   case 1:
 795     iLayer = AliGeomManager::kSPD2;
 796     break;
 797   case 2:
 798     iLayer = AliGeomManager::kSDD1;
 799     break;
 800   case 3:
 801     iLayer = AliGeomManager::kSDD2;
 802     break;
 803   case 4:
 804     iLayer = AliGeomManager::kSSD1;
 805     break;
 806   case 5:
 807     iLayer = AliGeomManager::kSSD2;
 808     break;
 809   default:
 810     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 811     break;
 812   };
 813   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 814   p.SetVolumeID((UShort_t)volid);
 815   return kTRUE;
 816 }
 817 //------------------------------------------------------------------------
 818 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
 819                         AliTrackPoint& p, const AliESDtrack *t) {
 820   //--------------------------------------------------------------------
 821   // Get track space point with index i
 822   // (assign error estimated during the tracking)
 823   //--------------------------------------------------------------------
 824
 825   Int_t l=(index & 0xf0000000) >> 28;
 826   Int_t c=(index & 0x0fffffff) >> 00;
 827   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 828   Int_t idet = cl->GetDetectorIndex();
 829
 830   const AliITSdetector &det=fgLayers[l].GetDetector(idet);
 831
 832   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 833   Float_t detxy[2];
 834   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 835   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 836   Double_t alpha = t->GetAlpha();
 837   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 838   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
 839   phi += alpha-det.GetPhi();
 840   Float_t tgphi = TMath::Tan(phi);
 841
 842   Float_t tgl = t->GetTgl(); // tgl about const along track
 843   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 844
 845   Float_t errlocalx,errlocalz;
 846   Bool_t addMisalErr=kFALSE;
 847   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
 848
 849   Float_t xyz[3];
 850   Float_t cov[6];
 851   cl->GetGlobalXYZ(xyz);
 852   //  cl->GetGlobalCov(cov);
 853   Float_t pos[3] = {0.,0.,0.};
 854   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
 855   tmpcl.GetGlobalCov(cov);
 856
 857   p.SetXYZ(xyz, cov);
 858   p.SetCharge(cl->GetQ());
 859   p.SetDriftTime(cl->GetDriftTime());
 860
 861   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 862   switch (l) {
 863   case 0:
 864     iLayer = AliGeomManager::kSPD1;
 865     break;
 866   case 1:
 867     iLayer = AliGeomManager::kSPD2;
 868     break;
 869   case 2:
 870     iLayer = AliGeomManager::kSDD1;
 871     break;
 872   case 3:
 873     iLayer = AliGeomManager::kSDD2;
 874     break;
 875   case 4:
 876     iLayer = AliGeomManager::kSSD1;
 877     break;
 878   case 5:
 879     iLayer = AliGeomManager::kSSD2;
 880     break;
 881   default:
 882     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 883     break;
 884   };
 885   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 886
 887   p.SetVolumeID((UShort_t)volid);
 888   return kTRUE;
 889 }
 890 //------------------------------------------------------------------------
 891 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
 892 {
 893   //--------------------------------------------------------------------
 894   // Follow prolongation tree
 895   //--------------------------------------------------------------------
 896   //
 897   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
 898   Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 899
 900
 901   AliESDtrack * esd = otrack->GetESDtrack();
 902   if (esd->GetV0Index(0)>0) {
 903     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
 904     //                      mapping of ESD track is different as ITS track in Containers
 905     //                      Need something more stable
 906     //                      Indexes are set back again to the ESD track indexes in UpdateTPCV0
 907     for (Int_t i=0;i<3;i++){
 908       Int_t  index = esd->GetV0Index(i);
 909       if (index==0) break;
 910       AliESDv0 * vertex = fEsd->GetV0(index);
 911       if (vertex->GetStatus()<0) continue;     // rejected V0
 912       //
 913       if (esd->GetSign()>0) {
 914         vertex->SetIndex(0,esdindex);
 915       } else {
 916         vertex->SetIndex(1,esdindex);
 917       }
 918     }
 919   }
 920   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
 921   if (!bestarray){
 922     bestarray = new TObjArray(5);
 923     fBestHypothesys.AddAt(bestarray,esdindex);
 924   }
 925
 926   //
 927   //setup tree of the prolongations
 928   //
 929   static AliITStrackMI tracks[7][100];
 930   AliITStrackMI *currenttrack;
 931   static AliITStrackMI currenttrack1;
 932   static AliITStrackMI currenttrack2;
 933   static AliITStrackMI backuptrack;
 934   Int_t ntracks[7];
 935   Int_t nindexes[7][100];
 936   Float_t normalizedchi2[100];
 937   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
 938   otrack->SetNSkipped(0);
 939   new (&(tracks[6][0])) AliITStrackMI(*otrack);
 940   ntracks[6]=1;
 941   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
 942   Int_t modstatus = 1; // found
 943   Float_t xloc,zloc;
 944   //
 945   //
 946   // follow prolongations
 947   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
 948     AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
 949     fI = ilayer;
 950     //
 951     AliITSlayer &layer=fgLayers[ilayer];
 952     Double_t r = layer.GetR();
 953     ntracks[ilayer]=0;
 954     //
 955     //
 956     Int_t nskipped=0;
 957     Float_t nused =0;
 958     for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
 959       //set current track
 960       if (ntracks[ilayer]>=100) break;
 961       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
 962       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
 963       if (ntracks[ilayer]>15+ilayer){
 964         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
 965         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
 966       }
 967
 968       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
 969
 970       // material between SSD and SDD, SDD and SPD
 971       if (ilayer==3)
 972         if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
 973       if (ilayer==1)
 974         if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
 975
 976       // detector number
 977       Double_t phi,z;
 978       if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
 979       Int_t idet=layer.FindDetectorIndex(phi,z);
 980
 981       Double_t trackGlobXYZ1[3];
 982       if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
 983
 984       // Get the budget to the primary vertex for the current track being prolonged
 985       Double_t budgetToPrimVertex = GetEffectiveThickness();
 986
 987       // check if we allow a prolongation without point
 988       Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
 989       if (skip) {
 990         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
 991         // propagate to the layer radius
 992         Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
 993         if(!vtrack->Propagate(xToGo)) continue;
 994         // apply correction for material of the current layer
 995         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
 996         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
 997         vtrack->SetClIndex(ilayer,0);
 998         modstatus = (skip==1 ? 3 : 4); // skipped : out in z
 999         if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1000           vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1001         }
1002         if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1003         ntracks[ilayer]++;
1004         continue;
1005       }
1006
1007       // track outside layer acceptance in z
1008       if (idet<0) continue;
1009
1010       //propagate to the intersection with the detector plane
1011       const AliITSdetector &det=layer.GetDetector(idet);
1012       new(&currenttrack2)  AliITStrackMI(currenttrack1);
1013       if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1014       if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1015       currenttrack1.SetDetectorIndex(idet);
1016       currenttrack2.SetDetectorIndex(idet);
1017       if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
1018
1019       //***************
1020       // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1021       //
1022       // road in global (rphi,z) [i.e. in tracking ref. system]
1023       Double_t zmin,zmax,ymin,ymax;
1024       if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1025
1026       // select clusters in road
1027       layer.SelectClusters(zmin,zmax,ymin,ymax);
1028       //********************
1029
1030       // Define criteria for track-cluster association
1031       Double_t msz = currenttrack1.GetSigmaZ2() +
1032         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1033         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1034         AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1035       Double_t msy = currenttrack1.GetSigmaY2() +
1036         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1037         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1038         AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1039       if (constrain) {
1040         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1041         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1042       }  else {
1043         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1044         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1045       }
1046       msz = 1./msz; // 1/RoadZ^2
1047       msy = 1./msy; // 1/RoadY^2
1048
1049       //
1050       //
1051       // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1052       //
1053       const AliITSRecPoint *cl=0;
1054       Int_t clidx=-1;
1055       Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1056       Bool_t deadzoneSPD=kFALSE;
1057       currenttrack = &currenttrack1;
1058
1059       // check if the road contains a dead zone
1060       Bool_t noClusters = kFALSE;
1061       if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1062       if (noClusters) AliDebug(2,"no clusters in road");
1063       Double_t dz=0.5*(zmax-zmin);
1064       Double_t dy=0.5*(ymax-ymin);
1065       Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters);
1066       if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1067       // create a prolongation without clusters (check also if there are no clusters in the road)
1068       if (dead ||
1069           (noClusters &&
1070            AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1071         AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1072         updatetrack->SetClIndex(ilayer,0);
1073         if (dead==0) {
1074           modstatus = 5; // no cls in road
1075         } else if (dead==1) {
1076           modstatus = 7; // holes in z in SPD
1077         } else if (dead==2 || dead==3) {
1078           modstatus = 2; // dead from OCDB
1079         }
1080         updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1081         // apply correction for material of the current layer
1082         CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1083         if (constrain) { // apply vertex constrain
1084           updatetrack->SetConstrain(constrain);
1085           Bool_t isPrim = kTRUE;
1086           if (ilayer<4) { // check that it's close to the vertex
1087             updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1088             if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1089                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1090                 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1091                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1092           }
1093           if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1094         }
1095         if (dead) {
1096           updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1097           if (dead==1) { // dead zone at z=0,+-7cm in SPD
1098             updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1099             deadzoneSPD=kTRUE;
1100           }
1101         }
1102         ntracks[ilayer]++;
1103       }
1104
1105       clidx=-1;
1106       // loop over clusters in the road
1107       while ((cl=layer.GetNextCluster(clidx))!=0) {
1108         if (ntracks[ilayer]>95) break; //space for skipped clusters
1109         Bool_t changedet =kFALSE;
1110         if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1111         Int_t idetc=cl->GetDetectorIndex();
1112
1113         if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1114           // take into account misalignment (bring track to real detector plane)
1115           Double_t xTrOrig = currenttrack->GetX();
1116           if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1117           // a first cut on track-cluster distance
1118           if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1119                (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1120             {  // cluster not associated to track
1121               AliDebug(2,"not associated");
1122               continue;
1123             }
1124           // bring track back to ideal detector plane
1125           if (!currenttrack->Propagate(xTrOrig)) continue;
1126         } else {                                      // have to move track to cluster's detector
1127           const AliITSdetector &detc=layer.GetDetector(idetc);
1128           // a first cut on track-cluster distance
1129           Double_t y;
1130           if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1131           if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1132                (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1133             continue; // cluster not associated to track
1134           //
1135           new (&backuptrack) AliITStrackMI(currenttrack2);
1136           changedet = kTRUE;
1137           currenttrack =&currenttrack2;
1138           if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1139             new (currenttrack) AliITStrackMI(backuptrack);
1140             changedet = kFALSE;
1141             continue;
1142           }
1143           currenttrack->SetDetectorIndex(idetc);
1144           // Get again the budget to the primary vertex
1145           // for the current track being prolonged, if had to change detector
1146           //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1147         }
1148
1149         // calculate track-clusters chi2
1150         chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1151         // chi2 cut
1152         AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1153         if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1154           if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1155           if (ntracks[ilayer]>=100) continue;
1156           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1157           updatetrack->SetClIndex(ilayer,0);
1158           if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
1159
1160           if (cl->GetQ()!=0) { // real cluster
1161             if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1162               AliDebug(2,"update failed");
1163               continue;
1164             }
1165             updatetrack->SetSampledEdx(cl->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
1166             modstatus = 1; // found
1167           } else {             // virtual cluster in dead zone
1168             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1169             updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1170             modstatus = 7; // holes in z in SPD
1171           }
1172
1173           if (changedet) {
1174             Float_t xlocnewdet,zlocnewdet;
1175             if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1176               updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1177             }
1178           } else {
1179             updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1180           }
1181           if (cl->IsUsed()) updatetrack->IncrementNUsed();
1182
1183           // apply correction for material of the current layer
1184           CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1185
1186           if (constrain) { // apply vertex constrain
1187             updatetrack->SetConstrain(constrain);
1188             Bool_t isPrim = kTRUE;
1189             if (ilayer<4) { // check that it's close to the vertex
1190               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1191               if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1192                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1193                   TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1194                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1195             }
1196             if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1197           } //apply vertex constrain
1198           ntracks[ilayer]++;
1199         }  // create new hypothesis
1200         else {
1201           AliDebug(2,"chi2 too large");
1202         }
1203
1204       } // loop over possible prolongations
1205
1206       // allow one prolongation without clusters
1207       if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1208         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1209         // apply correction for material of the current layer
1210         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1211         vtrack->SetClIndex(ilayer,0);
1212         modstatus = 3; // skipped
1213         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1214         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1215         vtrack->IncrementNSkipped();
1216         ntracks[ilayer]++;
1217       }
1218
1219       // allow one prolongation without clusters for tracks with |tgl|>1.1
1220       if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) {  //big theta - for low flux
1221         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1222         // apply correction for material of the current layer
1223         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1224         vtrack->SetClIndex(ilayer,0);
1225         modstatus = 3; // skipped
1226         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1227         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1228         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1229         ntracks[ilayer]++;
1230       }
1231
1232
1233     } // loop over tracks in layer ilayer+1
1234
1235     //loop over track candidates for the current layer
1236     //
1237     //
1238     Int_t accepted=0;
1239
1240     Int_t golden=0;
1241     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1242       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1243       if (normalizedchi2[itrack] <
1244           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1245       if (ilayer>4) {
1246         accepted++;
1247       } else {
1248         if (constrain) { // constrain
1249           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1250             accepted++;
1251         } else {        // !constrain
1252           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1253             accepted++;
1254         }
1255       }
1256     }
1257     // sort tracks by increasing normalized chi2
1258     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1259     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1260     if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1261     if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1262   } // end loop over layers
1263
1264
1265   //
1266   // Now select tracks to be kept
1267   //
1268   Int_t max = constrain ? 20 : 5;
1269
1270   // tracks that reach layer 0 (SPD inner)
1271   for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1272     AliITStrackMI & track= tracks[0][nindexes[0][i]];
1273     if (track.GetNumberOfClusters()<2) continue;
1274     if (!constrain && track.GetNormChi2(0) >
1275         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1276       continue;
1277     }
1278     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1279   }
1280
1281   // tracks that reach layer 1 (SPD outer)
1282   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1283     AliITStrackMI & track= tracks[1][nindexes[1][i]];
1284     if (track.GetNumberOfClusters()<4) continue;
1285     if (!constrain && track.GetNormChi2(1) >
1286         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1287     if (constrain) track.IncrementNSkipped();
1288     if (!constrain) {
1289       track.SetD(0,track.GetD(GetX(),GetY()));
1290       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1291       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1292         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1293       }
1294     }
1295     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1296   }
1297
1298   // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1299   if (!constrain){
1300     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1301       AliITStrackMI & track= tracks[2][nindexes[2][i]];
1302       if (track.GetNumberOfClusters()<3) continue;
1303       if (!constrain && track.GetNormChi2(2) >
1304           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1305       if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1306       if (!constrain){
1307         track.SetD(0,track.GetD(GetX(),GetY()));
1308         track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1309         if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1310           track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1311         }
1312       }
1313       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1314     }
1315   }
1316
1317   if (!constrain) {
1318     //
1319     // register best track of each layer - important for V0 finder
1320     //
1321     for (Int_t ilayer=0;ilayer<5;ilayer++){
1322       if (ntracks[ilayer]==0) continue;
1323       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1324       if (track.GetNumberOfClusters()<1) continue;
1325       CookLabel(&track,0);
1326       bestarray->AddAt(new AliITStrackMI(track),ilayer);
1327     }
1328   }
1329   //
1330   // update TPC V0 information
1331   //
1332   if (otrack->GetESDtrack()->GetV0Index(0)>0){
1333     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1334     for (Int_t i=0;i<3;i++){
1335       Int_t  index = otrack->GetESDtrack()->GetV0Index(i);
1336       if (index==0) break;
1337       AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1338       if (vertex->GetStatus()<0) continue;     // rejected V0
1339       //
1340       if (otrack->GetSign()>0) {
1341         vertex->SetIndex(0,esdindex);
1342       }
1343       else{
1344         vertex->SetIndex(1,esdindex);
1345       }
1346       //find nearest layer with track info
1347       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1348       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1349       Int_t nearest     = nearestold;
1350       for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1351         if (ntracks[nearest]==0){
1352           nearest = ilayer;
1353         }
1354       }
1355       //
1356       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1357       if (nearestold<5&&nearest<5){
1358         Bool_t accept = track.GetNormChi2(nearest)<10;
1359         if (accept){
1360           if (track.GetSign()>0) {
1361             vertex->SetParamP(track);
1362             vertex->Update(fprimvertex);
1363             //vertex->SetIndex(0,track.fESDtrack->GetID());
1364             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1365           }else{
1366             vertex->SetParamN(track);
1367             vertex->Update(fprimvertex);
1368             //vertex->SetIndex(1,track.fESDtrack->GetID());
1369             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1370           }
1371           vertex->SetStatus(vertex->GetStatus()+1);
1372         }else{
1373           //vertex->SetStatus(-2);  // reject V0  - not enough clusters
1374         }
1375       }
1376     }
1377   }
1378
1379 }
1380 //------------------------------------------------------------------------
1381 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1382 {
1383   //--------------------------------------------------------------------
1384   //
1385   //
1386   return fgLayers[layer];
1387 }
1388 //------------------------------------------------------------------------
1389 AliITStrackerMI::AliITSlayer::AliITSlayer():
1390 fR(0),
1391 fPhiOffset(0),
1392 fNladders(0),
1393 fZOffset(0),
1394 fNdetectors(0),
1395 fDetectors(0),
1396 fN(0),
1397 fDy5(0),
1398 fDy10(0),
1399 fDy20(0),
1400 fClustersCs(0),
1401 fClusterIndexCs(0),
1402 fYcs(0),
1403 fZcs(0),
1404 fNcs(0),
1405 fCurrentSlice(-1),
1406 fZmax(0),
1407 fYmin(0),
1408 fYmax(0),
1409 fI(0),
1410 fImax(0),
1411 fSkip(0),
1412 fAccepted(0),
1413 fRoad(0){
1414   //--------------------------------------------------------------------
1415   //default AliITSlayer constructor
1416   //--------------------------------------------------------------------
1417   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1418     fClusterWeight[i]=0;
1419     fClusterTracks[0][i]=-1;
1420     fClusterTracks[1][i]=-1;
1421     fClusterTracks[2][i]=-1;
1422     fClusterTracks[3][i]=-1;
1423   }
1424 }
1425 //------------------------------------------------------------------------
1426 AliITStrackerMI::AliITSlayer::
1427 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1428 fR(r),
1429 fPhiOffset(p),
1430 fNladders(nl),
1431 fZOffset(z),
1432 fNdetectors(nd),
1433 fDetectors(0),
1434 fN(0),
1435 fDy5(0),
1436 fDy10(0),
1437 fDy20(0),
1438 fClustersCs(0),
1439 fClusterIndexCs(0),
1440 fYcs(0),
1441 fZcs(0),
1442 fNcs(0),
1443 fCurrentSlice(-1),
1444 fZmax(0),
1445 fYmin(0),
1446 fYmax(0),
1447 fI(0),
1448 fImax(0),
1449 fSkip(0),
1450 fAccepted(0),
1451 fRoad(0) {
1452   //--------------------------------------------------------------------
1453   //main AliITSlayer constructor
1454   //--------------------------------------------------------------------
1455   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1456   fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1457 }
1458 //------------------------------------------------------------------------
1459 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1460 fR(layer.fR),
1461 fPhiOffset(layer.fPhiOffset),
1462 fNladders(layer.fNladders),
1463 fZOffset(layer.fZOffset),
1464 fNdetectors(layer.fNdetectors),
1465 fDetectors(layer.fDetectors),
1466 fN(layer.fN),
1467 fDy5(layer.fDy5),
1468 fDy10(layer.fDy10),
1469 fDy20(layer.fDy20),
1470 fClustersCs(layer.fClustersCs),
1471 fClusterIndexCs(layer.fClusterIndexCs),
1472 fYcs(layer.fYcs),
1473 fZcs(layer.fZcs),
1474 fNcs(layer.fNcs),
1475 fCurrentSlice(layer.fCurrentSlice),
1476 fZmax(layer.fZmax),
1477 fYmin(layer.fYmin),
1478 fYmax(layer.fYmax),
1479 fI(layer.fI),
1480 fImax(layer.fImax),
1481 fSkip(layer.fSkip),
1482 fAccepted(layer.fAccepted),
1483 fRoad(layer.fRoad){
1484   //Copy constructor
1485 }
1486 //------------------------------------------------------------------------
1487 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1488   //--------------------------------------------------------------------
1489   // AliITSlayer destructor
1490   //--------------------------------------------------------------------
1491   delete [] fDetectors;
1492   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1493   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1494     fClusterWeight[i]=0;
1495     fClusterTracks[0][i]=-1;
1496     fClusterTracks[1][i]=-1;
1497     fClusterTracks[2][i]=-1;
1498     fClusterTracks[3][i]=-1;
1499   }
1500 }
1501 //------------------------------------------------------------------------
1502 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1503   //--------------------------------------------------------------------
1504   // This function removes loaded clusters
1505   //--------------------------------------------------------------------
1506   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1507   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1508     fClusterWeight[i]=0;
1509     fClusterTracks[0][i]=-1;
1510     fClusterTracks[1][i]=-1;
1511     fClusterTracks[2][i]=-1;
1512     fClusterTracks[3][i]=-1;
1513   }
1514
1515   fN=0;
1516   fI=0;
1517 }
1518 //------------------------------------------------------------------------
1519 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1520   //--------------------------------------------------------------------
1521   // This function reset weights of the clusters
1522   //--------------------------------------------------------------------
1523   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1524     fClusterWeight[i]=0;
1525     fClusterTracks[0][i]=-1;
1526     fClusterTracks[1][i]=-1;
1527     fClusterTracks[2][i]=-1;
1528     fClusterTracks[3][i]=-1;
1529   }
1530   for (Int_t i=0; i<fN;i++) {
1531     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1532     if (cl&&cl->IsUsed()) cl->Use();
1533   }
1534
1535 }
1536 //------------------------------------------------------------------------
1537 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1538   //--------------------------------------------------------------------
1539   // This function calculates the road defined by the cluster density
1540   //--------------------------------------------------------------------
1541   Int_t n=0;
1542   for (Int_t i=0; i<fN; i++) {
1543      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1544   }
1545   if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1546 }
1547 //------------------------------------------------------------------------
1548 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1549   //--------------------------------------------------------------------
1550   //This function adds a cluster to this layer
1551   //--------------------------------------------------------------------
1552   if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1553     ::Error("InsertCluster","Too many clusters !\n");
1554     return 1;
1555   }
1556   fCurrentSlice=-1;
1557   fClusters[fN]=cl;
1558   fN++;
1559   AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1560   if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1561   if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1562   if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1563   if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1564
1565   return 0;
1566 }
1567 //------------------------------------------------------------------------
1568 void  AliITStrackerMI::AliITSlayer::SortClusters()
1569 {
1570   //
1571   //sort clusters
1572   //
1573   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1574   Float_t *z                = new Float_t[fN];
1575   Int_t   * index           = new Int_t[fN];
1576   //
1577   for (Int_t i=0;i<fN;i++){
1578     z[i] = fClusters[i]->GetZ();
1579   }
1580   TMath::Sort(fN,z,index,kFALSE);
1581   for (Int_t i=0;i<fN;i++){
1582     clusters[i] = fClusters[index[i]];
1583   }
1584   //
1585   for (Int_t i=0;i<fN;i++){
1586     fClusters[i] = clusters[i];
1587     fZ[i]        = fClusters[i]->GetZ();
1588     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1589     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1590     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1591     fY[i] = y;
1592   }
1593   delete[] index;
1594   delete[] z;
1595   delete[] clusters;
1596   //
1597
1598   fYB[0]=10000000;
1599   fYB[1]=-10000000;
1600   for (Int_t i=0;i<fN;i++){
1601     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1602     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1603     fClusterIndex[i] = i;
1604   }
1605   //
1606   // fill slices
1607   fDy5 = (fYB[1]-fYB[0])/5.;
1608   fDy10 = (fYB[1]-fYB[0])/10.;
1609   fDy20 = (fYB[1]-fYB[0])/20.;
1610   for (Int_t i=0;i<6;i++)  fN5[i] =0;
1611   for (Int_t i=0;i<11;i++) fN10[i]=0;
1612   for (Int_t i=0;i<21;i++) fN20[i]=0;
1613   //
1614   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;}
1615   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;}
1616   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;}
1617   //
1618   //
1619   for (Int_t i=0;i<fN;i++)
1620     for (Int_t irot=-1;irot<=1;irot++){
1621       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1622       // slice 5
1623       for (Int_t slice=0; slice<6;slice++){
1624         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1625           fClusters5[slice][fN5[slice]] = fClusters[i];
1626           fY5[slice][fN5[slice]] = curY;
1627           fZ5[slice][fN5[slice]] = fZ[i];
1628           fClusterIndex5[slice][fN5[slice]]=i;
1629           fN5[slice]++;
1630         }
1631       }
1632       // slice 10
1633       for (Int_t slice=0; slice<11;slice++){
1634         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1635           fClusters10[slice][fN10[slice]] = fClusters[i];
1636           fY10[slice][fN10[slice]] = curY;
1637           fZ10[slice][fN10[slice]] = fZ[i];
1638           fClusterIndex10[slice][fN10[slice]]=i;
1639           fN10[slice]++;
1640         }
1641       }
1642       // slice 20
1643       for (Int_t slice=0; slice<21;slice++){
1644         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1645           fClusters20[slice][fN20[slice]] = fClusters[i];
1646           fY20[slice][fN20[slice]] = curY;
1647           fZ20[slice][fN20[slice]] = fZ[i];
1648           fClusterIndex20[slice][fN20[slice]]=i;
1649           fN20[slice]++;
1650         }
1651       }
1652     }
1653
1654   //
1655   // consistency check
1656   //
1657   for (Int_t i=0;i<fN-1;i++){
1658     if (fZ[i]>fZ[i+1]){
1659       printf("Bug\n");
1660     }
1661   }
1662   //
1663   for (Int_t slice=0;slice<21;slice++)
1664   for (Int_t i=0;i<fN20[slice]-1;i++){
1665     if (fZ20[slice][i]>fZ20[slice][i+1]){
1666       printf("Bug\n");
1667     }
1668   }
1669
1670
1671 }
1672 //------------------------------------------------------------------------
1673 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1674   //--------------------------------------------------------------------
1675   // This function returns the index of the nearest cluster
1676   //--------------------------------------------------------------------
1677   Int_t ncl=0;
1678   const Float_t *zcl;
1679   if (fCurrentSlice<0) {
1680     ncl = fN;
1681     zcl   = fZ;
1682   }
1683   else{
1684     ncl   = fNcs;
1685     zcl   = fZcs;;
1686   }
1687
1688   if (ncl==0) return 0;
1689   Int_t b=0, e=ncl-1, m=(b+e)/2;
1690   for (; b<e; m=(b+e)/2) {
1691     //    if (z > fClusters[m]->GetZ()) b=m+1;
1692     if (z > zcl[m]) b=m+1;
1693     else e=m;
1694   }
1695   return m;
1696 }
1697 //------------------------------------------------------------------------
1698 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 {
1699   //--------------------------------------------------------------------
1700   // This function computes the rectangular road for this track
1701   //--------------------------------------------------------------------
1702
1703
1704   AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1705   // take into account the misalignment: propagate track to misaligned detector plane
1706   if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1707
1708   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1709                     TMath::Sqrt(track->GetSigmaZ2() +
1710                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1711                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1712                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1713   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1714                     TMath::Sqrt(track->GetSigmaY2() +
1715                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1716                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1717                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1718
1719   // track at boundary between detectors, enlarge road
1720   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1721   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1722        (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1723        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1724        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1725     Float_t tgl = TMath::Abs(track->GetTgl());
1726     if (tgl > 1.) tgl=1.;
1727     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1728     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1729     Float_t snp = TMath::Abs(track->GetSnp());
1730     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1731     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1732   } // boundary
1733
1734   // add to the road a term (up to 2-3 mm) to deal with misalignments
1735   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1736   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1737
1738   Double_t r = fgLayers[ilayer].GetR();
1739   zmin = track->GetZ() - dz;
1740   zmax = track->GetZ() + dz;
1741   ymin = track->GetY() + r*det.GetPhi() - dy;
1742   ymax = track->GetY() + r*det.GetPhi() + dy;
1743
1744   // bring track back to idead detector plane
1745   if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1746
1747   return kTRUE;
1748 }
1749 //------------------------------------------------------------------------
1750 void AliITStrackerMI::AliITSlayer::
1751 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1752   //--------------------------------------------------------------------
1753   // This function sets the "window"
1754   //--------------------------------------------------------------------
1755
1756   Double_t circle=2*TMath::Pi()*fR;
1757   fYmin = ymin; fYmax =ymax;
1758   Float_t ymiddle = (fYmin+fYmax)*0.5;
1759   if (ymiddle<fYB[0]) {
1760     fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1761   } else if (ymiddle>fYB[1]) {
1762     fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1763   }
1764
1765   //
1766   fCurrentSlice =-1;
1767   // defualt take all
1768   fClustersCs = fClusters;
1769   fClusterIndexCs = fClusterIndex;
1770   fYcs  = fY;
1771   fZcs  = fZ;
1772   fNcs  = fN;
1773   //
1774   //is in 20 slice?
1775   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1776     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1777     if (slice<0) slice=0;
1778     if (slice>20) slice=20;
1779     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1780     if (isOK) {
1781       fCurrentSlice=slice;
1782       fClustersCs = fClusters20[fCurrentSlice];
1783       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1784       fYcs  = fY20[fCurrentSlice];
1785       fZcs  = fZ20[fCurrentSlice];
1786       fNcs  = fN20[fCurrentSlice];
1787     }
1788   }
1789   //
1790   //is in 10 slice?
1791   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1792     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1793     if (slice<0) slice=0;
1794     if (slice>10) slice=10;
1795     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1796     if (isOK) {
1797       fCurrentSlice=slice;
1798       fClustersCs = fClusters10[fCurrentSlice];
1799       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1800       fYcs  = fY10[fCurrentSlice];
1801       fZcs  = fZ10[fCurrentSlice];
1802       fNcs  = fN10[fCurrentSlice];
1803     }
1804   }
1805   //
1806   //is in 5 slice?
1807   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1808     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1809     if (slice<0) slice=0;
1810     if (slice>5) slice=5;
1811     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1812     if (isOK) {
1813       fCurrentSlice=slice;
1814       fClustersCs = fClusters5[fCurrentSlice];
1815       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1816       fYcs  = fY5[fCurrentSlice];
1817       fZcs  = fZ5[fCurrentSlice];
1818       fNcs  = fN5[fCurrentSlice];
1819     }
1820   }
1821   //
1822   fI=FindClusterIndex(zmin); fZmax=zmax;
1823   fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1824   fSkip = 0;
1825   fAccepted =0;
1826
1827   return;
1828 }
1829 //------------------------------------------------------------------------
1830 Int_t AliITStrackerMI::AliITSlayer::
1831 FindDetectorIndex(Double_t phi, Double_t z) const {
1832   //--------------------------------------------------------------------
1833   //This function finds the detector crossed by the track
1834   //--------------------------------------------------------------------
1835   Double_t dphi;
1836   if (fZOffset<0)            // old geometry
1837     dphi = -(phi-fPhiOffset);
1838   else                       // new geometry
1839     dphi = phi-fPhiOffset;
1840
1841
1842   if      (dphi <  0) dphi += 2*TMath::Pi();
1843   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1844   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1845   if (np>=fNladders) np-=fNladders;
1846   if (np<0)          np+=fNladders;
1847
1848
1849   Double_t dz=fZOffset-z;
1850   Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1851   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1852   if (nz>=fNdetectors) return -1;
1853   if (nz<0)            return -1;
1854
1855   // ad hoc correction for 3rd ladder of SDD inner layer,
1856   // which is reversed (rotated by pi around local y)
1857   // this correction is OK only from AliITSv11Hybrid onwards
1858   if (GetR()>12. && GetR()<20.) { // SDD inner
1859     if(np==2) { // 3rd ladder
1860       nz = (fNdetectors-1) - nz;
1861     }
1862   }
1863   //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
1864
1865
1866   return np*fNdetectors + nz;
1867 }
1868 //------------------------------------------------------------------------
1869 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1870 {
1871   //--------------------------------------------------------------------
1872   // This function returns clusters within the "window"
1873   //--------------------------------------------------------------------
1874
1875   if (fCurrentSlice<0) {
1876     Double_t rpi2 = 2.*fR*TMath::Pi();
1877     for (Int_t i=fI; i<fImax; i++) {
1878       Double_t y = fY[i];
1879       if (fYmax<y) y -= rpi2;
1880       if (fYmin>y) y += rpi2;
1881       if (y<fYmin) continue;
1882       if (y>fYmax) continue;
1883       if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1884       ci=i;
1885       if (!test) fI=i+1;
1886       return fClusters[i];
1887     }
1888   } else {
1889     for (Int_t i=fI; i<fImax; i++) {
1890       if (fYcs[i]<fYmin) continue;
1891       if (fYcs[i]>fYmax) continue;
1892       if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1893       ci=fClusterIndexCs[i];
1894       if (!test) fI=i+1;
1895       return fClustersCs[i];
1896     }
1897   }
1898   return 0;
1899 }
1900 //------------------------------------------------------------------------
1901 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1902 const {
1903   //--------------------------------------------------------------------
1904   // This function returns the layer thickness at this point (units X0)
1905   //--------------------------------------------------------------------
1906   Double_t d=0.0085;
1907   x0=AliITSRecoParam::GetX0Air();
1908   if (43<fR&&fR<45) { //SSD2
1909      Double_t dd=0.0034;
1910      d=dd;
1911      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1912      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1913      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1914      for (Int_t i=0; i<12; i++) {
1915        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1916           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1917           d+=0.0034;
1918           break;
1919        }
1920        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1921           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1922           d+=0.0034;
1923           break;
1924        }
1925        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1926        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1927      }
1928   } else
1929   if (37<fR&&fR<41) { //SSD1
1930      Double_t dd=0.0034;
1931      d=dd;
1932      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1933      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1934      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1935      for (Int_t i=0; i<11; i++) {
1936        if (TMath::Abs(z-3.9*i)<0.15) {
1937           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1938           d+=dd;
1939           break;
1940        }
1941        if (TMath::Abs(z+3.9*i)<0.15) {
1942           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1943           d+=dd;
1944           break;
1945        }
1946        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1947        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1948      }
1949   } else
1950   if (13<fR&&fR<26) { //SDD
1951      Double_t dd=0.0033;
1952      d=dd;
1953      if (TMath::Abs(y-0.00)>3.30) d+=dd;
1954
1955      if (TMath::Abs(y-1.80)<0.55) {
1956         d+=0.016;
1957         for (Int_t j=0; j<20; j++) {
1958           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1959           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1960         }
1961      }
1962      if (TMath::Abs(y+1.80)<0.55) {
1963         d+=0.016;
1964         for (Int_t j=0; j<20; j++) {
1965           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1966           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1967         }
1968      }
1969
1970      for (Int_t i=0; i<4; i++) {
1971        if (TMath::Abs(z-7.3*i)<0.60) {
1972           d+=dd;
1973           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1974           break;
1975        }
1976        if (TMath::Abs(z+7.3*i)<0.60) {
1977           d+=dd;
1978           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1979           break;
1980        }
1981      }
1982   } else
1983   if (6<fR&&fR<8) {   //SPD2
1984      Double_t dd=0.0063; x0=21.5;
1985      d=dd;
1986      if (TMath::Abs(y-3.08)>0.5) d+=dd;
1987      if (TMath::Abs(y-3.03)<0.10) d+=0.014;
1988   } else
1989   if (3<fR&&fR<5) {   //SPD1
1990      Double_t dd=0.0063; x0=21.5;
1991      d=dd;
1992      if (TMath::Abs(y+0.21)>0.6) d+=dd;
1993      if (TMath::Abs(y+0.10)<0.10) d+=0.014;
1994   }
1995
1996   return d;
1997 }
1998 //------------------------------------------------------------------------
1999 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2000 fR(det.fR),
2001 fRmisal(det.fRmisal),
2002 fPhi(det.fPhi),
2003 fSinPhi(det.fSinPhi),
2004 fCosPhi(det.fCosPhi),
2005 fYmin(det.fYmin),
2006 fYmax(det.fYmax),
2007 fZmin(det.fZmin),
2008 fZmax(det.fZmax),
2009 fIsBad(det.fIsBad),
2010 fNChips(det.fNChips),
2011 fChipIsBad(det.fChipIsBad)
2012 {
2013   //Copy constructor
2014 }
2015 //------------------------------------------------------------------------
2016 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2017                                                  AliITSDetTypeRec *detTypeRec)
2018 {
2019   //--------------------------------------------------------------------
2020   // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2021   //--------------------------------------------------------------------
2022
2023   // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2024   // while in the tracker they start from 0 for each layer
2025   for(Int_t il=0; il<ilayer; il++)
2026     idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2027
2028   Int_t detType;
2029   if (ilayer==0 || ilayer==1) {        // ----------  SPD
2030     detType = 0;
2031   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
2032     detType = 1;
2033   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
2034     detType = 2;
2035   } else {
2036     printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2037     return;
2038   }
2039
2040   // Get calibration from AliITSDetTypeRec
2041   AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2042   calib->SetModuleIndex(idet);
2043   AliITSCalibration *calibSPDdead = 0;
2044   if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2045   if (calib->IsBad() ||
2046       (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2047     {
2048       SetBad();
2049       //      printf("lay %d bad %d\n",ilayer,idet);
2050     }
2051
2052   // Get segmentation from AliITSDetTypeRec
2053   AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2054
2055   // Read info about bad chips
2056   fNChips = segm->GetMaximumChipIndex()+1;
2057   //printf("ilayer %d  detType %d idet %d fNChips %d %d  GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2058   if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2059   fChipIsBad = new Bool_t[fNChips];
2060   for (Int_t iCh=0;iCh<fNChips;iCh++) {
2061     fChipIsBad[iCh] = calib->IsChipBad(iCh);
2062     if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2063     //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2064   }
2065
2066   return;
2067 }
2068 //------------------------------------------------------------------------
2069 Double_t AliITStrackerMI::GetEffectiveThickness()
2070 {
2071   //--------------------------------------------------------------------
2072   // Returns the thickness between the current layer and the vertex (units X0)
2073   //--------------------------------------------------------------------
2074
2075   if(fUseTGeo!=0) {
2076     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2077     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2078     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2079   }
2080
2081   // beam pipe
2082   Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2083   Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2084
2085   // layers
2086   Double_t x0=0;
2087   Double_t xn=fgLayers[fI].GetR();
2088   for (Int_t i=0; i<fI; i++) {
2089     Double_t xi=fgLayers[i].GetR();
2090     Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2091     d+=dLayer*xi*xi;
2092   }
2093
2094   // shields
2095   if (fI>1) {
2096     Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2097     d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2098   }
2099   if (fI>3) {
2100     Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2101     d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2102   }
2103   return d/(xn*xn);
2104 }
2105 //------------------------------------------------------------------------
2106 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2107   //-------------------------------------------------------------------
2108   // This function returns number of clusters within the "window"
2109   //--------------------------------------------------------------------
2110   Int_t ncl=0;
2111   for (Int_t i=fI; i<fN; i++) {
2112     const AliITSRecPoint *c=fClusters[i];
2113     if (c->GetZ() > fZmax) break;
2114     if (c->IsUsed()) continue;
2115     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2116     Double_t y=fR*det.GetPhi() + c->GetY();
2117
2118     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2119     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2120
2121     if (y<fYmin) continue;
2122     if (y>fYmax) continue;
2123     ncl++;
2124   }
2125   return ncl;
2126 }
2127 //------------------------------------------------------------------------
2128 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2129                                 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2130 {
2131   //--------------------------------------------------------------------
2132   // This function refits the track "track" at the position "x" using
2133   // the clusters from "clusters"
2134   // If "extra"==kTRUE,
2135   //    the clusters from overlapped modules get attached to "track"
2136   // If "planeff"==kTRUE,
2137   //    special approach for plane efficiency evaluation is applyed
2138   //--------------------------------------------------------------------
2139
2140   Int_t index[AliITSgeomTGeo::kNLayers];
2141   Int_t k;
2142   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2143   Int_t nc=clusters->GetNumberOfClusters();
2144   for (k=0; k<nc; k++) {
2145     Int_t idx=clusters->GetClusterIndex(k);
2146     Int_t ilayer=(idx&0xf0000000)>>28;
2147     index[ilayer]=idx;
2148   }
2149
2150   return RefitAt(xx,track,index,extra,planeeff); // call the method below
2151 }
2152 //------------------------------------------------------------------------
2153 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2154                                 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2155 {
2156   //--------------------------------------------------------------------
2157   // This function refits the track "track" at the position "x" using
2158   // the clusters from array
2159   // If "extra"==kTRUE,
2160   //    the clusters from overlapped modules get attached to "track"
2161   // If "planeff"==kTRUE,
2162   //    special approach for plane efficiency evaluation is applyed
2163   //--------------------------------------------------------------------
2164   Int_t index[AliITSgeomTGeo::kNLayers];
2165   Int_t k;
2166   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2167   //
2168   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2169     index[k]=clusters[k];
2170   }
2171
2172   // special for cosmics: check which the innermost layer crossed
2173   // by the track
2174   Int_t innermostlayer=5;
2175   Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2176   for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2177     if(drphi < fgLayers[innermostlayer].GetR()) break;
2178   }
2179   //printf(" drphi  %f  innermost %d\n",drphi,innermostlayer);
2180
2181   Int_t modstatus=1; // found
2182   Float_t xloc,zloc;
2183   Int_t from, to, step;
2184   if (xx > track->GetX()) {
2185       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2186       step=+1;
2187   } else {
2188       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2189       step=-1;
2190   }
2191   TString dir = (step>0 ? "outward" : "inward");
2192
2193   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2194      AliITSlayer &layer=fgLayers[ilayer];
2195      Double_t r=layer.GetR();
2196      if (step<0 && xx>r) break;
2197
2198      // material between SSD and SDD, SDD and SPD
2199      Double_t hI=ilayer-0.5*step;
2200      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2201        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2202      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2203        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2204
2205      // remember old position [SR, GSI 18.02.2003]
2206      Double_t oldX=0., oldY=0., oldZ=0.;
2207      if (track->IsStartedTimeIntegral() && step==1) {
2208         if (!track->GetGlobalXYZat(track->GetX(),oldX,oldY,oldZ)) return kFALSE;
2209      }
2210      //
2211
2212      Double_t oldGlobXYZ[3];
2213      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2214      //TMath::Sqrt(track->GetSigmaY2());
2215
2216      Double_t phi,z;
2217      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2218
2219      Int_t idet=layer.FindDetectorIndex(phi,z);
2220
2221      // check if we allow a prolongation without point for large-eta tracks
2222      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2223      if (skip==2) {
2224        // propagate to the layer radius
2225        Double_t xToGo; if (!track->GetLocalXat(r,xToGo)) return kFALSE;
2226        if (!track->Propagate(xToGo)) return kFALSE;
2227        // apply correction for material of the current layer
2228        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
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        // track time update [SR, GSI 17.02.2003]
2234        if (track->IsStartedTimeIntegral() && step==1) {
2235          Double_t newX, newY, newZ;
2236          if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2237          Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2238                         (oldZ-newZ)*(oldZ-newZ);
2239          track->AddTimeStep(TMath::Sqrt(dL2));
2240        }
2241        continue;
2242      }
2243
2244      if (idet<0) return kFALSE;
2245
2246      const AliITSdetector &det=layer.GetDetector(idet);
2247      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2248
2249      track->SetDetectorIndex(idet);
2250      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2251
2252      Double_t dz,zmin,zmax,dy,ymin,ymax;
2253
2254      const AliITSRecPoint *clAcc=0;
2255      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2256
2257      Int_t idx=index[ilayer];
2258      if (idx>=0) { // cluster in this layer
2259        modstatus = 6; // no refit
2260        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2261        if (cl) {
2262          if (idet != cl->GetDetectorIndex()) {
2263            idet=cl->GetDetectorIndex();
2264            const AliITSdetector &detc=layer.GetDetector(idet);
2265            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2266            track->SetDetectorIndex(idet);
2267            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2268          }
2269          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2270          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2271          if (chi2<maxchi2) {
2272            clAcc=cl;
2273            maxchi2=chi2;
2274            modstatus = 1; // found
2275          } else {
2276             return kFALSE; //
2277          }
2278        }
2279      } else { // no cluster in this layer
2280        if (skip==1) {
2281          modstatus = 3; // skipped
2282          // Plane Eff determination:
2283          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2284            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2285               UseTrackForPlaneEff(track,ilayer);
2286          }
2287        } else {
2288          modstatus = 5; // no cls in road
2289          // check dead
2290          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2291          dz = 0.5*(zmax-zmin);
2292          dy = 0.5*(ymax-ymin);
2293          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2294          if (dead==1) modstatus = 7; // holes in z in SPD
2295          if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2296        }
2297      }
2298
2299      if (clAcc) {
2300        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2301        track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2302      }
2303      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2304
2305
2306      if (extra) { // search for extra clusters in overlapped modules
2307        AliITStrackV2 tmp(*track);
2308        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2309        layer.SelectClusters(zmin,zmax,ymin,ymax);
2310
2311        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2312        Int_t idetExtra=-1;
2313        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2314        Double_t tolerance=0.1;
2315        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2316          // only clusters in another module! (overlaps)
2317          idetExtra = clExtra->GetDetectorIndex();
2318          if (idet == idetExtra) continue;
2319
2320          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2321
2322          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2323          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2324          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2325          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2326
2327          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2328          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2329        }
2330        if (cci>=0) {
2331          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2332          track->SetExtraModule(ilayer,idetExtra);
2333        }
2334      } // end search for extra clusters in overlapped modules
2335
2336      // Correct for material of the current layer
2337      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2338
2339      // track time update [SR, GSI 17.02.2003]
2340      if (track->IsStartedTimeIntegral() && step==1) {
2341         Double_t newX, newY, newZ;
2342         if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2343         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2344                        (oldZ-newZ)*(oldZ-newZ);
2345         track->AddTimeStep(TMath::Sqrt(dL2));
2346      }
2347      //
2348
2349   } // end loop on layers
2350
2351   if (!track->Propagate(xx)) return kFALSE;
2352
2353   return kTRUE;
2354 }
2355 //------------------------------------------------------------------------
2356 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2357 {
2358   //
2359   // calculate normalized chi2
2360   //  return NormalizedChi2(track,0);
2361   Float_t chi2 = 0;
2362   Float_t sum=0;
2363   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2364   //  track->fdEdxMismatch=0;
2365   Float_t dedxmismatch =0;
2366   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2367   if (mode<100){
2368     for (Int_t i = 0;i<6;i++){
2369       if (track->GetClIndex(i)>0){
2370         Float_t cerry, cerrz;
2371         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2372         else
2373           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2374         cerry*=cerry;
2375         cerrz*=cerrz;
2376         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2377         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2378           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2379           if (ratio<0.5) {
2380             cchi2+=(0.5-ratio)*10.;
2381             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2382             dedxmismatch+=(0.5-ratio)*10.;
2383           }
2384         }
2385         if (i<2 ||i>3){
2386           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2387           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2388           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2389           if (i<2) chi2+=2*cl->GetDeltaProbability();
2390         }
2391         chi2+=cchi2;
2392         sum++;
2393       }
2394     }
2395     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2396       track->SetdEdxMismatch(dedxmismatch);
2397     }
2398   }
2399   else{
2400     for (Int_t i = 0;i<4;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         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2408         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2409         sum++;
2410       }
2411     }
2412     for (Int_t i = 4;i<6;i++){
2413       if (track->GetClIndex(i)>0){
2414         Float_t cerry, cerrz;
2415         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2416         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2417         cerry*=cerry;
2418         cerrz*=cerrz;
2419         Float_t cerryb, cerrzb;
2420         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2421         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2422         cerryb*=cerryb;
2423         cerrzb*=cerrzb;
2424         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2425         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2426         sum++;
2427       }
2428     }
2429   }
2430   if (track->GetESDtrack()->GetTPCsignal()>85){
2431     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2432     if (ratio<0.5) {
2433       chi2+=(0.5-ratio)*5.;
2434     }
2435     if (ratio>2){
2436       chi2+=(ratio-2.0)*3;
2437     }
2438   }
2439   //
2440   Double_t match = TMath::Sqrt(track->GetChi22());
2441   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2442   if (!track->GetConstrain()) {
2443     if (track->GetNumberOfClusters()>2) {
2444       match/=track->GetNumberOfClusters()-2.;
2445     } else {
2446       match=0;
2447     }
2448   }
2449   if (match<0) match=0;
2450   Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2451   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2452     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2453                                 1./(1.+track->GetNSkipped()));
2454
2455  return normchi2;
2456 }
2457 //------------------------------------------------------------------------
2458 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2459 {
2460   //
2461   // return matching chi2 between two tracks
2462   Double_t largeChi2=1000.;
2463
2464   AliITStrackMI track3(*track2);
2465   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2466   TMatrixD vec(5,1);
2467   vec(0,0)=track1->GetY()   - track3.GetY();
2468   vec(1,0)=track1->GetZ()   - track3.GetZ();
2469   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2470   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2471   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2472   //
2473   TMatrixD cov(5,5);
2474   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2475   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2476   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2477   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2478   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2479
2480   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2481   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2482   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2483   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2484   //
2485   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2486   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2487   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2488   //
2489   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2490   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2491   //
2492   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2493
2494   cov.Invert();
2495   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2496   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2497   return chi2(0,0);
2498 }
2499 //------------------------------------------------------------------------
2500 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr)
2501 {
2502   //
2503   //  return probability that given point (characterized by z position and error)
2504   //  is in SPD dead zone
2505   //
2506   Double_t probability = 0.;
2507   Double_t absz = TMath::Abs(zpos);
2508   Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2509                                   0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2510   if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2511   Double_t zmin, zmax;
2512   if (zpos<-6.) { // dead zone at z = -7
2513     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2514     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2515   } else if (zpos>6.) { // dead zone at z = +7
2516     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2517     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2518   } else if (absz<2.) { // dead zone at z = 0
2519     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2520     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2521   } else {
2522     zmin = 0.;
2523     zmax = 0.;
2524   }
2525   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2526   // dead zone)
2527   probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2528                       TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2529   return probability;
2530 }
2531 //------------------------------------------------------------------------
2532 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2533 {
2534   //
2535   // calculate normalized chi2
2536   Float_t chi2[6];
2537   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2538   Float_t ncl = 0;
2539   for (Int_t i = 0;i<6;i++){
2540     if (TMath::Abs(track->GetDy(i))>0){
2541       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2542       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2543       ncl++;
2544     }
2545     else{chi2[i]=10000;}
2546   }
2547   Int_t index[6];
2548   TMath::Sort(6,chi2,index,kFALSE);
2549   Float_t max = float(ncl)*fac-1.;
2550   Float_t sumchi=0, sumweight=0;
2551   for (Int_t i=0;i<max+1;i++){
2552     Float_t weight = (i<max)?1.:(max+1.-i);
2553     sumchi+=weight*chi2[index[i]];
2554     sumweight+=weight;
2555   }
2556   Double_t normchi2 = sumchi/sumweight;
2557   return normchi2;
2558 }
2559 //------------------------------------------------------------------------
2560 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2561 {
2562   //
2563   // calculate normalized chi2
2564   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2565   Int_t npoints = 0;
2566   Double_t res =0;
2567   for (Int_t i=0;i<6;i++){
2568     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2569     Double_t sy1 = forwardtrack->GetSigmaY(i);
2570     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2571     Double_t sy2 = backtrack->GetSigmaY(i);
2572     Double_t sz2 = backtrack->GetSigmaZ(i);
2573     if (i<2){ sy2=1000.;sz2=1000;}
2574     //
2575     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2576     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2577     //
2578     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2579     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2580     //
2581     res+= nz0*nz0+ny0*ny0;
2582     npoints++;
2583   }
2584   if (npoints>1) return
2585                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2586                    //2*forwardtrack->fNUsed+
2587                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2588                                   1./(1.+forwardtrack->GetNSkipped()));
2589   return 1000;
2590 }
2591 //------------------------------------------------------------------------
2592 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2593   //--------------------------------------------------------------------
2594   //       Return pointer to a given cluster
2595   //--------------------------------------------------------------------
2596   Int_t l=(index & 0xf0000000) >> 28;
2597   Int_t c=(index & 0x0fffffff) >> 00;
2598   return fgLayers[l].GetWeight(c);
2599 }
2600 //------------------------------------------------------------------------
2601 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2602 {
2603   //---------------------------------------------
2604   // register track to the list
2605   //
2606   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2607   //
2608   //
2609   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2610     Int_t index = track->GetClusterIndex(icluster);
2611     Int_t l=(index & 0xf0000000) >> 28;
2612     Int_t c=(index & 0x0fffffff) >> 00;
2613     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2614     for (Int_t itrack=0;itrack<4;itrack++){
2615       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2616         fgLayers[l].SetClusterTracks(itrack,c,id);
2617         break;
2618       }
2619     }
2620   }
2621 }
2622 //------------------------------------------------------------------------
2623 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2624 {
2625   //---------------------------------------------
2626   // unregister track from the list
2627   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2628     Int_t index = track->GetClusterIndex(icluster);
2629     Int_t l=(index & 0xf0000000) >> 28;
2630     Int_t c=(index & 0x0fffffff) >> 00;
2631     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2632     for (Int_t itrack=0;itrack<4;itrack++){
2633       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2634         fgLayers[l].SetClusterTracks(itrack,c,-1);
2635       }
2636     }
2637   }
2638 }
2639 //------------------------------------------------------------------------
2640 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2641 {
2642   //-------------------------------------------------------------
2643   //get number of shared clusters
2644   //-------------------------------------------------------------
2645   Float_t shared=0;
2646   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2647   // mean  number of clusters
2648   Float_t *ny = GetNy(id), *nz = GetNz(id);
2649
2650
2651   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2652     Int_t index = track->GetClusterIndex(icluster);
2653     Int_t l=(index & 0xf0000000) >> 28;
2654     Int_t c=(index & 0x0fffffff) >> 00;
2655     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2656     if (ny[l]==0){
2657       printf("problem\n");
2658     }
2659     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2660     Float_t weight=1;
2661     //
2662     Float_t deltan = 0;
2663     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2664     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2665       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2666     if (l<2 || l>3){
2667       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2668     }
2669     else{
2670       deltan = (cl->GetNz()-nz[l]);
2671     }
2672     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2673     weight = 2./TMath::Max(3.+deltan,2.);
2674     //
2675     for (Int_t itrack=0;itrack<4;itrack++){
2676       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2677         list[l]=index;
2678         clist[l] = (AliITSRecPoint*)GetCluster(index);
2679         shared+=weight;
2680         break;
2681       }
2682     }
2683   }
2684   track->SetNUsed(shared);
2685   return shared;
2686 }
2687 //------------------------------------------------------------------------
2688 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2689 {
2690   //
2691   // find first shared track
2692   //
2693   // mean  number of clusters
2694   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2695   //
2696   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2697   Int_t sharedtrack=100000;
2698   Int_t tracks[24],trackindex=0;
2699   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2700   //
2701   for (Int_t icluster=0;icluster<6;icluster++){
2702     if (clusterlist[icluster]<0) continue;
2703     Int_t index = clusterlist[icluster];
2704     Int_t l=(index & 0xf0000000) >> 28;
2705     Int_t c=(index & 0x0fffffff) >> 00;
2706     if (ny[l]==0){
2707       printf("problem\n");
2708     }
2709     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2710     //if (l>3) continue;
2711     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2712     //
2713     Float_t deltan = 0;
2714     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2715     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2716       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2717     if (l<2 || l>3){
2718       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2719     }
2720     else{
2721       deltan = (cl->GetNz()-nz[l]);
2722     }
2723     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2724     //
2725     for (Int_t itrack=3;itrack>=0;itrack--){
2726       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2727       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2728        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2729        trackindex++;
2730       }
2731     }
2732   }
2733   if (trackindex==0) return -1;
2734   if (trackindex==1){
2735     sharedtrack = tracks[0];
2736   }else{
2737     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2738     else{
2739       //
2740       Int_t tracks2[24], cluster[24];
2741       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2742       Int_t index =0;
2743       //
2744       for (Int_t i=0;i<trackindex;i++){
2745         if (tracks[i]<0) continue;
2746         tracks2[index] = tracks[i];
2747         cluster[index]++;
2748         for (Int_t j=i+1;j<trackindex;j++){
2749           if (tracks[j]<0) continue;
2750           if (tracks[j]==tracks[i]){
2751             cluster[index]++;
2752             tracks[j]=-1;
2753           }
2754         }
2755         index++;
2756       }
2757       Int_t max=0;
2758       for (Int_t i=0;i<index;i++){
2759         if (cluster[index]>max) {
2760           sharedtrack=tracks2[index];
2761           max=cluster[index];
2762         }
2763       }
2764     }
2765   }
2766
2767   if (sharedtrack>=100000) return -1;
2768   //
2769   // make list of overlaps
2770   shared =0;
2771   for (Int_t icluster=0;icluster<6;icluster++){
2772     if (clusterlist[icluster]<0) continue;
2773     Int_t index = clusterlist[icluster];
2774     Int_t l=(index & 0xf0000000) >> 28;
2775     Int_t c=(index & 0x0fffffff) >> 00;
2776     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2777     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2778     if (l==0 || l==1){
2779       if (cl->GetNy()>2) continue;
2780       if (cl->GetNz()>2) continue;
2781     }
2782     if (l==4 || l==5){
2783       if (cl->GetNy()>3) continue;
2784       if (cl->GetNz()>3) continue;
2785     }
2786     //
2787     for (Int_t itrack=3;itrack>=0;itrack--){
2788       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2789       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2790         overlist[l]=index;
2791         shared++;
2792       }
2793     }
2794   }
2795   return sharedtrack;
2796 }
2797 //------------------------------------------------------------------------
2798 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2799   //
2800   // try to find track hypothesys without conflicts
2801   // with minimal chi2;
2802   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2803   Int_t entries1 = arr1->GetEntriesFast();
2804   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2805   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2806   Int_t entries2 = arr2->GetEntriesFast();
2807   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2808   //
2809   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2810   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2811   if (track10->Pt()>0.5+track20->Pt()) return track10;
2812
2813   for (Int_t itrack=0;itrack<entries1;itrack++){
2814     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2815     UnRegisterClusterTracks(track,trackID1);
2816   }
2817   //
2818   for (Int_t itrack=0;itrack<entries2;itrack++){
2819     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2820     UnRegisterClusterTracks(track,trackID2);
2821   }
2822   Int_t index1=0;
2823   Int_t index2=0;
2824   Float_t maxconflicts=6;
2825   Double_t maxchi2 =1000.;
2826   //
2827   // get weight of hypothesys - using best hypothesy
2828   Double_t w1,w2;
2829
2830   Int_t list1[6],list2[6];
2831   AliITSRecPoint *clist1[6], *clist2[6] ;
2832   RegisterClusterTracks(track10,trackID1);
2833   RegisterClusterTracks(track20,trackID2);
2834   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2835   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2836   UnRegisterClusterTracks(track10,trackID1);
2837   UnRegisterClusterTracks(track20,trackID2);
2838   //
2839   // normalized chi2
2840   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2841   Float_t nerry[6],nerrz[6];
2842   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2843   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2844   for (Int_t i=0;i<6;i++){
2845      if ( (erry1[i]>0) && (erry2[i]>0)) {
2846        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2847        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2848      }else{
2849        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2850        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2851      }
2852      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2853        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2854        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2855        ncl1++;
2856      }
2857      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2858        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2859        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2860        ncl2++;
2861      }
2862   }
2863   chi21/=ncl1;
2864   chi22/=ncl2;
2865   //
2866   //
2867   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2868   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2869   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2870   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2871   //
2872   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2873         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2874         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2875         );
2876   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2877         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2878         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2879         );
2880
2881   Double_t sumw = w1+w2;
2882   w1/=sumw;
2883   w2/=sumw;
2884   if (w1<w2*0.5) {
2885     w1 = (d2+0.5)/(d1+d2+1.);
2886     w2 = (d1+0.5)/(d1+d2+1.);
2887   }
2888   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2889   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2890   //
2891   // get pair of "best" hypothesys
2892   //
2893   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2894   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2895
2896   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2897     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2898     //if (track1->fFakeRatio>0) continue;
2899     RegisterClusterTracks(track1,trackID1);
2900     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2901       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2902
2903       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2904       //if (track2->fFakeRatio>0) continue;
2905       Float_t nskipped=0;
2906       RegisterClusterTracks(track2,trackID2);
2907       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2908       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2909       UnRegisterClusterTracks(track2,trackID2);
2910       //
2911       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2912       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2913       if (nskipped>0.5) continue;
2914       //
2915       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2916       if (conflict1+1<cconflict1) continue;
2917       if (conflict2+1<cconflict2) continue;
2918       Float_t conflict=0;
2919       Float_t sumchi2=0;
2920       Float_t sum=0;
2921       for (Int_t i=0;i<6;i++){
2922         //
2923         Float_t c1 =0.; // conflict coeficients
2924         Float_t c2 =0.;
2925         if (clist1[i]&&clist2[i]){
2926           Float_t deltan = 0;
2927           if (i<2 || i>3){
2928             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2929           }
2930           else{
2931             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2932           }
2933           c1  = 2./TMath::Max(3.+deltan,2.);
2934           c2  = 2./TMath::Max(3.+deltan,2.);
2935         }
2936         else{
2937           if (clist1[i]){
2938             Float_t deltan = 0;
2939             if (i<2 || i>3){
2940               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2941             }
2942             else{
2943               deltan = (clist1[i]->GetNz()-nz1[i]);
2944             }
2945             c1  = 2./TMath::Max(3.+deltan,2.);
2946             c2  = 0;
2947           }
2948
2949           if (clist2[i]){
2950             Float_t deltan = 0;
2951             if (i<2 || i>3){
2952               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2953             }
2954             else{
2955               deltan = (clist2[i]->GetNz()-nz2[i]);
2956             }
2957             c2  = 2./TMath::Max(3.+deltan,2.);
2958             c1  = 0;
2959           }
2960         }
2961         //
2962         chi21=0;chi22=0;
2963         if (TMath::Abs(track1->GetDy(i))>0.) {
2964           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2965             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2966           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2967           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2968         }else{
2969           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2970         }
2971         //
2972         if (TMath::Abs(track2->GetDy(i))>0.) {
2973           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2974             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2975           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2976           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2977         }
2978         else{
2979           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2980         }
2981         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2982         if (chi21>0) sum+=w1;
2983         if (chi22>0) sum+=w2;
2984         conflict+=(c1+c2);
2985       }
2986       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2987       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2988       Double_t normchi2 = 2*conflict+sumchi2/sum;
2989       if ( normchi2 <maxchi2 ){
2990         index1 = itrack1;
2991         index2 = itrack2;
2992         maxconflicts = conflict;
2993         maxchi2 = normchi2;
2994       }
2995     }
2996     UnRegisterClusterTracks(track1,trackID1);
2997   }
2998   //
2999   //  if (maxconflicts<4 && maxchi2<th0){
3000   if (maxchi2<th0*2.){
3001     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3002     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3003     track1->SetChi2MIP(5,maxconflicts);
3004     track1->SetChi2MIP(6,maxchi2);
3005     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3006     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3007     track1->SetChi2MIP(8,index1);
3008     fBestTrackIndex[trackID1] =index1;
3009     UpdateESDtrack(track1, AliESDtrack::kITSin);
3010   }
3011   else if (track10->GetChi2MIP(0)<th1){
3012     track10->SetChi2MIP(5,maxconflicts);
3013     track10->SetChi2MIP(6,maxchi2);
3014     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3015     UpdateESDtrack(track10,AliESDtrack::kITSin);
3016   }
3017
3018   for (Int_t itrack=0;itrack<entries1;itrack++){
3019     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3020     UnRegisterClusterTracks(track,trackID1);
3021   }
3022   //
3023   for (Int_t itrack=0;itrack<entries2;itrack++){
3024     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3025     UnRegisterClusterTracks(track,trackID2);
3026   }
3027
3028   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3029       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3030     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3031   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3032     RegisterClusterTracks(track10,trackID1);
3033   }
3034   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3035       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3036     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3037     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3038     RegisterClusterTracks(track20,trackID2);
3039   }
3040   return track10;
3041
3042 }
3043 //------------------------------------------------------------------------
3044 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3045   //--------------------------------------------------------------------
3046   // This function marks clusters assigned to the track
3047   //--------------------------------------------------------------------
3048   AliTracker::UseClusters(t,from);
3049
3050   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3051   //if (c->GetQ()>2) c->Use();
3052   if (c->GetSigmaZ2()>0.1) c->Use();
3053   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3054   //if (c->GetQ()>2) c->Use();
3055   if (c->GetSigmaZ2()>0.1) c->Use();
3056
3057 }
3058 //------------------------------------------------------------------------
3059 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3060 {
3061   //------------------------------------------------------------------
3062   // add track to the list of hypothesys
3063   //------------------------------------------------------------------
3064
3065   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3066     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3067   //
3068   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3069   if (!array) {
3070     array = new TObjArray(10);
3071     fTrackHypothesys.AddAt(array,esdindex);
3072   }
3073   array->AddLast(track);
3074 }
3075 //------------------------------------------------------------------------
3076 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3077 {
3078   //-------------------------------------------------------------------
3079   // compress array of track hypothesys
3080   // keep only maxsize best hypothesys
3081   //-------------------------------------------------------------------
3082   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3083   if (! (fTrackHypothesys.At(esdindex)) ) return;
3084   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3085   Int_t entries = array->GetEntriesFast();
3086   //
3087   //- find preliminary besttrack as a reference
3088   Float_t minchi2=10000;
3089   Int_t maxn=0;
3090   AliITStrackMI * besttrack=0;
3091   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3092     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3093     if (!track) continue;
3094     Float_t chi2 = NormalizedChi2(track,0);
3095     //
3096     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3097     track->SetLabel(tpcLabel);
3098     CookdEdx(track);
3099     track->SetFakeRatio(1.);
3100     CookLabel(track,0.); //For comparison only
3101     //
3102     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3103     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3104       if (track->GetNumberOfClusters()<maxn) continue;
3105       maxn = track->GetNumberOfClusters();
3106       if (chi2<minchi2){
3107         minchi2=chi2;
3108         besttrack=track;
3109       }
3110     }
3111     else{
3112       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3113         delete array->RemoveAt(itrack);
3114       }
3115     }
3116   }
3117   if (!besttrack) return;
3118   //
3119   //
3120   //take errors of best track as a reference
3121   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3122   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3123   for (Int_t j=0;j<6;j++) {
3124     if (besttrack->GetClIndex(j)>0){
3125       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3126       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3127       ny[j]   = besttrack->GetNy(j);
3128       nz[j]   = besttrack->GetNz(j);
3129     }
3130   }
3131   //
3132   // calculate normalized chi2
3133   //
3134   Float_t * chi2        = new Float_t[entries];
3135   Int_t * index         = new Int_t[entries];
3136   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3137   for (Int_t itrack=0;itrack<entries;itrack++){
3138     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3139     if (track){
3140       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3141       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3142         chi2[itrack] = track->GetChi2MIP(0);
3143       else{
3144         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3145           delete array->RemoveAt(itrack);
3146         }
3147       }
3148     }
3149   }
3150   //
3151   TMath::Sort(entries,chi2,index,kFALSE);
3152   besttrack = (AliITStrackMI*)array->At(index[0]);
3153   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3154     for (Int_t j=0;j<6;j++){
3155       if (besttrack->GetClIndex(j)>0){
3156         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3157         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3158         ny[j]   = besttrack->GetNy(j);
3159         nz[j]   = besttrack->GetNz(j);
3160       }
3161     }
3162   }
3163   //
3164   // calculate one more time with updated normalized errors
3165   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3166   for (Int_t itrack=0;itrack<entries;itrack++){
3167     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3168     if (track){
3169       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3170       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3171         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3172       else
3173         {
3174           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3175             delete array->RemoveAt(itrack);
3176           }
3177         }
3178     }
3179   }
3180   entries = array->GetEntriesFast();
3181   //
3182   //
3183   if (entries>0){
3184     TObjArray * newarray = new TObjArray();
3185     TMath::Sort(entries,chi2,index,kFALSE);
3186     besttrack = (AliITStrackMI*)array->At(index[0]);
3187     if (besttrack){
3188       //
3189       for (Int_t j=0;j<6;j++){
3190         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3191           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3192           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3193           ny[j]   = besttrack->GetNy(j);
3194           nz[j]   = besttrack->GetNz(j);
3195         }
3196       }
3197       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3198       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3199       Float_t minn = besttrack->GetNumberOfClusters()-3;
3200       Int_t accepted=0;
3201       for (Int_t i=0;i<entries;i++){
3202         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3203         if (!track) continue;
3204         if (accepted>maxcut) break;
3205         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3206         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3207           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3208             delete array->RemoveAt(index[i]);
3209             continue;
3210           }
3211         }
3212         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3213         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3214           if (!shortbest) accepted++;
3215           //
3216           newarray->AddLast(array->RemoveAt(index[i]));
3217           for (Int_t j=0;j<6;j++){
3218             if (nz[j]==0){
3219               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3220               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3221               ny[j]   = track->GetNy(j);
3222               nz[j]   = track->GetNz(j);
3223             }
3224           }
3225         }
3226         else{
3227           delete array->RemoveAt(index[i]);
3228         }
3229       }
3230       array->Delete();
3231       delete fTrackHypothesys.RemoveAt(esdindex);
3232       fTrackHypothesys.AddAt(newarray,esdindex);
3233     }
3234     else{
3235       array->Delete();
3236       delete fTrackHypothesys.RemoveAt(esdindex);
3237     }
3238   }
3239   delete [] chi2;
3240   delete [] index;
3241 }
3242 //------------------------------------------------------------------------
3243 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3244 {
3245   //-------------------------------------------------------------
3246   // try to find best hypothesy
3247   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3248   //-------------------------------------------------------------
3249   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3250   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3251   if (!array) return 0;
3252   Int_t entries = array->GetEntriesFast();
3253   if (!entries) return 0;
3254   Float_t minchi2 = 100000;
3255   AliITStrackMI * besttrack=0;
3256   //
3257   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3258   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3259   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3260   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3261   //
3262   for (Int_t i=0;i<entries;i++){
3263     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3264     if (!track) continue;
3265     Float_t sigmarfi,sigmaz;
3266     GetDCASigma(track,sigmarfi,sigmaz);
3267     track->SetDnorm(0,sigmarfi);
3268     track->SetDnorm(1,sigmaz);
3269     //
3270     track->SetChi2MIP(1,1000000);
3271     track->SetChi2MIP(2,1000000);
3272     track->SetChi2MIP(3,1000000);
3273     //
3274     // backtrack
3275     backtrack = new(backtrack) AliITStrackMI(*track);
3276     if (track->GetConstrain()) {
3277       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3278       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;
3279       backtrack->ResetCovariance(10.);
3280     }else{
3281       backtrack->ResetCovariance(10.);
3282     }
3283     backtrack->ResetClusters();
3284
3285     Double_t x = original->GetX();
3286     if (!RefitAt(x,backtrack,track)) continue;
3287     //
3288     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3289     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3290     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3291     track->SetChi22(GetMatchingChi2(backtrack,original));
3292
3293     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3294     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3295     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3296
3297
3298     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3299     //
3300     //forward track - without constraint
3301     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3302     forwardtrack->ResetClusters();
3303     x = track->GetX();
3304     RefitAt(x,forwardtrack,track);
3305     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3306     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3307     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3308
3309     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3310     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3311     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3312     forwardtrack->SetD(0,track->GetD(0));
3313     forwardtrack->SetD(1,track->GetD(1));
3314     {
3315       Int_t list[6];
3316       AliITSRecPoint* clist[6];
3317       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3318       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3319     }
3320
3321     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3322     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3323     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3324       track->SetChi2MIP(3,1000);
3325       continue;
3326     }
3327     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3328     //
3329     for (Int_t ichi=0;ichi<5;ichi++){
3330       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3331     }
3332     if (chi2 < minchi2){
3333       //besttrack = new AliITStrackMI(*forwardtrack);
3334       besttrack = track;
3335       besttrack->SetLabel(track->GetLabel());
3336       besttrack->SetFakeRatio(track->GetFakeRatio());
3337       minchi2   = chi2;
3338       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3339       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3340       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3341     }
3342   }
3343   delete backtrack;
3344   delete forwardtrack;
3345   Int_t accepted=0;
3346   for (Int_t i=0;i<entries;i++){
3347     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3348
3349     if (!track) continue;
3350
3351     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3352         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3353         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3354       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3355         delete array->RemoveAt(i);
3356         continue;
3357       }
3358     }
3359     else{
3360       accepted++;
3361     }
3362   }
3363   //
3364   array->Compress();
3365   SortTrackHypothesys(esdindex,checkmax,1);
3366
3367   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3368   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3369   besttrack = (AliITStrackMI*)array->At(0);
3370   if (!besttrack)  return 0;
3371   besttrack->SetChi2MIP(8,0);
3372   fBestTrackIndex[esdindex]=0;
3373   entries = array->GetEntriesFast();
3374   AliITStrackMI *longtrack =0;
3375   minchi2 =1000;
3376   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3377   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3378     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3379     if (!track->GetConstrain()) continue;
3380     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3381     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3382     if (track->GetChi2MIP(0)>4.) continue;
3383     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3384     longtrack =track;
3385   }
3386   //if (longtrack) besttrack=longtrack;
3387
3388   Int_t list[6];
3389   AliITSRecPoint * clist[6];
3390   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3391   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3392       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3393     RegisterClusterTracks(besttrack,esdindex);
3394   }
3395   //
3396   //
3397   if (shared>0.0){
3398     Int_t nshared;
3399     Int_t overlist[6];
3400     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3401     if (sharedtrack>=0){
3402       //
3403       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3404       if (besttrack){
3405         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3406       }
3407       else return 0;
3408     }
3409   }
3410
3411   if (shared>2.5) return 0;
3412   if (shared>1.0) return besttrack;
3413   //
3414   // Don't sign clusters if not gold track
3415   //
3416   if (!besttrack->IsGoldPrimary()) return besttrack;
3417   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3418   //
3419   if (fConstraint[fPass]){
3420     //
3421     // sign clusters
3422     //
3423     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3424     for (Int_t i=0;i<6;i++){
3425       Int_t index = besttrack->GetClIndex(i);
3426       if (index<=0) continue;
3427       Int_t ilayer =  (index & 0xf0000000) >> 28;
3428       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3429       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3430       if (!c) continue;
3431       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3432       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3433       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3434       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3435         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3436       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3437
3438       Bool_t cansign = kTRUE;
3439       for (Int_t itrack=0;itrack<entries; itrack++){
3440         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3441         if (!track) continue;
3442         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3443         if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3444           cansign = kFALSE;
3445           break;
3446         }
3447       }
3448       if (cansign){
3449         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3450         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3451         if (!c->IsUsed()) c->Use();
3452       }
3453     }
3454   }
3455   return besttrack;
3456 }
3457 //------------------------------------------------------------------------
3458 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3459 {
3460   //
3461   // get "best" hypothesys
3462   //
3463
3464   Int_t nentries = itsTracks.GetEntriesFast();
3465   for (Int_t i=0;i<nentries;i++){
3466     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3467     if (!track) continue;
3468     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3469     if (!array) continue;
3470     if (array->GetEntriesFast()<=0) continue;
3471     //
3472     AliITStrackMI* longtrack=0;
3473     Float_t minn=0;
3474     Float_t maxchi2=1000;
3475     for (Int_t j=0;j<array->GetEntriesFast();j++){
3476       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3477       if (!trackHyp) continue;
3478       if (trackHyp->GetGoldV0()) {
3479         longtrack = trackHyp;   //gold V0 track taken
3480         break;
3481       }
3482       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3483       Float_t chi2 = trackHyp->GetChi2MIP(0);
3484       if (fAfterV0){
3485         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3486       }
3487       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3488       //
3489       if (chi2 > maxchi2) continue;
3490       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3491       maxchi2 = chi2;
3492       longtrack=trackHyp;
3493     }
3494     //
3495     //
3496     //
3497     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3498     if (!longtrack) {longtrack = besttrack;}
3499     else besttrack= longtrack;
3500     //
3501     if (besttrack) {
3502       Int_t list[6];
3503       AliITSRecPoint * clist[6];
3504       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3505       //
3506       track->SetNUsed(shared);
3507       track->SetNSkipped(besttrack->GetNSkipped());
3508       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3509       if (shared>0) {
3510         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3511         Int_t nshared;
3512         Int_t overlist[6];
3513         //
3514         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3515         //if (sharedtrack==-1) sharedtrack=0;
3516         if (sharedtrack>=0) {
3517           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3518         }