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