ITS/AliITStrackerMI.cxx

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