ITS/AliITStrackerMI.cxx

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