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,AliTracker::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->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) 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->PropagateTo(xTrOrig+cl->GetX(),0.,0.)) 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->PropagateTo(xTrOrig,0.,0.)) 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
2215      Double_t phi,z;
2216      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2217
2218      Int_t idet=layer.FindDetectorIndex(phi,z);
2219
2220      // check if we allow a prolongation without point for large-eta tracks
2221      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2222      if (skip==2) {
2223        // propagate to the layer radius
2224        Double_t xToGo; if (!track->GetLocalXat(r,xToGo)) return kFALSE;
2225        if (!track->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return kFALSE;
2226        // apply correction for material of the current layer
2227        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2228        modstatus = 4; // out in z
2229        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2230          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2231        }
2232        // track time update [SR, GSI 17.02.2003]
2233        if (track->IsStartedTimeIntegral() && step==1) {
2234          Double_t newX, newY, newZ;
2235          if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2236          Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2237                         (oldZ-newZ)*(oldZ-newZ);
2238          track->AddTimeStep(TMath::Sqrt(dL2));
2239        }
2240        continue;
2241      }
2242
2243      if (idet<0) return kFALSE;
2244
2245      const AliITSdetector &det=layer.GetDetector(idet);
2246      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2247
2248      track->SetDetectorIndex(idet);
2249      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2250
2251      Double_t dz,zmin,zmax,dy,ymin,ymax;
2252
2253      const AliITSRecPoint *clAcc=0;
2254      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2255
2256      Int_t idx=index[ilayer];
2257      if (idx>=0) { // cluster in this layer
2258        modstatus = 6; // no refit
2259        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2260        if (cl) {
2261          if (idet != cl->GetDetectorIndex()) {
2262            idet=cl->GetDetectorIndex();
2263            const AliITSdetector &detc=layer.GetDetector(idet);
2264            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2265            track->SetDetectorIndex(idet);
2266            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2267          }
2268          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2269          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2270          if (chi2<maxchi2) {
2271            clAcc=cl;
2272            maxchi2=chi2;
2273            modstatus = 1; // found
2274          } else {
2275             return kFALSE; //
2276          }
2277        }
2278      } else { // no cluster in this layer
2279        if (skip==1) {
2280          modstatus = 3; // skipped
2281          // Plane Eff determination:
2282          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2283            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2284               UseTrackForPlaneEff(track,ilayer);
2285          }
2286        } else {
2287          modstatus = 5; // no cls in road
2288          // check dead
2289          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2290          dz = 0.5*(zmax-zmin);
2291          dy = 0.5*(ymax-ymin);
2292          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2293          if (dead==1) modstatus = 7; // holes in z in SPD
2294          if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2295        }
2296      }
2297
2298      if (clAcc) {
2299        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2300        track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2301      }
2302      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2303
2304
2305      if (extra) { // search for extra clusters in overlapped modules
2306        AliITStrackV2 tmp(*track);
2307        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2308        layer.SelectClusters(zmin,zmax,ymin,ymax);
2309
2310        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2311        Int_t idetExtra=-1;
2312        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2313        Double_t tolerance=0.1;
2314        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2315          // only clusters in another module! (overlaps)
2316          idetExtra = clExtra->GetDetectorIndex();
2317          if (idet == idetExtra) continue;
2318
2319          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2320
2321          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2322          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2323          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2324          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2325
2326          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2327          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2328        }
2329        if (cci>=0) {
2330          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2331          track->SetExtraModule(ilayer,idetExtra);
2332        }
2333      } // end search for extra clusters in overlapped modules
2334
2335      // Correct for material of the current layer
2336      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2337
2338      // track time update [SR, GSI 17.02.2003]
2339      if (track->IsStartedTimeIntegral() && step==1) {
2340         Double_t newX, newY, newZ;
2341         if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2342         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2343                        (oldZ-newZ)*(oldZ-newZ);
2344         track->AddTimeStep(TMath::Sqrt(dL2));
2345      }
2346      //
2347
2348   } // end loop on layers
2349
2350   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2351
2352   return kTRUE;
2353 }
2354 //------------------------------------------------------------------------
2355 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2356 {
2357   //
2358   // calculate normalized chi2
2359   //  return NormalizedChi2(track,0);
2360   Float_t chi2 = 0;
2361   Float_t sum=0;
2362   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2363   //  track->fdEdxMismatch=0;
2364   Float_t dedxmismatch =0;
2365   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2366   if (mode<100){
2367     for (Int_t i = 0;i<6;i++){
2368       if (track->GetClIndex(i)>0){
2369         Float_t cerry, cerrz;
2370         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2371         else
2372           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2373         cerry*=cerry;
2374         cerrz*=cerrz;
2375         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2376         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2377           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2378           if (ratio<0.5) {
2379             cchi2+=(0.5-ratio)*10.;
2380             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2381             dedxmismatch+=(0.5-ratio)*10.;
2382           }
2383         }
2384         if (i<2 ||i>3){
2385           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2386           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2387           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2388           if (i<2) chi2+=2*cl->GetDeltaProbability();
2389         }
2390         chi2+=cchi2;
2391         sum++;
2392       }
2393     }
2394     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2395       track->SetdEdxMismatch(dedxmismatch);
2396     }
2397   }
2398   else{
2399     for (Int_t i = 0;i<4;i++){
2400       if (track->GetClIndex(i)>0){
2401         Float_t cerry, cerrz;
2402         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2403         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2404         cerry*=cerry;
2405         cerrz*=cerrz;
2406         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2407         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2408         sum++;
2409       }
2410     }
2411     for (Int_t i = 4;i<6;i++){
2412       if (track->GetClIndex(i)>0){
2413         Float_t cerry, cerrz;
2414         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2415         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2416         cerry*=cerry;
2417         cerrz*=cerrz;
2418         Float_t cerryb, cerrzb;
2419         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2420         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2421         cerryb*=cerryb;
2422         cerrzb*=cerrzb;
2423         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2424         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2425         sum++;
2426       }
2427     }
2428   }
2429   if (track->GetESDtrack()->GetTPCsignal()>85){
2430     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2431     if (ratio<0.5) {
2432       chi2+=(0.5-ratio)*5.;
2433     }
2434     if (ratio>2){
2435       chi2+=(ratio-2.0)*3;
2436     }
2437   }
2438   //
2439   Double_t match = TMath::Sqrt(track->GetChi22());
2440   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2441   if (!track->GetConstrain()) {
2442     if (track->GetNumberOfClusters()>2) {
2443       match/=track->GetNumberOfClusters()-2.;
2444     } else {
2445       match=0;
2446     }
2447   }
2448   if (match<0) match=0;
2449   Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2450   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2451     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2452                                 1./(1.+track->GetNSkipped()));
2453
2454  return normchi2;
2455 }
2456 //------------------------------------------------------------------------
2457 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2458 {
2459   //
2460   // return matching chi2 between two tracks
2461   Double_t largeChi2=1000.;
2462
2463   AliITStrackMI track3(*track2);
2464   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2465   TMatrixD vec(5,1);
2466   vec(0,0)=track1->GetY()   - track3.GetY();
2467   vec(1,0)=track1->GetZ()   - track3.GetZ();
2468   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2469   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2470   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2471   //
2472   TMatrixD cov(5,5);
2473   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2474   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2475   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2476   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2477   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2478
2479   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2480   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2481   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2482   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2483   //
2484   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2485   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2486   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2487   //
2488   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2489   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2490   //
2491   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2492
2493   cov.Invert();
2494   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2495   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2496   return chi2(0,0);
2497 }
2498 //------------------------------------------------------------------------
2499 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr)
2500 {
2501   //
2502   //  return probability that given point (characterized by z position and error)
2503   //  is in SPD dead zone
2504   //
2505   Double_t probability = 0.;
2506   Double_t absz = TMath::Abs(zpos);
2507   Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2508                                   0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2509   if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2510   Double_t zmin, zmax;
2511   if (zpos<-6.) { // dead zone at z = -7
2512     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2513     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2514   } else if (zpos>6.) { // dead zone at z = +7
2515     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2516     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2517   } else if (absz<2.) { // dead zone at z = 0
2518     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2519     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2520   } else {
2521     zmin = 0.;
2522     zmax = 0.;
2523   }
2524   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2525   // dead zone)
2526   probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2527                       TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2528   return probability;
2529 }
2530 //------------------------------------------------------------------------
2531 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2532 {
2533   //
2534   // calculate normalized chi2
2535   Float_t chi2[6];
2536   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2537   Float_t ncl = 0;
2538   for (Int_t i = 0;i<6;i++){
2539     if (TMath::Abs(track->GetDy(i))>0){
2540       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2541       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2542       ncl++;
2543     }
2544     else{chi2[i]=10000;}
2545   }
2546   Int_t index[6];
2547   TMath::Sort(6,chi2,index,kFALSE);
2548   Float_t max = float(ncl)*fac-1.;
2549   Float_t sumchi=0, sumweight=0;
2550   for (Int_t i=0;i<max+1;i++){
2551     Float_t weight = (i<max)?1.:(max+1.-i);
2552     sumchi+=weight*chi2[index[i]];
2553     sumweight+=weight;
2554   }
2555   Double_t normchi2 = sumchi/sumweight;
2556   return normchi2;
2557 }
2558 //------------------------------------------------------------------------
2559 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2560 {
2561   //
2562   // calculate normalized chi2
2563   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2564   Int_t npoints = 0;
2565   Double_t res =0;
2566   for (Int_t i=0;i<6;i++){
2567     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2568     Double_t sy1 = forwardtrack->GetSigmaY(i);
2569     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2570     Double_t sy2 = backtrack->GetSigmaY(i);
2571     Double_t sz2 = backtrack->GetSigmaZ(i);
2572     if (i<2){ sy2=1000.;sz2=1000;}
2573     //
2574     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2575     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2576     //
2577     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2578     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2579     //
2580     res+= nz0*nz0+ny0*ny0;
2581     npoints++;
2582   }
2583   if (npoints>1) return
2584                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2585                    //2*forwardtrack->fNUsed+
2586                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2587                                   1./(1.+forwardtrack->GetNSkipped()));
2588   return 1000;
2589 }
2590 //------------------------------------------------------------------------
2591 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2592   //--------------------------------------------------------------------
2593   //       Return pointer to a given cluster
2594   //--------------------------------------------------------------------
2595   Int_t l=(index & 0xf0000000) >> 28;
2596   Int_t c=(index & 0x0fffffff) >> 00;
2597   return fgLayers[l].GetWeight(c);
2598 }
2599 //------------------------------------------------------------------------
2600 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2601 {
2602   //---------------------------------------------
2603   // register track to the list
2604   //
2605   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2606   //
2607   //
2608   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2609     Int_t index = track->GetClusterIndex(icluster);
2610     Int_t l=(index & 0xf0000000) >> 28;
2611     Int_t c=(index & 0x0fffffff) >> 00;
2612     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2613     for (Int_t itrack=0;itrack<4;itrack++){
2614       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2615         fgLayers[l].SetClusterTracks(itrack,c,id);
2616         break;
2617       }
2618     }
2619   }
2620 }
2621 //------------------------------------------------------------------------
2622 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2623 {
2624   //---------------------------------------------
2625   // unregister track from the list
2626   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2627     Int_t index = track->GetClusterIndex(icluster);
2628     Int_t l=(index & 0xf0000000) >> 28;
2629     Int_t c=(index & 0x0fffffff) >> 00;
2630     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2631     for (Int_t itrack=0;itrack<4;itrack++){
2632       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2633         fgLayers[l].SetClusterTracks(itrack,c,-1);
2634       }
2635     }
2636   }
2637 }
2638 //------------------------------------------------------------------------
2639 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2640 {
2641   //-------------------------------------------------------------
2642   //get number of shared clusters
2643   //-------------------------------------------------------------
2644   Float_t shared=0;
2645   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2646   // mean  number of clusters
2647   Float_t *ny = GetNy(id), *nz = GetNz(id);
2648
2649
2650   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2651     Int_t index = track->GetClusterIndex(icluster);
2652     Int_t l=(index & 0xf0000000) >> 28;
2653     Int_t c=(index & 0x0fffffff) >> 00;
2654     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2655     if (ny[l]==0){
2656       printf("problem\n");
2657     }
2658     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2659     Float_t weight=1;
2660     //
2661     Float_t deltan = 0;
2662     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2663     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2664       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2665     if (l<2 || l>3){
2666       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2667     }
2668     else{
2669       deltan = (cl->GetNz()-nz[l]);
2670     }
2671     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2672     weight = 2./TMath::Max(3.+deltan,2.);
2673     //
2674     for (Int_t itrack=0;itrack<4;itrack++){
2675       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2676         list[l]=index;
2677         clist[l] = (AliITSRecPoint*)GetCluster(index);
2678         shared+=weight;
2679         break;
2680       }
2681     }
2682   }
2683   track->SetNUsed(shared);
2684   return shared;
2685 }
2686 //------------------------------------------------------------------------
2687 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2688 {
2689   //
2690   // find first shared track
2691   //
2692   // mean  number of clusters
2693   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2694   //
2695   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2696   Int_t sharedtrack=100000;
2697   Int_t tracks[24],trackindex=0;
2698   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2699   //
2700   for (Int_t icluster=0;icluster<6;icluster++){
2701     if (clusterlist[icluster]<0) continue;
2702     Int_t index = clusterlist[icluster];
2703     Int_t l=(index & 0xf0000000) >> 28;
2704     Int_t c=(index & 0x0fffffff) >> 00;
2705     if (ny[l]==0){
2706       printf("problem\n");
2707     }
2708     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2709     //if (l>3) continue;
2710     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2711     //
2712     Float_t deltan = 0;
2713     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2714     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2715       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2716     if (l<2 || l>3){
2717       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2718     }
2719     else{
2720       deltan = (cl->GetNz()-nz[l]);
2721     }
2722     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2723     //
2724     for (Int_t itrack=3;itrack>=0;itrack--){
2725       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2726       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2727        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2728        trackindex++;
2729       }
2730     }
2731   }
2732   if (trackindex==0) return -1;
2733   if (trackindex==1){
2734     sharedtrack = tracks[0];
2735   }else{
2736     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2737     else{
2738       //
2739       Int_t tracks2[24], cluster[24];
2740       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2741       Int_t index =0;
2742       //
2743       for (Int_t i=0;i<trackindex;i++){
2744         if (tracks[i]<0) continue;
2745         tracks2[index] = tracks[i];
2746         cluster[index]++;
2747         for (Int_t j=i+1;j<trackindex;j++){
2748           if (tracks[j]<0) continue;
2749           if (tracks[j]==tracks[i]){
2750             cluster[index]++;
2751             tracks[j]=-1;
2752           }
2753         }
2754         index++;
2755       }
2756       Int_t max=0;
2757       for (Int_t i=0;i<index;i++){
2758         if (cluster[index]>max) {
2759           sharedtrack=tracks2[index];
2760           max=cluster[index];
2761         }
2762       }
2763     }
2764   }
2765
2766   if (sharedtrack>=100000) return -1;
2767   //
2768   // make list of overlaps
2769   shared =0;
2770   for (Int_t icluster=0;icluster<6;icluster++){
2771     if (clusterlist[icluster]<0) continue;
2772     Int_t index = clusterlist[icluster];
2773     Int_t l=(index & 0xf0000000) >> 28;
2774     Int_t c=(index & 0x0fffffff) >> 00;
2775     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2776     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2777     if (l==0 || l==1){
2778       if (cl->GetNy()>2) continue;
2779       if (cl->GetNz()>2) continue;
2780     }
2781     if (l==4 || l==5){
2782       if (cl->GetNy()>3) continue;
2783       if (cl->GetNz()>3) continue;
2784     }
2785     //
2786     for (Int_t itrack=3;itrack>=0;itrack--){
2787       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2788       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2789         overlist[l]=index;
2790         shared++;
2791       }
2792     }
2793   }
2794   return sharedtrack;
2795 }
2796 //------------------------------------------------------------------------
2797 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2798   //
2799   // try to find track hypothesys without conflicts
2800   // with minimal chi2;
2801   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2802   Int_t entries1 = arr1->GetEntriesFast();
2803   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2804   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2805   Int_t entries2 = arr2->GetEntriesFast();
2806   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2807   //
2808   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2809   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2810   if (track10->Pt()>0.5+track20->Pt()) return track10;
2811
2812   for (Int_t itrack=0;itrack<entries1;itrack++){
2813     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2814     UnRegisterClusterTracks(track,trackID1);
2815   }
2816   //
2817   for (Int_t itrack=0;itrack<entries2;itrack++){
2818     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2819     UnRegisterClusterTracks(track,trackID2);
2820   }
2821   Int_t index1=0;
2822   Int_t index2=0;
2823   Float_t maxconflicts=6;
2824   Double_t maxchi2 =1000.;
2825   //
2826   // get weight of hypothesys - using best hypothesy
2827   Double_t w1,w2;
2828
2829   Int_t list1[6],list2[6];
2830   AliITSRecPoint *clist1[6], *clist2[6] ;
2831   RegisterClusterTracks(track10,trackID1);
2832   RegisterClusterTracks(track20,trackID2);
2833   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2834   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2835   UnRegisterClusterTracks(track10,trackID1);
2836   UnRegisterClusterTracks(track20,trackID2);
2837   //
2838   // normalized chi2
2839   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2840   Float_t nerry[6],nerrz[6];
2841   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2842   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2843   for (Int_t i=0;i<6;i++){
2844      if ( (erry1[i]>0) && (erry2[i]>0)) {
2845        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2846        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2847      }else{
2848        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2849        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2850      }
2851      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2852        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2853        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2854        ncl1++;
2855      }
2856      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2857        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2858        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2859        ncl2++;
2860      }
2861   }
2862   chi21/=ncl1;
2863   chi22/=ncl2;
2864   //
2865   //
2866   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2867   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2868   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2869   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2870   //
2871   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2872         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2873         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2874         );
2875   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2876         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2877         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2878         );
2879
2880   Double_t sumw = w1+w2;
2881   w1/=sumw;
2882   w2/=sumw;
2883   if (w1<w2*0.5) {
2884     w1 = (d2+0.5)/(d1+d2+1.);
2885     w2 = (d1+0.5)/(d1+d2+1.);
2886   }
2887   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2888   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2889   //
2890   // get pair of "best" hypothesys
2891   //
2892   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2893   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2894
2895   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2896     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2897     //if (track1->fFakeRatio>0) continue;
2898     RegisterClusterTracks(track1,trackID1);
2899     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2900       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2901
2902       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2903       //if (track2->fFakeRatio>0) continue;
2904       Float_t nskipped=0;
2905       RegisterClusterTracks(track2,trackID2);
2906       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2907       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2908       UnRegisterClusterTracks(track2,trackID2);
2909       //
2910       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2911       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2912       if (nskipped>0.5) continue;
2913       //
2914       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2915       if (conflict1+1<cconflict1) continue;
2916       if (conflict2+1<cconflict2) continue;
2917       Float_t conflict=0;
2918       Float_t sumchi2=0;
2919       Float_t sum=0;
2920       for (Int_t i=0;i<6;i++){
2921         //
2922         Float_t c1 =0.; // conflict coeficients
2923         Float_t c2 =0.;
2924         if (clist1[i]&&clist2[i]){
2925           Float_t deltan = 0;
2926           if (i<2 || i>3){
2927             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2928           }
2929           else{
2930             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2931           }
2932           c1  = 2./TMath::Max(3.+deltan,2.);
2933           c2  = 2./TMath::Max(3.+deltan,2.);
2934         }
2935         else{
2936           if (clist1[i]){
2937             Float_t deltan = 0;
2938             if (i<2 || i>3){
2939               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2940             }
2941             else{
2942               deltan = (clist1[i]->GetNz()-nz1[i]);
2943             }
2944             c1  = 2./TMath::Max(3.+deltan,2.);
2945             c2  = 0;
2946           }
2947
2948           if (clist2[i]){
2949             Float_t deltan = 0;
2950             if (i<2 || i>3){
2951               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2952             }
2953             else{
2954               deltan = (clist2[i]->GetNz()-nz2[i]);
2955             }
2956             c2  = 2./TMath::Max(3.+deltan,2.);
2957             c1  = 0;
2958           }
2959         }
2960         //
2961         chi21=0;chi22=0;
2962         if (TMath::Abs(track1->GetDy(i))>0.) {
2963           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2964             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2965           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2966           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2967         }else{
2968           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2969         }
2970         //
2971         if (TMath::Abs(track2->GetDy(i))>0.) {
2972           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2973             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2974           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2975           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2976         }
2977         else{
2978           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2979         }
2980         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2981         if (chi21>0) sum+=w1;
2982         if (chi22>0) sum+=w2;
2983         conflict+=(c1+c2);
2984       }
2985       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2986       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2987       Double_t normchi2 = 2*conflict+sumchi2/sum;
2988       if ( normchi2 <maxchi2 ){
2989         index1 = itrack1;
2990         index2 = itrack2;
2991         maxconflicts = conflict;
2992         maxchi2 = normchi2;
2993       }
2994     }
2995     UnRegisterClusterTracks(track1,trackID1);
2996   }
2997   //
2998   //  if (maxconflicts<4 && maxchi2<th0){
2999   if (maxchi2<th0*2.){
3000     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3001     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3002     track1->SetChi2MIP(5,maxconflicts);
3003     track1->SetChi2MIP(6,maxchi2);
3004     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3005     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3006     track1->SetChi2MIP(8,index1);
3007     fBestTrackIndex[trackID1] =index1;
3008     UpdateESDtrack(track1, AliESDtrack::kITSin);
3009   }
3010   else if (track10->GetChi2MIP(0)<th1){
3011     track10->SetChi2MIP(5,maxconflicts);
3012     track10->SetChi2MIP(6,maxchi2);
3013     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3014     UpdateESDtrack(track10,AliESDtrack::kITSin);
3015   }
3016
3017   for (Int_t itrack=0;itrack<entries1;itrack++){
3018     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3019     UnRegisterClusterTracks(track,trackID1);
3020   }
3021   //
3022   for (Int_t itrack=0;itrack<entries2;itrack++){
3023     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3024     UnRegisterClusterTracks(track,trackID2);
3025   }
3026
3027   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3028       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3029     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3030   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3031     RegisterClusterTracks(track10,trackID1);
3032   }
3033   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3034       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3035     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3036     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3037     RegisterClusterTracks(track20,trackID2);
3038   }
3039   return track10;
3040
3041 }
3042 //------------------------------------------------------------------------
3043 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3044   //--------------------------------------------------------------------
3045   // This function marks clusters assigned to the track
3046   //--------------------------------------------------------------------
3047   AliTracker::UseClusters(t,from);
3048
3049   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3050   //if (c->GetQ()>2) c->Use();
3051   if (c->GetSigmaZ2()>0.1) c->Use();
3052   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3053   //if (c->GetQ()>2) c->Use();
3054   if (c->GetSigmaZ2()>0.1) c->Use();
3055
3056 }
3057 //------------------------------------------------------------------------
3058 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3059 {
3060   //------------------------------------------------------------------
3061   // add track to the list of hypothesys
3062   //------------------------------------------------------------------
3063
3064   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3065     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3066   //
3067   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3068   if (!array) {
3069     array = new TObjArray(10);
3070     fTrackHypothesys.AddAt(array,esdindex);
3071   }
3072   array->AddLast(track);
3073 }
3074 //------------------------------------------------------------------------
3075 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3076 {
3077   //-------------------------------------------------------------------
3078   // compress array of track hypothesys
3079   // keep only maxsize best hypothesys
3080   //-------------------------------------------------------------------
3081   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3082   if (! (fTrackHypothesys.At(esdindex)) ) return;
3083   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3084   Int_t entries = array->GetEntriesFast();
3085   //
3086   //- find preliminary besttrack as a reference
3087   Float_t minchi2=10000;
3088   Int_t maxn=0;
3089   AliITStrackMI * besttrack=0;
3090   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3091     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3092     if (!track) continue;
3093     Float_t chi2 = NormalizedChi2(track,0);
3094     //
3095     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3096     track->SetLabel(tpcLabel);
3097     CookdEdx(track);
3098     track->SetFakeRatio(1.);
3099     CookLabel(track,0.); //For comparison only
3100     //
3101     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3102     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3103       if (track->GetNumberOfClusters()<maxn) continue;
3104       maxn = track->GetNumberOfClusters();
3105       if (chi2<minchi2){
3106         minchi2=chi2;
3107         besttrack=track;
3108       }
3109     }
3110     else{
3111       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3112         delete array->RemoveAt(itrack);
3113       }
3114     }
3115   }
3116   if (!besttrack) return;
3117   //
3118   //
3119   //take errors of best track as a reference
3120   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3121   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3122   for (Int_t j=0;j<6;j++) {
3123     if (besttrack->GetClIndex(j)>0){
3124       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3125       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3126       ny[j]   = besttrack->GetNy(j);
3127       nz[j]   = besttrack->GetNz(j);
3128     }
3129   }
3130   //
3131   // calculate normalized chi2
3132   //
3133   Float_t * chi2        = new Float_t[entries];
3134   Int_t * index         = new Int_t[entries];
3135   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3136   for (Int_t itrack=0;itrack<entries;itrack++){
3137     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3138     if (track){
3139       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3140       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3141         chi2[itrack] = track->GetChi2MIP(0);
3142       else{
3143         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3144           delete array->RemoveAt(itrack);
3145         }
3146       }
3147     }
3148   }
3149   //
3150   TMath::Sort(entries,chi2,index,kFALSE);
3151   besttrack = (AliITStrackMI*)array->At(index[0]);
3152   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3153     for (Int_t j=0;j<6;j++){
3154       if (besttrack->GetClIndex(j)>0){
3155         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3156         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3157         ny[j]   = besttrack->GetNy(j);
3158         nz[j]   = besttrack->GetNz(j);
3159       }
3160     }
3161   }
3162   //
3163   // calculate one more time with updated normalized errors
3164   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3165   for (Int_t itrack=0;itrack<entries;itrack++){
3166     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3167     if (track){
3168       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3169       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3170         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3171       else
3172         {
3173           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3174             delete array->RemoveAt(itrack);
3175           }
3176         }
3177     }
3178   }
3179   entries = array->GetEntriesFast();
3180   //
3181   //
3182   if (entries>0){
3183     TObjArray * newarray = new TObjArray();
3184     TMath::Sort(entries,chi2,index,kFALSE);
3185     besttrack = (AliITStrackMI*)array->At(index[0]);
3186     if (besttrack){
3187       //
3188       for (Int_t j=0;j<6;j++){
3189         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3190           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3191           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3192           ny[j]   = besttrack->GetNy(j);
3193           nz[j]   = besttrack->GetNz(j);
3194         }
3195       }
3196       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3197       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3198       Float_t minn = besttrack->GetNumberOfClusters()-3;
3199       Int_t accepted=0;
3200       for (Int_t i=0;i<entries;i++){
3201         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3202         if (!track) continue;
3203         if (accepted>maxcut) break;
3204         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3205         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3206           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3207             delete array->RemoveAt(index[i]);
3208             continue;
3209           }
3210         }
3211         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3212         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3213           if (!shortbest) accepted++;
3214           //
3215           newarray->AddLast(array->RemoveAt(index[i]));
3216           for (Int_t j=0;j<6;j++){
3217             if (nz[j]==0){
3218               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3219               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3220               ny[j]   = track->GetNy(j);
3221               nz[j]   = track->GetNz(j);
3222             }
3223           }
3224         }
3225         else{
3226           delete array->RemoveAt(index[i]);
3227         }
3228       }
3229       array->Delete();
3230       delete fTrackHypothesys.RemoveAt(esdindex);
3231       fTrackHypothesys.AddAt(newarray,esdindex);
3232     }
3233     else{
3234       array->Delete();
3235       delete fTrackHypothesys.RemoveAt(esdindex);
3236     }
3237   }
3238   delete [] chi2;
3239   delete [] index;
3240 }
3241 //------------------------------------------------------------------------
3242 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3243 {
3244   //-------------------------------------------------------------
3245   // try to find best hypothesy
3246   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3247   //-------------------------------------------------------------
3248   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3249   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3250   if (!array) return 0;
3251   Int_t entries = array->GetEntriesFast();
3252   if (!entries) return 0;
3253   Float_t minchi2 = 100000;
3254   AliITStrackMI * besttrack=0;
3255   //
3256   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3257   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3258   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3259   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3260   //
3261   for (Int_t i=0;i<entries;i++){
3262     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3263     if (!track) continue;
3264     Float_t sigmarfi,sigmaz;
3265     GetDCASigma(track,sigmarfi,sigmaz);
3266     track->SetDnorm(0,sigmarfi);
3267     track->SetDnorm(1,sigmaz);
3268     //
3269     track->SetChi2MIP(1,1000000);
3270     track->SetChi2MIP(2,1000000);
3271     track->SetChi2MIP(3,1000000);
3272     //
3273     // backtrack
3274     backtrack = new(backtrack) AliITStrackMI(*track);
3275     if (track->GetConstrain()) {
3276       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3277       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;
3278       backtrack->ResetCovariance(10.);
3279     }else{
3280       backtrack->ResetCovariance(10.);
3281     }
3282     backtrack->ResetClusters();
3283
3284     Double_t x = original->GetX();
3285     if (!RefitAt(x,backtrack,track)) continue;
3286     //
3287     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3288     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3289     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3290     track->SetChi22(GetMatchingChi2(backtrack,original));
3291
3292     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3293     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3294     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3295
3296
3297     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3298     //
3299     //forward track - without constraint
3300     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3301     forwardtrack->ResetClusters();
3302     x = track->GetX();
3303     RefitAt(x,forwardtrack,track);
3304     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3305     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3306     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3307
3308     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3309     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3310     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3311     forwardtrack->SetD(0,track->GetD(0));
3312     forwardtrack->SetD(1,track->GetD(1));
3313     {
3314       Int_t list[6];
3315       AliITSRecPoint* clist[6];
3316       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3317       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3318     }
3319
3320     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3321     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3322     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3323       track->SetChi2MIP(3,1000);
3324       continue;
3325     }
3326     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3327     //
3328     for (Int_t ichi=0;ichi<5;ichi++){
3329       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3330     }
3331     if (chi2 < minchi2){
3332       //besttrack = new AliITStrackMI(*forwardtrack);
3333       besttrack = track;
3334       besttrack->SetLabel(track->GetLabel());
3335       besttrack->SetFakeRatio(track->GetFakeRatio());
3336       minchi2   = chi2;
3337       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3338       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3339       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3340     }
3341   }
3342   delete backtrack;
3343   delete forwardtrack;
3344   Int_t accepted=0;
3345   for (Int_t i=0;i<entries;i++){
3346     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3347
3348     if (!track) continue;
3349
3350     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3351         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3352         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3353       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3354         delete array->RemoveAt(i);
3355         continue;
3356       }
3357     }
3358     else{
3359       accepted++;
3360     }
3361   }
3362   //
3363   array->Compress();
3364   SortTrackHypothesys(esdindex,checkmax,1);
3365
3366   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3367   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3368   besttrack = (AliITStrackMI*)array->At(0);
3369   if (!besttrack)  return 0;
3370   besttrack->SetChi2MIP(8,0);
3371   fBestTrackIndex[esdindex]=0;
3372   entries = array->GetEntriesFast();
3373   AliITStrackMI *longtrack =0;
3374   minchi2 =1000;
3375   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3376   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3377     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3378     if (!track->GetConstrain()) continue;
3379     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3380     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3381     if (track->GetChi2MIP(0)>4.) continue;
3382     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3383     longtrack =track;
3384   }
3385   //if (longtrack) besttrack=longtrack;
3386
3387   Int_t list[6];
3388   AliITSRecPoint * clist[6];
3389   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3390   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3391       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3392     RegisterClusterTracks(besttrack,esdindex);
3393   }
3394   //
3395   //
3396   if (shared>0.0){
3397     Int_t nshared;
3398     Int_t overlist[6];
3399     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3400     if (sharedtrack>=0){
3401       //
3402       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3403       if (besttrack){
3404         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3405       }
3406       else return 0;
3407     }
3408   }
3409
3410   if (shared>2.5) return 0;
3411   if (shared>1.0) return besttrack;
3412   //
3413   // Don't sign clusters if not gold track
3414   //
3415   if (!besttrack->IsGoldPrimary()) return besttrack;
3416   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3417   //
3418   if (fConstraint[fPass]){
3419     //
3420     // sign clusters
3421     //
3422     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3423     for (Int_t i=0;i<6;i++){
3424       Int_t index = besttrack->GetClIndex(i);
3425       if (index<=0) continue;
3426       Int_t ilayer =  (index & 0xf0000000) >> 28;
3427       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3428       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3429       if (!c) continue;
3430       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3431       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3432       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3433       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3434         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3435       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3436
3437       Bool_t cansign = kTRUE;
3438       for (Int_t itrack=0;itrack<entries; itrack++){
3439         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3440         if (!track) continue;
3441         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3442         if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3443           cansign = kFALSE;
3444           break;
3445         }
3446       }
3447       if (cansign){
3448         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3449         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3450         if (!c->IsUsed()) c->Use();
3451       }
3452     }
3453   }
3454   return besttrack;
3455 }
3456 //------------------------------------------------------------------------
3457 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3458 {
3459   //
3460   // get "best" hypothesys
3461   //
3462
3463   Int_t nentries = itsTracks.GetEntriesFast();
3464   for (Int_t i=0;i<nentries;i++){
3465     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3466     if (!track) continue;
3467     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3468     if (!array) continue;
3469     if (array->GetEntriesFast()<=0) continue;
3470     //
3471     AliITStrackMI* longtrack=0;
3472     Float_t minn=0;
3473     Float_t maxchi2=1000;
3474     for (Int_t j=0;j<array->GetEntriesFast();j++){
3475       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3476       if (!trackHyp) continue;
3477       if (trackHyp->GetGoldV0()) {
3478         longtrack = trackHyp;   //gold V0 track taken
3479         break;
3480       }
3481       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3482       Float_t chi2 = trackHyp->GetChi2MIP(0);
3483       if (fAfterV0){
3484         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3485       }
3486       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3487       //
3488       if (chi2 > maxchi2) continue;
3489       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3490       maxchi2 = chi2;
3491       longtrack=trackHyp;
3492     }
3493     //
3494     //
3495     //
3496     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3497     if (!longtrack) {longtrack = besttrack;}
3498     else besttrack= longtrack;
3499     //
3500     if (besttrack) {
3501       Int_t list[6];
3502       AliITSRecPoint * clist[6];
3503       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3504       //
3505       track->SetNUsed(shared);
3506       track->SetNSkipped(besttrack->GetNSkipped());
3507       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3508       if (shared>0) {
3509         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3510         Int_t nshared;
3511         Int_t overlist[6];
3512         //
3513         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3514         //if (sharedtrack==-1) sharedtrack=0;
3515         if (sharedtrack>=0) {
3516           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3517         }
3518       }