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