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