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