ITS/AliITStrackerMI.cxx

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