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