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