ITS/AliITStrackerMI.cxx

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