ITS/AliITStrackerMI.cxx

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