ITS/AliITStrackerMI.cxx

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