ITS/AliITStrackerMI.cxx

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