ITS/AliITStrackerMI.cxx

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