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