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          Double_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,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->Propagate(xToGo)) 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->Propagate(xTrOrig+cl->GetX())) 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->Propagate(xTrOrig)) 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      //TMath::Sqrt(track->GetSigmaY2());
2215
2216      Double_t phi,z;
2217      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2218
2219      Int_t idet=layer.FindDetectorIndex(phi,z);
2220
2221      // check if we allow a prolongation without point for large-eta tracks
2222      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2223      if (skip==2) {
2224        // propagate to the layer radius
2225        Double_t xToGo; if (!track->GetLocalXat(r,xToGo)) return kFALSE;
2226        if (!track->Propagate(xToGo)) return kFALSE;
2227        // apply correction for material of the current layer
2228        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2229        modstatus = 4; // out in z
2230        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2231          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2232        }
2233        // track time update [SR, GSI 17.02.2003]
2234        if (track->IsStartedTimeIntegral() && step==1) {
2235          Double_t newX, newY, newZ;
2236          if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2237          Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2238                         (oldZ-newZ)*(oldZ-newZ);
2239          track->AddTimeStep(TMath::Sqrt(dL2));
2240        }
2241        continue;
2242      }
2243
2244      if (idet<0) return kFALSE;
2245
2246      const AliITSdetector &det=layer.GetDetector(idet);
2247      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2248
2249      track->SetDetectorIndex(idet);
2250      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2251
2252      Double_t dz,zmin,zmax,dy,ymin,ymax;
2253
2254      const AliITSRecPoint *clAcc=0;
2255      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2256
2257      Int_t idx=index[ilayer];
2258      if (idx>=0) { // cluster in this layer
2259        modstatus = 6; // no refit
2260        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2261        if (cl) {
2262          if (idet != cl->GetDetectorIndex()) {
2263            idet=cl->GetDetectorIndex();
2264            const AliITSdetector &detc=layer.GetDetector(idet);
2265            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2266            track->SetDetectorIndex(idet);
2267            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2268          }
2269          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2270          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2271          if (chi2<maxchi2) {
2272            clAcc=cl;
2273            maxchi2=chi2;
2274            modstatus = 1; // found
2275          } else {
2276             return kFALSE; //
2277          }
2278        }
2279      } else { // no cluster in this layer
2280        if (skip==1) {
2281          modstatus = 3; // skipped
2282          // Plane Eff determination:
2283          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2284            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2285               UseTrackForPlaneEff(track,ilayer);
2286          }
2287        } else {
2288          modstatus = 5; // no cls in road
2289          // check dead
2290          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2291          dz = 0.5*(zmax-zmin);
2292          dy = 0.5*(ymax-ymin);
2293          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2294          if (dead==1) modstatus = 7; // holes in z in SPD
2295          if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2296        }
2297      }
2298
2299      if (clAcc) {
2300        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2301        track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2302      }
2303      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2304
2305
2306      if (extra) { // search for extra clusters in overlapped modules
2307        AliITStrackV2 tmp(*track);
2308        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2309        layer.SelectClusters(zmin,zmax,ymin,ymax);
2310
2311        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2312        Int_t idetExtra=-1;
2313        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2314        Double_t tolerance=0.1;
2315        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2316          // only clusters in another module! (overlaps)
2317          idetExtra = clExtra->GetDetectorIndex();
2318          if (idet == idetExtra) continue;
2319
2320          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2321
2322          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2323          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2324          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2325          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2326
2327          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2328          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2329        }
2330        if (cci>=0) {
2331          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2332          track->SetExtraModule(ilayer,idetExtra);
2333        }
2334      } // end search for extra clusters in overlapped modules
2335
2336      // Correct for material of the current layer
2337      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2338
2339      // track time update [SR, GSI 17.02.2003]
2340      if (track->IsStartedTimeIntegral() && step==1) {
2341         Double_t newX, newY, newZ;
2342         if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2343         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2344                        (oldZ-newZ)*(oldZ-newZ);
2345         track->AddTimeStep(TMath::Sqrt(dL2));
2346      }
2347      //
2348
2349   } // end loop on layers
2350
2351   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2352
2353   return kTRUE;
2354 }
2355 //------------------------------------------------------------------------
2356 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2357 {
2358   //
2359   // calculate normalized chi2
2360   //  return NormalizedChi2(track,0);
2361   Float_t chi2 = 0;
2362   Float_t sum=0;
2363   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2364   //  track->fdEdxMismatch=0;
2365   Float_t dedxmismatch =0;
2366   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2367   if (mode<100){
2368     for (Int_t i = 0;i<6;i++){
2369       if (track->GetClIndex(i)>0){
2370         Float_t cerry, cerrz;
2371         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2372         else
2373           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2374         cerry*=cerry;
2375         cerrz*=cerrz;
2376         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2377         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2378           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2379           if (ratio<0.5) {
2380             cchi2+=(0.5-ratio)*10.;
2381             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2382             dedxmismatch+=(0.5-ratio)*10.;
2383           }
2384         }
2385         if (i<2 ||i>3){
2386           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2387           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2388           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2389           if (i<2) chi2+=2*cl->GetDeltaProbability();
2390         }
2391         chi2+=cchi2;
2392         sum++;
2393       }
2394     }
2395     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2396       track->SetdEdxMismatch(dedxmismatch);
2397     }
2398   }
2399   else{
2400     for (Int_t i = 0;i<4;i++){
2401       if (track->GetClIndex(i)>0){
2402         Float_t cerry, cerrz;
2403         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2404         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2405         cerry*=cerry;
2406         cerrz*=cerrz;
2407         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2408         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2409         sum++;
2410       }
2411     }
2412     for (Int_t i = 4;i<6;i++){
2413       if (track->GetClIndex(i)>0){
2414         Float_t cerry, cerrz;
2415         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2416         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2417         cerry*=cerry;
2418         cerrz*=cerrz;
2419         Float_t cerryb, cerrzb;
2420         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2421         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2422         cerryb*=cerryb;
2423         cerrzb*=cerrzb;
2424         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2425         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2426         sum++;
2427       }
2428     }
2429   }
2430   if (track->GetESDtrack()->GetTPCsignal()>85){
2431     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2432     if (ratio<0.5) {
2433       chi2+=(0.5-ratio)*5.;
2434     }
2435     if (ratio>2){
2436       chi2+=(ratio-2.0)*3;
2437     }
2438   }
2439   //
2440   Double_t match = TMath::Sqrt(track->GetChi22());
2441   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2442   if (!track->GetConstrain()) {
2443     if (track->GetNumberOfClusters()>2) {
2444       match/=track->GetNumberOfClusters()-2.;
2445     } else {
2446       match=0;
2447     }
2448   }
2449   if (match<0) match=0;
2450   Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2451   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2452     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2453                                 1./(1.+track->GetNSkipped()));
2454
2455  return normchi2;
2456 }
2457 //------------------------------------------------------------------------
2458 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2459 {
2460   //
2461   // return matching chi2 between two tracks
2462   Double_t largeChi2=1000.;
2463
2464   AliITStrackMI track3(*track2);
2465   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2466   TMatrixD vec(5,1);
2467   vec(0,0)=track1->GetY()   - track3.GetY();
2468   vec(1,0)=track1->GetZ()   - track3.GetZ();
2469   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2470   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2471   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2472   //
2473   TMatrixD cov(5,5);
2474   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2475   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2476   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2477   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2478   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2479
2480   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2481   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2482   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2483   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2484   //
2485   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2486   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2487   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2488   //
2489   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2490   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2491   //
2492   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2493
2494   cov.Invert();
2495   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2496   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2497   return chi2(0,0);
2498 }
2499 //------------------------------------------------------------------------
2500 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr)
2501 {
2502   //
2503   //  return probability that given point (characterized by z position and error)
2504   //  is in SPD dead zone
2505   //
2506   Double_t probability = 0.;
2507   Double_t absz = TMath::Abs(zpos);
2508   Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2509                                   0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2510   if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2511   Double_t zmin, zmax;
2512   if (zpos<-6.) { // dead zone at z = -7
2513     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2514     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2515   } else if (zpos>6.) { // dead zone at z = +7
2516     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2517     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2518   } else if (absz<2.) { // dead zone at z = 0
2519     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2520     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2521   } else {
2522     zmin = 0.;
2523     zmax = 0.;
2524   }
2525   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2526   // dead zone)
2527   probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2528                       TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2529   return probability;
2530 }
2531 //------------------------------------------------------------------------
2532 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2533 {
2534   //
2535   // calculate normalized chi2
2536   Float_t chi2[6];
2537   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2538   Float_t ncl = 0;
2539   for (Int_t i = 0;i<6;i++){
2540     if (TMath::Abs(track->GetDy(i))>0){
2541       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2542       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2543       ncl++;
2544     }
2545     else{chi2[i]=10000;}
2546   }
2547   Int_t index[6];
2548   TMath::Sort(6,chi2,index,kFALSE);
2549   Float_t max = float(ncl)*fac-1.;
2550   Float_t sumchi=0, sumweight=0;
2551   for (Int_t i=0;i<max+1;i++){
2552     Float_t weight = (i<max)?1.:(max+1.-i);
2553     sumchi+=weight*chi2[index[i]];
2554     sumweight+=weight;
2555   }
2556   Double_t normchi2 = sumchi/sumweight;
2557   return normchi2;
2558 }
2559 //------------------------------------------------------------------------
2560 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2561 {
2562   //
2563   // calculate normalized chi2
2564   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2565   Int_t npoints = 0;
2566   Double_t res =0;
2567   for (Int_t i=0;i<6;i++){
2568     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2569     Double_t sy1 = forwardtrack->GetSigmaY(i);
2570     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2571     Double_t sy2 = backtrack->GetSigmaY(i);
2572     Double_t sz2 = backtrack->GetSigmaZ(i);
2573     if (i<2){ sy2=1000.;sz2=1000;}
2574     //
2575     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2576     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2577     //
2578     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2579     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2580     //
2581     res+= nz0*nz0+ny0*ny0;
2582     npoints++;
2583   }
2584   if (npoints>1) return
2585                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2586                    //2*forwardtrack->fNUsed+
2587                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2588                                   1./(1.+forwardtrack->GetNSkipped()));
2589   return 1000;
2590 }
2591 //------------------------------------------------------------------------
2592 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2593   //--------------------------------------------------------------------
2594   //       Return pointer to a given cluster
2595   //--------------------------------------------------------------------
2596   Int_t l=(index & 0xf0000000) >> 28;
2597   Int_t c=(index & 0x0fffffff) >> 00;
2598   return fgLayers[l].GetWeight(c);
2599 }
2600 //------------------------------------------------------------------------
2601 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2602 {
2603   //---------------------------------------------
2604   // register track to the list
2605   //
2606   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2607   //
2608   //
2609   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2610     Int_t index = track->GetClusterIndex(icluster);
2611     Int_t l=(index & 0xf0000000) >> 28;
2612     Int_t c=(index & 0x0fffffff) >> 00;
2613     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2614     for (Int_t itrack=0;itrack<4;itrack++){
2615       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2616         fgLayers[l].SetClusterTracks(itrack,c,id);
2617         break;
2618       }
2619     }
2620   }
2621 }
2622 //------------------------------------------------------------------------
2623 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2624 {
2625   //---------------------------------------------
2626   // unregister track from the list
2627   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2628     Int_t index = track->GetClusterIndex(icluster);
2629     Int_t l=(index & 0xf0000000) >> 28;
2630     Int_t c=(index & 0x0fffffff) >> 00;
2631     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2632     for (Int_t itrack=0;itrack<4;itrack++){
2633       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2634         fgLayers[l].SetClusterTracks(itrack,c,-1);
2635       }
2636     }
2637   }
2638 }
2639 //------------------------------------------------------------------------
2640 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2641 {
2642   //-------------------------------------------------------------
2643   //get number of shared clusters
2644   //-------------------------------------------------------------
2645   Float_t shared=0;
2646   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2647   // mean  number of clusters
2648   Float_t *ny = GetNy(id), *nz = GetNz(id);
2649
2650
2651   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2652     Int_t index = track->GetClusterIndex(icluster);
2653     Int_t l=(index & 0xf0000000) >> 28;
2654     Int_t c=(index & 0x0fffffff) >> 00;
2655     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2656     if (ny[l]==0){
2657       printf("problem\n");
2658     }
2659     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2660     Float_t weight=1;
2661     //
2662     Float_t deltan = 0;
2663     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2664     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2665       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2666     if (l<2 || l>3){
2667       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2668     }
2669     else{
2670       deltan = (cl->GetNz()-nz[l]);
2671     }
2672     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2673     weight = 2./TMath::Max(3.+deltan,2.);
2674     //
2675     for (Int_t itrack=0;itrack<4;itrack++){
2676       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2677         list[l]=index;
2678         clist[l] = (AliITSRecPoint*)GetCluster(index);
2679         shared+=weight;
2680         break;
2681       }
2682     }
2683   }
2684   track->SetNUsed(shared);
2685   return shared;
2686 }
2687 //------------------------------------------------------------------------
2688 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2689 {
2690   //
2691   // find first shared track
2692   //
2693   // mean  number of clusters
2694   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2695   //
2696   for (Int_t i=0;i<6;i++) overlist[i]=-1;
2697   Int_t sharedtrack=100000;
2698   Int_t tracks[24],trackindex=0;
2699   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2700   //
2701   for (Int_t icluster=0;icluster<6;icluster++){
2702     if (clusterlist[icluster]<0) continue;
2703     Int_t index = clusterlist[icluster];
2704     Int_t l=(index & 0xf0000000) >> 28;
2705     Int_t c=(index & 0x0fffffff) >> 00;
2706     if (ny[l]==0){
2707       printf("problem\n");
2708     }
2709     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2710     //if (l>3) continue;
2711     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2712     //
2713     Float_t deltan = 0;
2714     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2715     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2716       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2717     if (l<2 || l>3){
2718       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2719     }
2720     else{
2721       deltan = (cl->GetNz()-nz[l]);
2722     }
2723     if (deltan>2.0) continue;  // extended - highly probable shared cluster
2724     //
2725     for (Int_t itrack=3;itrack>=0;itrack--){
2726       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2727       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2728        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
2729        trackindex++;
2730       }
2731     }
2732   }
2733   if (trackindex==0) return -1;
2734   if (trackindex==1){
2735     sharedtrack = tracks[0];
2736   }else{
2737     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2738     else{
2739       //
2740       Int_t tracks2[24], cluster[24];
2741       for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2742       Int_t index =0;
2743       //
2744       for (Int_t i=0;i<trackindex;i++){
2745         if (tracks[i]<0) continue;
2746         tracks2[index] = tracks[i];
2747         cluster[index]++;
2748         for (Int_t j=i+1;j<trackindex;j++){
2749           if (tracks[j]<0) continue;
2750           if (tracks[j]==tracks[i]){
2751             cluster[index]++;
2752             tracks[j]=-1;
2753           }
2754         }
2755         index++;
2756       }
2757       Int_t max=0;
2758       for (Int_t i=0;i<index;i++){
2759         if (cluster[index]>max) {
2760           sharedtrack=tracks2[index];
2761           max=cluster[index];
2762         }
2763       }
2764     }
2765   }
2766
2767   if (sharedtrack>=100000) return -1;
2768   //
2769   // make list of overlaps
2770   shared =0;
2771   for (Int_t icluster=0;icluster<6;icluster++){
2772     if (clusterlist[icluster]<0) continue;
2773     Int_t index = clusterlist[icluster];
2774     Int_t l=(index & 0xf0000000) >> 28;
2775     Int_t c=(index & 0x0fffffff) >> 00;
2776     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2777     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2778     if (l==0 || l==1){
2779       if (cl->GetNy()>2) continue;
2780       if (cl->GetNz()>2) continue;
2781     }
2782     if (l==4 || l==5){
2783       if (cl->GetNy()>3) continue;
2784       if (cl->GetNz()>3) continue;
2785     }
2786     //
2787     for (Int_t itrack=3;itrack>=0;itrack--){
2788       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2789       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2790         overlist[l]=index;
2791         shared++;
2792       }
2793     }
2794   }
2795   return sharedtrack;
2796 }
2797 //------------------------------------------------------------------------
2798 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2799   //
2800   // try to find track hypothesys without conflicts
2801   // with minimal chi2;
2802   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2803   Int_t entries1 = arr1->GetEntriesFast();
2804   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2805   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2806   Int_t entries2 = arr2->GetEntriesFast();
2807   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2808   //
2809   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2810   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2811   if (track10->Pt()>0.5+track20->Pt()) return track10;
2812
2813   for (Int_t itrack=0;itrack<entries1;itrack++){
2814     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2815     UnRegisterClusterTracks(track,trackID1);
2816   }
2817   //
2818   for (Int_t itrack=0;itrack<entries2;itrack++){
2819     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2820     UnRegisterClusterTracks(track,trackID2);
2821   }
2822   Int_t index1=0;
2823   Int_t index2=0;
2824   Float_t maxconflicts=6;
2825   Double_t maxchi2 =1000.;
2826   //
2827   // get weight of hypothesys - using best hypothesy
2828   Double_t w1,w2;
2829
2830   Int_t list1[6],list2[6];
2831   AliITSRecPoint *clist1[6], *clist2[6] ;
2832   RegisterClusterTracks(track10,trackID1);
2833   RegisterClusterTracks(track20,trackID2);
2834   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2835   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2836   UnRegisterClusterTracks(track10,trackID1);
2837   UnRegisterClusterTracks(track20,trackID2);
2838   //
2839   // normalized chi2
2840   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2841   Float_t nerry[6],nerrz[6];
2842   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2843   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2844   for (Int_t i=0;i<6;i++){
2845      if ( (erry1[i]>0) && (erry2[i]>0)) {
2846        nerry[i] = TMath::Min(erry1[i],erry2[i]);
2847        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2848      }else{
2849        nerry[i] = TMath::Max(erry1[i],erry2[i]);
2850        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2851      }
2852      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2853        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2854        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2855        ncl1++;
2856      }
2857      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2858        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2859        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2860        ncl2++;
2861      }
2862   }
2863   chi21/=ncl1;
2864   chi22/=ncl2;
2865   //
2866   //
2867   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2868   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2869   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2870   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2871   //
2872   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2873         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2874         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2875         );
2876   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2877         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2878         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2879         );
2880
2881   Double_t sumw = w1+w2;
2882   w1/=sumw;
2883   w2/=sumw;
2884   if (w1<w2*0.5) {
2885     w1 = (d2+0.5)/(d1+d2+1.);
2886     w2 = (d1+0.5)/(d1+d2+1.);
2887   }
2888   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2889   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2890   //
2891   // get pair of "best" hypothesys
2892   //
2893   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2894   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2895
2896   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2897     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2898     //if (track1->fFakeRatio>0) continue;
2899     RegisterClusterTracks(track1,trackID1);
2900     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2901       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2902
2903       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2904       //if (track2->fFakeRatio>0) continue;
2905       Float_t nskipped=0;
2906       RegisterClusterTracks(track2,trackID2);
2907       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2908       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2909       UnRegisterClusterTracks(track2,trackID2);
2910       //
2911       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2912       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2913       if (nskipped>0.5) continue;
2914       //
2915       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2916       if (conflict1+1<cconflict1) continue;
2917       if (conflict2+1<cconflict2) continue;
2918       Float_t conflict=0;
2919       Float_t sumchi2=0;
2920       Float_t sum=0;
2921       for (Int_t i=0;i<6;i++){
2922         //
2923         Float_t c1 =0.; // conflict coeficients
2924         Float_t c2 =0.;
2925         if (clist1[i]&&clist2[i]){
2926           Float_t deltan = 0;
2927           if (i<2 || i>3){
2928             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2929           }
2930           else{
2931             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2932           }
2933           c1  = 2./TMath::Max(3.+deltan,2.);
2934           c2  = 2./TMath::Max(3.+deltan,2.);
2935         }
2936         else{
2937           if (clist1[i]){
2938             Float_t deltan = 0;
2939             if (i<2 || i>3){
2940               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2941             }
2942             else{
2943               deltan = (clist1[i]->GetNz()-nz1[i]);
2944             }
2945             c1  = 2./TMath::Max(3.+deltan,2.);
2946             c2  = 0;
2947           }
2948
2949           if (clist2[i]){
2950             Float_t deltan = 0;
2951             if (i<2 || i>3){
2952               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2953             }
2954             else{
2955               deltan = (clist2[i]->GetNz()-nz2[i]);
2956             }
2957             c2  = 2./TMath::Max(3.+deltan,2.);
2958             c1  = 0;
2959           }
2960         }
2961         //
2962         chi21=0;chi22=0;
2963         if (TMath::Abs(track1->GetDy(i))>0.) {
2964           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2965             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2966           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2967           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2968         }else{
2969           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2970         }
2971         //
2972         if (TMath::Abs(track2->GetDy(i))>0.) {
2973           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2974             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2975           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2976           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2977         }
2978         else{
2979           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2980         }
2981         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2982         if (chi21>0) sum+=w1;
2983         if (chi22>0) sum+=w2;
2984         conflict+=(c1+c2);
2985       }
2986       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2987       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2988       Double_t normchi2 = 2*conflict+sumchi2/sum;
2989       if ( normchi2 <maxchi2 ){
2990         index1 = itrack1;
2991         index2 = itrack2;
2992         maxconflicts = conflict;
2993         maxchi2 = normchi2;
2994       }
2995     }
2996     UnRegisterClusterTracks(track1,trackID1);
2997   }
2998   //
2999   //  if (maxconflicts<4 && maxchi2<th0){
3000   if (maxchi2<th0*2.){
3001     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3002     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3003     track1->SetChi2MIP(5,maxconflicts);
3004     track1->SetChi2MIP(6,maxchi2);
3005     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3006     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3007     track1->SetChi2MIP(8,index1);
3008     fBestTrackIndex[trackID1] =index1;
3009     UpdateESDtrack(track1, AliESDtrack::kITSin);
3010   }
3011   else if (track10->GetChi2MIP(0)<th1){
3012     track10->SetChi2MIP(5,maxconflicts);
3013     track10->SetChi2MIP(6,maxchi2);
3014     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3015     UpdateESDtrack(track10,AliESDtrack::kITSin);
3016   }
3017
3018   for (Int_t itrack=0;itrack<entries1;itrack++){
3019     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3020     UnRegisterClusterTracks(track,trackID1);
3021   }
3022   //
3023   for (Int_t itrack=0;itrack<entries2;itrack++){
3024     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3025     UnRegisterClusterTracks(track,trackID2);
3026   }
3027
3028   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3029       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3030     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3031   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3032     RegisterClusterTracks(track10,trackID1);
3033   }
3034   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3035       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3036     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3037     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3038     RegisterClusterTracks(track20,trackID2);
3039   }
3040   return track10;
3041
3042 }
3043 //------------------------------------------------------------------------
3044 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3045   //--------------------------------------------------------------------
3046   // This function marks clusters assigned to the track
3047   //--------------------------------------------------------------------
3048   AliTracker::UseClusters(t,from);
3049
3050   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3051   //if (c->GetQ()>2) c->Use();
3052   if (c->GetSigmaZ2()>0.1) c->Use();
3053   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3054   //if (c->GetQ()>2) c->Use();
3055   if (c->GetSigmaZ2()>0.1) c->Use();
3056
3057 }
3058 //------------------------------------------------------------------------
3059 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3060 {
3061   //------------------------------------------------------------------
3062   // add track to the list of hypothesys
3063   //------------------------------------------------------------------
3064
3065   if (esdindex>=fTrackHypothesys.GetEntriesFast())
3066     fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3067   //
3068   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3069   if (!array) {
3070     array = new TObjArray(10);
3071     fTrackHypothesys.AddAt(array,esdindex);
3072   }
3073   array->AddLast(track);
3074 }
3075 //------------------------------------------------------------------------
3076 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3077 {
3078   //-------------------------------------------------------------------
3079   // compress array of track hypothesys
3080   // keep only maxsize best hypothesys
3081   //-------------------------------------------------------------------
3082   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3083   if (! (fTrackHypothesys.At(esdindex)) ) return;
3084   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3085   Int_t entries = array->GetEntriesFast();
3086   //
3087   //- find preliminary besttrack as a reference
3088   Float_t minchi2=10000;
3089   Int_t maxn=0;
3090   AliITStrackMI * besttrack=0;
3091   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3092     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3093     if (!track) continue;
3094     Float_t chi2 = NormalizedChi2(track,0);
3095     //
3096     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3097     track->SetLabel(tpcLabel);
3098     CookdEdx(track);
3099     track->SetFakeRatio(1.);
3100     CookLabel(track,0.); //For comparison only
3101     //
3102     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3103     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3104       if (track->GetNumberOfClusters()<maxn) continue;
3105       maxn = track->GetNumberOfClusters();
3106       if (chi2<minchi2){
3107         minchi2=chi2;
3108         besttrack=track;
3109       }
3110     }
3111     else{
3112       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3113         delete array->RemoveAt(itrack);
3114       }
3115     }
3116   }
3117   if (!besttrack) return;
3118   //
3119   //
3120   //take errors of best track as a reference
3121   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3122   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3123   for (Int_t j=0;j<6;j++) {
3124     if (besttrack->GetClIndex(j)>0){
3125       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3126       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3127       ny[j]   = besttrack->GetNy(j);
3128       nz[j]   = besttrack->GetNz(j);
3129     }
3130   }
3131   //
3132   // calculate normalized chi2
3133   //
3134   Float_t * chi2        = new Float_t[entries];
3135   Int_t * index         = new Int_t[entries];
3136   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3137   for (Int_t itrack=0;itrack<entries;itrack++){
3138     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3139     if (track){
3140       track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3141       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3142         chi2[itrack] = track->GetChi2MIP(0);
3143       else{
3144         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3145           delete array->RemoveAt(itrack);
3146         }
3147       }
3148     }
3149   }
3150   //
3151   TMath::Sort(entries,chi2,index,kFALSE);
3152   besttrack = (AliITStrackMI*)array->At(index[0]);
3153   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3154     for (Int_t j=0;j<6;j++){
3155       if (besttrack->GetClIndex(j)>0){
3156         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3157         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3158         ny[j]   = besttrack->GetNy(j);
3159         nz[j]   = besttrack->GetNz(j);
3160       }
3161     }
3162   }
3163   //
3164   // calculate one more time with updated normalized errors
3165   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3166   for (Int_t itrack=0;itrack<entries;itrack++){
3167     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3168     if (track){
3169       track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3170       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3171         chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3172       else
3173         {
3174           if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3175             delete array->RemoveAt(itrack);
3176           }
3177         }
3178     }
3179   }
3180   entries = array->GetEntriesFast();
3181   //
3182   //
3183   if (entries>0){
3184     TObjArray * newarray = new TObjArray();
3185     TMath::Sort(entries,chi2,index,kFALSE);
3186     besttrack = (AliITStrackMI*)array->At(index[0]);
3187     if (besttrack){
3188       //
3189       for (Int_t j=0;j<6;j++){
3190         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3191           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3192           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3193           ny[j]   = besttrack->GetNy(j);
3194           nz[j]   = besttrack->GetNz(j);
3195         }
3196       }
3197       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3198       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3199       Float_t minn = besttrack->GetNumberOfClusters()-3;
3200       Int_t accepted=0;
3201       for (Int_t i=0;i<entries;i++){
3202         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3203         if (!track) continue;
3204         if (accepted>maxcut) break;
3205         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3206         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3207           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3208             delete array->RemoveAt(index[i]);
3209             continue;
3210           }
3211         }
3212         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3213         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3214           if (!shortbest) accepted++;
3215           //
3216           newarray->AddLast(array->RemoveAt(index[i]));
3217           for (Int_t j=0;j<6;j++){
3218             if (nz[j]==0){
3219               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3220               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3221               ny[j]   = track->GetNy(j);
3222               nz[j]   = track->GetNz(j);
3223             }
3224           }
3225         }
3226         else{
3227           delete array->RemoveAt(index[i]);
3228         }
3229       }
3230       array->Delete();
3231       delete fTrackHypothesys.RemoveAt(esdindex);
3232       fTrackHypothesys.AddAt(newarray,esdindex);
3233     }
3234     else{
3235       array->Delete();
3236       delete fTrackHypothesys.RemoveAt(esdindex);
3237     }
3238   }
3239   delete [] chi2;
3240   delete [] index;
3241 }
3242 //------------------------------------------------------------------------
3243 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3244 {
3245   //-------------------------------------------------------------
3246   // try to find best hypothesy
3247   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3248   //-------------------------------------------------------------
3249   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3250   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3251   if (!array) return 0;
3252   Int_t entries = array->GetEntriesFast();
3253   if (!entries) return 0;
3254   Float_t minchi2 = 100000;
3255   AliITStrackMI * besttrack=0;
3256   //
3257   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3258   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3259   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3260   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3261   //
3262   for (Int_t i=0;i<entries;i++){
3263     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3264     if (!track) continue;
3265     Float_t sigmarfi,sigmaz;
3266     GetDCASigma(track,sigmarfi,sigmaz);
3267     track->SetDnorm(0,sigmarfi);
3268     track->SetDnorm(1,sigmaz);
3269     //
3270     track->SetChi2MIP(1,1000000);
3271     track->SetChi2MIP(2,1000000);
3272     track->SetChi2MIP(3,1000000);
3273     //
3274     // backtrack
3275     backtrack = new(backtrack) AliITStrackMI(*track);
3276     if (track->GetConstrain()) {
3277       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3278       if (!backtrack->Improve(0,xyzVtx,ersVtx))         continue;
3279       backtrack->ResetCovariance(10.);
3280     }else{
3281       backtrack->ResetCovariance(10.);
3282     }
3283     backtrack->ResetClusters();
3284
3285     Double_t x = original->GetX();
3286     if (!RefitAt(x,backtrack,track)) continue;
3287     //
3288     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3289     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3290     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3291     track->SetChi22(GetMatchingChi2(backtrack,original));
3292
3293     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3294     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3295     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3296
3297
3298     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3299     //
3300     //forward track - without constraint
3301     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3302     forwardtrack->ResetClusters();
3303     x = track->GetX();
3304     RefitAt(x,forwardtrack,track);
3305     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3306     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3307     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3308
3309     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3310     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3311     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3312     forwardtrack->SetD(0,track->GetD(0));
3313     forwardtrack->SetD(1,track->GetD(1));
3314     {
3315       Int_t list[6];
3316       AliITSRecPoint* clist[6];
3317       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3318       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3319     }
3320
3321     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3322     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3323     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3324       track->SetChi2MIP(3,1000);
3325       continue;
3326     }
3327     Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3328     //
3329     for (Int_t ichi=0;ichi<5;ichi++){
3330       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3331     }
3332     if (chi2 < minchi2){
3333       //besttrack = new AliITStrackMI(*forwardtrack);
3334       besttrack = track;
3335       besttrack->SetLabel(track->GetLabel());
3336       besttrack->SetFakeRatio(track->GetFakeRatio());
3337       minchi2   = chi2;
3338       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3339       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3340       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3341     }
3342   }
3343   delete backtrack;
3344   delete forwardtrack;
3345   Int_t accepted=0;
3346   for (Int_t i=0;i<entries;i++){
3347     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3348
3349     if (!track) continue;
3350
3351     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3352         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3353         track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3354       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3355         delete array->RemoveAt(i);
3356         continue;
3357       }
3358     }
3359     else{
3360       accepted++;
3361     }
3362   }
3363   //
3364   array->Compress();
3365   SortTrackHypothesys(esdindex,checkmax,1);
3366
3367   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3368   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3369   besttrack = (AliITStrackMI*)array->At(0);
3370   if (!besttrack)  return 0;
3371   besttrack->SetChi2MIP(8,0);
3372   fBestTrackIndex[esdindex]=0;
3373   entries = array->GetEntriesFast();
3374   AliITStrackMI *longtrack =0;
3375   minchi2 =1000;
3376   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3377   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3378     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3379     if (!track->GetConstrain()) continue;
3380     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3381     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3382     if (track->GetChi2MIP(0)>4.) continue;
3383     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3384     longtrack =track;
3385   }
3386   //if (longtrack) besttrack=longtrack;
3387
3388   Int_t list[6];
3389   AliITSRecPoint * clist[6];
3390   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3391   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3392       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3393     RegisterClusterTracks(besttrack,esdindex);
3394   }
3395   //
3396   //
3397   if (shared>0.0){
3398     Int_t nshared;
3399     Int_t overlist[6];
3400     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3401     if (sharedtrack>=0){
3402       //
3403       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3404       if (besttrack){
3405         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3406       }
3407       else return 0;
3408     }
3409   }
3410
3411   if (shared>2.5) return 0;
3412   if (shared>1.0) return besttrack;
3413   //
3414   // Don't sign clusters if not gold track
3415   //
3416   if (!besttrack->IsGoldPrimary()) return besttrack;
3417   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3418   //
3419   if (fConstraint[fPass]){
3420     //
3421     // sign clusters
3422     //
3423     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3424     for (Int_t i=0;i<6;i++){
3425       Int_t index = besttrack->GetClIndex(i);
3426       if (index<=0) continue;
3427       Int_t ilayer =  (index & 0xf0000000) >> 28;
3428       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3429       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3430       if (!c) continue;
3431       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3432       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3433       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3434       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3435         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3436       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3437
3438       Bool_t cansign = kTRUE;
3439       for (Int_t itrack=0;itrack<entries; itrack++){
3440         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3441         if (!track) continue;
3442         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3443         if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3444           cansign = kFALSE;
3445           break;
3446         }
3447       }
3448       if (cansign){
3449         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3450         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3451         if (!c->IsUsed()) c->Use();
3452       }
3453     }
3454   }
3455   return besttrack;
3456 }
3457 //------------------------------------------------------------------------
3458 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3459 {
3460   //
3461   // get "best" hypothesys
3462   //
3463
3464   Int_t nentries = itsTracks.GetEntriesFast();
3465   for (Int_t i=0;i<nentries;i++){
3466     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3467     if (!track) continue;
3468     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3469     if (!array) continue;
3470     if (array->GetEntriesFast()<=0) continue;
3471     //
3472     AliITStrackMI* longtrack=0;
3473     Float_t minn=0;
3474     Float_t maxchi2=1000;
3475     for (Int_t j=0;j<array->GetEntriesFast();j++){
3476       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3477       if (!trackHyp) continue;
3478       if (trackHyp->GetGoldV0()) {
3479         longtrack = trackHyp;   //gold V0 track taken
3480         break;
3481       }
3482       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3483       Float_t chi2 = trackHyp->GetChi2MIP(0);
3484       if (fAfterV0){
3485         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3486       }
3487       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3488       //
3489       if (chi2 > maxchi2) continue;
3490       minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3491       maxchi2 = chi2;
3492       longtrack=trackHyp;
3493     }
3494     //
3495     //
3496     //
3497     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3498     if (!longtrack) {longtrack = besttrack;}
3499     else besttrack= longtrack;
3500     //
3501     if (besttrack) {
3502       Int_t list[6];
3503       AliITSRecPoint * clist[6];
3504       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3505       //
3506       track->SetNUsed(shared);
3507       track->SetNSkipped(besttrack->GetNSkipped());
3508       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3509       if (shared>0) {
3510         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3511         Int_t nshared;
3512         Int_t overlist[6];
3513         //
3514         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3515         //if (sharedtrack==-1) sharedtrack=0;
3516         if (sharedtrack>=0) {
3517           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3518         }
3519       }
3520       if (besttrack&&fAfterV0) {
3521         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3522       }
3523       if (besttrack&&fConstraint[fPass])
3524         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3525       if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3526         if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3527              TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3528       }
3529
3530     }
3531   }
3532 }
3533 //------------------------------------------------------------------------
3534 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3535   //--------------------------------------------------------------------
3536   //This function "cooks" a track label. If label<0, this track is fake.
3537   //--------------------------------------------------------------------
3538   Int_t tpcLabel=-1;
3539
3540   if ( track->GetESDtrack())   tpcLabel =  TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3541
3542    track->SetChi2MIP(9,0);
3543    Int_t nwrong=0;
3544    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3545      Int_t cindex = track->GetClusterIndex(i);
3546      Int_t l=(cindex & 0xf0000000) >> 28;
3547      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3548      Int_t isWrong=1;
3549      for (Int_t ind=0;ind<3;ind++){
3550        if (tpcLabel>0)
3551          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3552        AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3553      }
3554      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3555      nwrong+=isWrong;
3556    }
3557    Int_t nclusters = track->GetNumberOfClusters();
3558    if (nclusters > 0) //PH Some tracks don't have any cluster
3559      track->SetFakeRatio(double(nwrong)/double(nclusters));
3560    if (tpcLabel>0){
3561      if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3562      else
3563        track->SetLabel(tpcLabel);
3564    }
3565    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3566
3567 }
3568 //------------------------------------------------------------------------
3569 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3570 {
3571   //
3572   //
3573   //  Int_t list[6];
3574   //AliITSRecPoint * clist[6];
3575   //  Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3576   Float_t dedx[4];
3577   Int_t accepted=0;
3578   track->SetChi2MIP(9,0);
3579   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3580     Int_t cindex = track->GetClusterIndex(i);
3581     Int_t l=(cindex & 0xf0000000) >> 28;
3582     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3583     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3584     Int_t isWrong=1;
3585     for (Int_t ind=0;ind<3;ind++){
3586       if (cl->GetLabel(ind)==lab) isWrong=0;
3587     }
3588     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3589     if (l<2) continue;
3590     //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue;  //shared track
3591     //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3592     //if (l<4&& !(cl->GetType()==1)) continue;
3593     dedx[accepted]= track->GetSampledEdx(i);
3594     //dedx[accepted]= track->fNormQ[l];
3595     accepted++;
3596   }
3597   if (accepted<1) {
3598     track->SetdEdx(0);
3599     return;
3600   }
3601   Int_t indexes[4];
3602   TMath::Sort(accepted,dedx,indexes,kFALSE);
3603   Double_t low=0.;
3604   Double_t up=0.51;
3605   Double_t nl=low*accepted, nu =up*accepted;
3606   Float_t sumamp = 0;
3607   Float_t sumweight =0;
3608   for (Int_t i=0; i<accepted; i++) {
3609     Float_t weight =1;
3610     if (i<nl+0.1)     weight = TMath::Max(1.-(nl-i),0.);
3611     if (i>nu-1)     weight = TMath::Max(nu-i,0.);
3612     sumamp+= dedx[indexes[i]]*weight;
3613     sumweight+=weight;
3614   }
3615   track->SetdEdx(sumamp/sumweight);
3616 }
3617 //------------------------------------------------------------------------
3618 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3619   //
3620   //
3621   if (fCoefficients) delete []fCoefficients;
3622   fCoefficients = new Float_t[ntracks*48];
3623   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3624 }
3625 //------------------------------------------------------------------------
3626 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3627 {
3628   //
3629   //
3630   //
3631   Float_t erry,errz;
3632   Float_t theta = track->GetTgl();
3633   Float_t phi   = track->GetSnp();
3634   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3635   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3636   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()));
3637   // Take into account the mis-alignment (bring track to cluster plane)
3638   Double_t xTrOrig=track->GetX();
3639   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3640   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()));
3641   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3642   // Bring the track back to detector plane in ideal geometry
3643   // [mis-alignment will be accounted for in UpdateMI()]
3644   if (!track->Propagate(xTrOrig)) return 1000.;
3645   Float_t ny,nz;
3646   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3647   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3648   if (delta>1){
3649     chi2+=0.5*TMath::Min(delta/2,2.);
3650     chi2+=2.*cluster->GetDeltaProbability();
3651   }
3652   //
3653   track->SetNy(layer,ny);
3654   track->SetNz(layer,nz);
3655   track->SetSigmaY(layer,erry);
3656   track->SetSigmaZ(layer, errz);
3657   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3658   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3659   return chi2;
3660
3661 }
3662 //------------------------------------------------------------------------
3663 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3664 {
3665   //
3666   //
3667   //
3668   Int_t layer = (index & 0xf0000000) >> 28;
3669   track->SetClIndex(layer, index);
3670   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3671     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3672       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3673       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3674     }
3675   }
3676
3677   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
3678
3679
3680   // Take into account the mis-alignment (bring track to cluster plane)
3681   Double_t xTrOrig=track->GetX();
3682   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3683   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3684   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3685   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
3686
3687   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3688
3689
3690   AliCluster c(*cl);
3691   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3692   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3693
3694
3695   Int_t updated = track->UpdateMI(&c,chi2,index);
3696
3697   // Bring the track back to detector plane in ideal geometry
3698   if (!track->Propagate(xTrOrig)) return 0;
3699
3700   if(!updated) AliDebug(2,"update failed");
3701   return updated;
3702 }
3703
3704 //------------------------------------------------------------------------
3705 void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3706 {
3707   //
3708   //DCA sigmas parameterization
3709   //to be paramterized using external parameters in future
3710   //
3711   //
3712   sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3713   sigmaz   = 0.0110+4.37*TMath::Abs(track->GetC());
3714 }
3715 //------------------------------------------------------------------------
3716 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3717 {
3718   //
3719   //
3720   Int_t entries = ClusterArray->GetEntriesFast();
3721   if (entries<4) return;
3722   AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3723   Int_t layer = cluster->GetLayer();
3724   if (layer>1) return;
3725   Int_t index[10000];
3726   Int_t ncandidates=0;
3727   Float_t r = (layer>0)? 7:4;
3728   //
3729   for (Int_t i=0;i<entries;i++){
3730     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3731     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3732     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3733     index[ncandidates] = i;  //candidate to belong to delta electron track
3734     ncandidates++;
3735     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3736       cl0->SetDeltaProbability(1);
3737     }
3738   }
3739   //
3740   //
3741   //
3742   for (Int_t i=0;i<ncandidates;i++){
3743     AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3744     if (cl0->GetDeltaProbability()>0.8) continue;
3745     //
3746     Int_t ncl = 0;
3747     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3748     sumy=sumz=sumy2=sumyz=sumw=0.0;
3749     for (Int_t j=0;j<ncandidates;j++){
3750       if (i==j) continue;
3751       AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3752       //
3753       Float_t dz = cl0->GetZ()-cl1->GetZ();
3754       Float_t dy = cl0->GetY()-cl1->GetY();
3755       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3756         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3757         y[ncl] = cl1->GetY();
3758         z[ncl] = cl1->GetZ();
3759         sumy+= y[ncl]*weight;
3760         sumz+= z[ncl]*weight;
3761         sumy2+=y[ncl]*y[ncl]*weight;
3762         sumyz+=y[ncl]*z[ncl]*weight;
3763         sumw+=weight;
3764         ncl++;
3765       }
3766     }
3767     if (ncl<4) continue;
3768     Float_t det = sumw*sumy2  - sumy*sumy;
3769     Float_t delta=1000;
3770     if (TMath::Abs(det)>0.01){
3771       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3772       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3773       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3774     }
3775     else{
3776       Float_t z0  = sumyz/sumy;
3777       delta = TMath::Abs(cl0->GetZ()-z0);
3778     }
3779     if ( delta<0.05) {
3780       cl0->SetDeltaProbability(1-20.*delta);
3781     }
3782   }
3783 }
3784 //------------------------------------------------------------------------
3785 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3786 {
3787   //
3788   //
3789   track->UpdateESDtrack(flags);
3790   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3791   if (oldtrack) delete oldtrack;
3792   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3793   if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3794     printf("Problem\n");
3795   }
3796 }
3797 //------------------------------------------------------------------------
3798 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3799   //
3800   // Get nearest upper layer close to the point xr.
3801   // rough approximation
3802   //
3803   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3804   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3805   Int_t res =6;
3806   for (Int_t i=0;i<6;i++){
3807     if (radius<kRadiuses[i]){
3808       res =i;
3809       break;
3810     }
3811   }
3812   return res;
3813 }
3814 //------------------------------------------------------------------------
3815 void AliITStrackerMI::UpdateTPCV0(AliESDEvent *event){
3816   //
3817   //try to update, or reject TPC  V0s
3818   //
3819   Int_t nv0s = event->GetNumberOfV0s();
3820   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3821
3822   for (Int_t i=0;i<nv0s;i++){
3823     AliESDv0 * vertex = event->GetV0(i);
3824     Int_t ip = vertex->GetIndex(0);
3825     Int_t im = vertex->GetIndex(1);
3826     //
3827     TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3828     TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3829     AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3830     AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3831     //
3832     //
3833     if (trackp){
3834       if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3835         if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3836         if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3837       }
3838     }
3839
3840     if (trackm){
3841       if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3842         if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3843         if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3844       }
3845     }
3846     if (vertex->GetStatus()==-100) continue;
3847     //
3848     Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
3849     Int_t clayer = GetNearestLayer(xrp);                    //I.B.
3850     vertex->SetNBefore(clayer);        //
3851     vertex->SetChi2Before(9*clayer);   //
3852     vertex->SetNAfter(6-clayer);       //
3853     vertex->SetChi2After(0);           //
3854     //
3855     if (clayer >1 ){ // calculate chi2 before vertex
3856       Float_t chi2p = 0, chi2m=0;
3857       //
3858       if (trackp){
3859         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3860           if (trackp->GetClIndex(ilayer)>0){
3861             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3862               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3863           }
3864           else{
3865             chi2p+=9;
3866           }
3867         }
3868       }else{
3869         chi2p = 9*clayer;
3870       }
3871       //
3872       if (trackm){
3873         for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3874           if (trackm->GetClIndex(ilayer)>0){
3875             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3876               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3877           }
3878           else{
3879             chi2m+=9;
3880           }
3881         }
3882       }else{
3883         chi2m = 9*clayer;
3884       }
3885       vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3886       if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10);  // track exist before vertex
3887     }
3888
3889     if (clayer < 5 ){ // calculate chi2 after vertex
3890       Float_t chi2p = 0, chi2m=0;
3891       //
3892       if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3893         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3894           if (trackp->GetClIndex(ilayer)>0){
3895             chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3896               trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3897           }
3898           else{
3899             chi2p+=9;
3900           }
3901         }
3902       }else{
3903         chi2p = 0;
3904       }
3905       //
3906       if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3907         for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3908           if (trackm->GetClIndex(ilayer)>0){
3909             chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3910               trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3911           }
3912           else{
3913             chi2m+=9;
3914           }
3915         }
3916       }else{
3917         chi2m = 0;
3918       }
3919       vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3920       if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20);  // track not found in ITS
3921     }
3922   }
3923   //
3924 }
3925 //------------------------------------------------------------------------
3926 void AliITStrackerMI::FindV02(AliESDEvent *event)
3927 {
3928   //
3929   // V0 finder
3930   //
3931   //  Cuts on DCA -  R dependent
3932   //          max distance DCA between 2 tracks cut
3933   //          maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3934   //
3935   const Float_t kMaxDist0      = 0.1;
3936   const Float_t kMaxDist1      = 0.1;
3937   const Float_t kMaxDist       = 1;
3938   const Float_t kMinPointAngle  = 0.85;
3939   const Float_t kMinPointAngle2 = 0.99;
3940   const Float_t kMinR           = 0.5;
3941   const Float_t kMaxR           = 220;
3942   //const Float_t kCausality0Cut   = 0.19;
3943   //const Float_t kLikelihood01Cut = 0.25;
3944   //const Float_t kPointAngleCut   = 0.9996;
3945   const Float_t kCausality0Cut   = 0.19;
3946   const Float_t kLikelihood01Cut = 0.45;
3947   const Float_t kLikelihood1Cut  = 0.5;
3948   const Float_t kCombinedCut     = 0.55;
3949
3950   //
3951   //
3952   TTreeSRedirector &cstream = *fDebugStreamer;
3953   Int_t ntracks    = event->GetNumberOfTracks();
3954   Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3955   fOriginal.Expand(ntracks);
3956   fTrackHypothesys.Expand(ntracks);
3957   fBestHypothesys.Expand(ntracks);
3958   //
3959   AliHelix * helixes   = new AliHelix[ntracks+2];
3960   TObjArray trackarray(ntracks+2);     //array with tracks - with vertex constrain
3961   TObjArray trackarrayc(ntracks+2);    //array of "best    tracks" - without vertex constrain
3962   TObjArray trackarrayl(ntracks+2);    //array of "longest tracks" - without vertex constrain
3963   Bool_t * forbidden   = new Bool_t [ntracks+2];
3964   Int_t   *itsmap      = new Int_t  [ntracks+2];
3965   Float_t *dist        = new Float_t[ntracks+2];
3966   Float_t *normdist0   = new Float_t[ntracks+2];
3967   Float_t *normdist1   = new Float_t[ntracks+2];
3968   Float_t *normdist    = new Float_t[ntracks+2];
3969   Float_t *norm        = new Float_t[ntracks+2];
3970   Float_t *maxr        = new Float_t[ntracks+2];
3971   Float_t *minr        = new Float_t[ntracks+2];
3972   Float_t *minPointAngle= new Float_t[ntracks+2];
3973   //
3974   AliV0 *pvertex      = new AliV0;
3975   AliITStrackMI * dummy= new AliITStrackMI;
3976   dummy->SetLabel(0);
3977   AliITStrackMI  trackat0;    //temporary track for DCA calculation
3978   //
3979   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3980   //
3981   // make ITS -  ESD map
3982   //
3983   for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3984     itsmap[itrack]        = -1;
3985     forbidden[itrack]     = kFALSE;
3986     maxr[itrack]          = kMaxR;
3987     minr[itrack]          = kMinR;
3988     minPointAngle[itrack] = kMinPointAngle;
3989   }
3990   for (Int_t itrack=0;itrack<nitstracks;itrack++){
3991     AliITStrackMI * original =   (AliITStrackMI*)(fOriginal.At(itrack));
3992     Int_t           esdindex =   original->GetESDtrack()->GetID();
3993     itsmap[esdindex]         =   itrack;
3994   }
3995   //
3996   // create ITS tracks from ESD tracks if not done before
3997   //
3998   for (Int_t itrack=0;itrack<ntracks;itrack++){
3999     if (itsmap[itrack]>=0) continue;
4000     AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
4001     //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
4002     //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
4003     tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP());   //I.B.
4004     if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
4005       // tracks which can reach inner part of ITS
4006       // propagate track to outer its volume - with correction for material
4007       CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
4008     }
4009     itsmap[itrack] = nitstracks;
4010     fOriginal.AddAt(tpctrack,nitstracks);
4011     nitstracks++;
4012   }
4013   //
4014   // fill temporary arrays
4015   //
4016   for (Int_t itrack=0;itrack<ntracks;itrack++){
4017     AliESDtrack *  esdtrack = event->GetTrack(itrack);
4018     Int_t          itsindex = itsmap[itrack];
4019     AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
4020     if (!original) continue;
4021     AliITStrackMI *bestConst  = 0;
4022     AliITStrackMI *bestLong   = 0;
4023     AliITStrackMI *best       = 0;
4024     //
4025     //
4026     TObjArray * array    = (TObjArray*)  fTrackHypothesys.At(itsindex);
4027     Int_t       hentries = (array==0) ?  0 : array->GetEntriesFast();
4028     // Get best track with vertex constrain
4029     for (Int_t ih=0;ih<hentries;ih++){
4030       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4031       if (!trackh->GetConstrain()) continue;
4032       if (!bestConst) bestConst = trackh;
4033       if (trackh->GetNumberOfClusters()>5.0){
4034         bestConst  = trackh;                         // full track -  with minimal chi2
4035         break;
4036       }
4037       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone())  continue;
4038       bestConst = trackh;
4039       break;
4040     }
4041     // Get best long track without vertex constrain and best track without vertex constrain
4042     for (Int_t ih=0;ih<hentries;ih++){
4043       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4044       if (trackh->GetConstrain()) continue;
4045       if (!best)     best     = trackh;
4046       if (!bestLong) bestLong = trackh;
4047       if (trackh->GetNumberOfClusters()>5.0){
4048         bestLong  = trackh;                         // full track -  with minimal chi2
4049         break;
4050       }
4051       if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone())  continue;
4052       bestLong = trackh;
4053     }
4054     if (!best) {
4055       best     = original;
4056       bestLong = original;
4057     }
4058     //I.B. trackat0 = *bestLong;
4059     new (&trackat0) AliITStrackMI(*bestLong);
4060     Double_t xx,yy,zz,alpha;
4061     if (!bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz)) continue;
4062     alpha = TMath::ATan2(yy,xx);
4063     if (!trackat0.Propagate(alpha,0)) continue;
4064     // calculate normalized distances to the vertex
4065     //
4066     Float_t ptfac  = (1.+100.*TMath::Abs(trackat0.GetC()));
4067     if ( bestLong->GetNumberOfClusters()>3 ){
4068       dist[itsindex]      = trackat0.GetY();
4069       norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4070       normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4071       normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4072       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4073       if (!bestConst){
4074         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
4075         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
4076         if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
4077       }else{
4078         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
4079         if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
4080       }
4081     }
4082     else{
4083       if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
4084         dist[itsindex] = bestConst->GetD(0);
4085         norm[itsindex] = bestConst->GetDnorm(0);
4086         normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4087         normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4088         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4089       }else{
4090         dist[itsindex]      = trackat0.GetY();
4091         norm[itsindex]      = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4092         normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4093         normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4094         normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4095         if (TMath::Abs(trackat0.GetTgl())>1.05){
4096           if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
4097           if (normdist[itsindex]>3) {
4098             minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
4099           }
4100         }
4101       }
4102     }
4103     //
4104     //-----------------------------------------------------------
4105     //Forbid primary track candidates -
4106     //
4107     //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
4108     //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
4109     //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
4110     //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
4111     //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
4112     //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
4113     //-----------------------------------------------------------
4114     if (bestConst){
4115       if (bestLong->GetNumberOfClusters()<4       && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5)               forbidden[itsindex]=kTRUE;
4116       if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5)               forbidden[itsindex]=kTRUE;
4117       if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
4118       if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0)                              forbidden[itsindex]=kTRUE;
4119       if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2)                             forbidden[itsindex]=kTRUE;
4120       if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1)                             forbidden[itsindex]=kTRUE;
4121       if (bestConst->GetNormChi2(0)<2.5) {
4122         minPointAngle[itsindex]= 0.9999;
4123         maxr[itsindex]         = 10;
4124       }
4125     }
4126     //
4127     //forbid daughter kink candidates
4128     //
4129     if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
4130     Bool_t isElectron = kTRUE;
4131     Bool_t isProton   = kTRUE;
4132     Double_t pid[5];
4133     esdtrack->GetESDpid(pid);
4134     for (Int_t i=1;i<5;i++){
4135       if (pid[0]<pid[i]) isElectron= kFALSE;
4136       if (pid[4]<pid[i]) isProton= kFALSE;
4137     }
4138     if (isElectron){
4139       forbidden[itsindex]=kFALSE;
4140       normdist[itsindex]*=-1;
4141     }
4142     if (isProton){
4143       if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
4144       normdist[itsindex]*=-1;
4145     }
4146
4147     //
4148     // Causality cuts in TPC volume
4149     //
4150     if (esdtrack->GetTPCdensity(0,10) >0.6)  maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
4151     if (esdtrack->GetTPCdensity(10,30)>0.6)  maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
4152     if (esdtrack->GetTPCdensity(20,40)>0.6)  maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
4153     if (esdtrack->GetTPCdensity(30,50)>0.6)  maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
4154     //
4155     if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
4156     //
4157     //
4158     if (kFALSE){
4159       cstream<<"Track"<<
4160         "Tr0.="<<best<<
4161         "Tr1.="<<((bestConst)? bestConst:dummy)<<
4162         "Tr2.="<<bestLong<<
4163         "Tr3.="<<&trackat0<<
4164         "Esd.="<<esdtrack<<
4165         "Dist="<<dist[itsindex]<<
4166         "ND0="<<normdist0[itsindex]<<
4167         "ND1="<<normdist1[itsindex]<<
4168         "ND="<<normdist[itsindex]<<
4169         "Pz="<<primvertex[2]<<
4170         "Forbid="<<forbidden[itsindex]<<
4171         "\n";
4172       //
4173     }
4174     trackarray.AddAt(best,itsindex);
4175     trackarrayc.AddAt(bestConst,itsindex);
4176     trackarrayl.AddAt(bestLong,itsindex);
4177     new (&helixes[itsindex]) AliHelix(*best);
4178   }
4179   //
4180   //
4181   //
4182   // first iterration of V0 finder
4183   //
4184   for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
4185     Int_t itrack0 = itsmap[iesd0];
4186     if (forbidden[itrack0]) continue;
4187     AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
4188     if (!btrack0) continue;
4189     if (btrack0->GetSign()>0) continue;
4190     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
4191     //
4192     for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
4193       Int_t itrack1 = itsmap[iesd1];
4194       if (forbidden[itrack1]) continue;
4195
4196       AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
4197       if (!btrack1) continue;
4198       if (btrack1->GetSign()<0) continue;
4199       Bool_t isGold = kFALSE;
4200       if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
4201         isGold = kTRUE;
4202       }
4203       AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
4204       AliHelix &h1 = helixes[itrack0];
4205       AliHelix &h2 = helixes[itrack1];
4206       //
4207       // find linear distance
4208       Double_t rmin =0;
4209       //
4210       //
4211       //
4212       Double_t phase[2][2],radius[2];
4213       Int_t  points = h1.GetRPHIintersections(h2, phase, radius);
4214       if    (points==0)  continue;
4215       Double_t delta[2]={1000000,1000000};
4216       rmin = radius[0];
4217       h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
4218       if (points==2){
4219         if (radius[1]<rmin) rmin = radius[1];
4220         h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
4221       }
4222       rmin = TMath::Sqrt(rmin);
4223       Double_t distance = 0;
4224       Double_t radiusC  = 0;
4225       Int_t    iphase   = 0;
4226       if (points==1 || delta[0]<delta[1]){
4227         distance = TMath::Sqrt(delta[0]);
4228         radiusC  = TMath::Sqrt(radius[0]);
4229       }else{
4230         distance = TMath::Sqrt(delta[1]);
4231         radiusC  = TMath::Sqrt(radius[1]);
4232         iphase=1;
4233       }
4234       if (radiusC<TMath::Max(minr[itrack0],minr[itrack1]))    continue;
4235       if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1]))     continue;
4236       Float_t maxDist  = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
4237       if (distance>maxDist) continue;
4238       Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
4239       if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
4240       //
4241       //
4242       //      Double_t distance = TestV0(h1,h2,pvertex,rmin);
4243       //
4244       //       if (distance>maxDist)           continue;
4245       //       if (pvertex->GetRr()<kMinR)     continue;
4246       //       if (pvertex->GetRr()>kMaxR)     continue;
4247       AliITStrackMI * track0=btrack0;
4248       AliITStrackMI * track1=btrack1;
4249       //      if (pvertex->GetRr()<3.5){
4250       if (radiusC<3.5){
4251         //use longest tracks inside the pipe
4252         track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4253         track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4254       }
4255       //
4256       //
4257       pvertex->SetParamN(*track0);
4258       pvertex->SetParamP(*track1);
4259       pvertex->Update(primvertex);
4260       pvertex->SetClusters(track0->ClIndex(),track1->ClIndex());  // register clusters
4261
4262       if (pvertex->GetRr()<kMinR) continue;
4263       if (pvertex->GetRr()>kMaxR) continue;
4264       if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4265 //Bo:      if (pvertex->GetDist2()>maxDist) continue;
4266       if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4267 //Bo:        pvertex->SetLab(0,track0->GetLabel());
4268 //Bo:        pvertex->SetLab(1,track1->GetLabel());
4269       pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4270       pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4271       //
4272       AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4273       AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4274
4275       //
4276       //
4277       TObjArray * array0b     = (TObjArray*)fBestHypothesys.At(itrack0);
4278       if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4279         fCurrentEsdTrack = itrack0;
4280         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4281       }
4282       TObjArray * array1b    = (TObjArray*)fBestHypothesys.At(itrack1);
4283       if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4284         fCurrentEsdTrack = itrack1;
4285         FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4286       }
4287       //
4288       AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4289       AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4290       AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4291       AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4292
4293       Float_t minchi2before0=16;
4294       Float_t minchi2before1=16;
4295       Float_t minchi2after0 =16;
4296       Float_t minchi2after1 =16;
4297       Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
4298       Int_t maxLayer = GetNearestLayer(xrp);                   //I.B.
4299
4300       if (array0b) for (Int_t i=0;i<5;i++){
4301         // best track after vertex
4302         AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4303         if (!btrack) continue;
4304         if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4305         //      if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4306         if (btrack->GetX()<pvertex->GetRr()-2.) {
4307           if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4308             Float_t sumchi2= 0;
4309             Float_t sumn   = 0;
4310             if (maxLayer<3){   // take prim vertex as additional measurement
4311               if (normdist[itrack0]>0 && htrackc0){
4312                 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4313               }else{
4314                 sumchi2 +=  TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4315               }
4316               sumn    +=  3-maxLayer;
4317             }
4318             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4319               sumn+=1.;
4320               if (!btrack->GetClIndex(ilayer)){
4321                 sumchi2+=25;
4322                 continue;
4323               }else{
4324                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4325                 for (Int_t itrack=0;itrack<4;itrack++){
4326                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4327                     sumchi2+=18.;  //shared cluster
4328                     break;
4329                   }
4330                 }
4331                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4332                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4333               }
4334             }
4335             sumchi2/=sumn;
4336             if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4337           }
4338           continue;   //safety space - Geo manager will give exact layer
4339         }
4340         track0b       = btrack;
4341         minchi2after0 = btrack->GetNormChi2(i);
4342         break;
4343       }
4344       if (array1b) for (Int_t i=0;i<5;i++){
4345         // best track after vertex
4346         AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4347         if (!btrack) continue;
4348         if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4349         //      if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4350         if (btrack->GetX()<pvertex->GetRr()-2){
4351           if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4352             Float_t sumchi2= 0;
4353             Float_t sumn   = 0;
4354             if (maxLayer<3){   // take prim vertex as additional measurement
4355               if (normdist[itrack1]>0 && htrackc1){
4356                 sumchi2 +=  TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4357               }else{
4358                 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4359               }
4360               sumn    +=  3-maxLayer;
4361             }
4362             for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4363               sumn+=1.;
4364               if (!btrack->GetClIndex(ilayer)){
4365                 sumchi2+=30;
4366                 continue;
4367               }else{
4368                 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4369                 for (Int_t itrack=0;itrack<4;itrack++){
4370                   if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4371                     sumchi2+=18.;  //shared cluster
4372                     break;
4373                   }
4374                 }
4375                 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4376                 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4377               }
4378             }
4379             sumchi2/=sumn;
4380             if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4381           }
4382           continue;   //safety space - Geo manager will give exact layer
4383         }
4384         track1b = btrack;
4385         minchi2after1 = btrack->GetNormChi2(i);
4386         break;
4387       }
4388       //
4389       // position resolution - used for DCA cut
4390       Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4391         (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4392         (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4393       sigmad =TMath::Sqrt(sigmad)+0.04;
4394       if (pvertex->GetRr()>50){
4395         Double_t cov0[15],cov1[15];
4396         track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4397         track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4398         sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4399           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4400           (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4401         sigmad =TMath::Sqrt(sigmad)+0.05;
4402       }
4403       //
4404       AliV0 vertex2;
4405       vertex2.SetParamN(*track0b);
4406       vertex2.SetParamP(*track1b);
4407       vertex2.Update(primvertex);
4408       //Bo:      if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4409       if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4410         pvertex->SetParamN(*track0b);
4411         pvertex->SetParamP(*track1b);
4412         pvertex->Update(primvertex);
4413         pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex());  // register clusters
4414         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4415         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4416       }
4417       pvertex->SetDistSigma(sigmad);
4418       //Bo:      pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4419       pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4420       //
4421       // define likelihhod and causalities
4422       //
4423       Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4424       if (maxLayer<1){
4425         Float_t fnorm0 = normdist[itrack0];
4426         if (fnorm0<0) fnorm0*=-3;
4427         Float_t fnorm1 = normdist[itrack1];
4428         if (fnorm1<0) fnorm1*=-3;
4429         if ((pvertex->GetAnglep()[2]>0.1) || ( (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 ) || (pvertex->GetRr()<3)){
4430           pb0    =  TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4431           pb1    =  TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4432         }
4433         pvertex->SetChi2Before(normdist[itrack0]);
4434         pvertex->SetChi2After(normdist[itrack1]);
4435         pvertex->SetNAfter(0);
4436         pvertex->SetNBefore(0);
4437       }else{
4438         pvertex->SetChi2Before(minchi2before0);
4439         pvertex->SetChi2After(minchi2before1);
4440          if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4441            pb0    =  TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4442            pb1    =  TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4443          }
4444          pvertex->SetNAfter(maxLayer);
4445          pvertex->SetNBefore(maxLayer);
4446       }
4447       if (pvertex->GetRr()<90){
4448         pa0  *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4449         pa1  *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4450       }
4451       if (pvertex->GetRr()<20){
4452         pa0  *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4453         pa1  *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4454       }
4455       //
4456       pvertex->SetCausality(pb0,pb1,pa0,pa1);
4457       //
4458       //  Likelihood calculations  - apply cuts
4459       //
4460       Bool_t v0OK = kTRUE;
4461       Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4462       p12        += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4463       p12         = TMath::Sqrt(p12);                             // "mean" momenta
4464       Float_t    sigmap0   = 0.0001+0.001/(0.1+pvertex->GetRr());
4465       Float_t    sigmap    = 0.5*sigmap0*(0.6+0.4*p12);           // "resolution: of point angle - as a function of radius and momenta
4466
4467       Float_t causalityA  = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4468       Float_t causalityB  = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4469                                         TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4470       //
4471       //Bo:      Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4472       Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4473       Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4474
4475       Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4476         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4477         0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4478         0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4479       //
4480       if (causalityA<kCausality0Cut)                                          v0OK = kFALSE;
4481       if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut)              v0OK = kFALSE;
4482       if (likelihood1<kLikelihood1Cut)                                        v0OK = kFALSE;
4483       if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4484
4485       //
4486       //
4487       if (kFALSE){
4488         Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4489         cstream<<"It0"<<
4490           "Tr0.="<<track0<<                       //best without constrain
4491           "Tr1.="<<track1<<                       //best without constrain
4492           "Tr0B.="<<track0b<<                     //best without constrain  after vertex
4493           "Tr1B.="<<track1b<<                     //best without constrain  after vertex
4494           "Tr0C.="<<htrackc0<<                    //best with constrain     if exist
4495           "Tr1C.="<<htrackc1<<                    //best with constrain     if exist
4496           "Tr0L.="<<track0l<<                     //longest best
4497           "Tr1L.="<<track1l<<                     //longest best
4498           "Esd0.="<<track0->GetESDtrack()<<           // esd track0 params
4499           "Esd1.="<<track1->GetESDtrack()<<           // esd track1 params
4500           "V0.="<<pvertex<<                       //vertex properties
4501           "V0b.="<<&vertex2<<                       //vertex properties at "best" track
4502           "ND0="<<normdist[itrack0]<<             //normalize distance for track0
4503           "ND1="<<normdist[itrack1]<<             //normalize distance for track1
4504           "Gold="<<gold<<                         //
4505           //      "RejectBase="<<rejectBase<<             //rejection in First itteration
4506           "OK="<<v0OK<<
4507           "rmin="<<rmin<<
4508           "sigmad="<<sigmad<<
4509           "\n";
4510       }
4511       //if (rejectBase) continue;
4512       //
4513       pvertex->SetStatus(0);
4514       //      if (rejectBase) {
4515       //        pvertex->SetStatus(-100);
4516       //}
4517       if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4518         //Bo:     pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4519         pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4520         pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4521         if (v0OK){
4522           //      AliV0vertex vertexjuri(*track0,*track1);
4523           //      vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4524           //      event->AddV0(&vertexjuri);
4525           pvertex->SetStatus(100);
4526         }
4527         pvertex->SetOnFlyStatus(kTRUE);
4528         pvertex->ChangeMassHypothesis(kK0Short);
4529         event->AddV0(pvertex);
4530       }
4531     }
4532   }
4533   //
4534   //
4535   // delete temporary arrays
4536   //
4537   delete[] forbidden;
4538   delete[] minPointAngle;
4539   delete[] maxr;
4540   delete[] minr;
4541   delete[] norm;
4542   delete[] normdist;
4543   delete[] normdist1;
4544   delete[] normdist0;
4545   delete[] dist;
4546   delete[] itsmap;
4547   delete[] helixes;
4548   delete   pvertex;
4549 }
4550 //------------------------------------------------------------------------
4551 void AliITStrackerMI::RefitV02(AliESDEvent *event)
4552 {
4553   //
4554   //try to refit  V0s in the third path of the reconstruction
4555   //
4556   TTreeSRedirector &cstream = *fDebugStreamer;
4557   //
4558   Int_t  nv0s = event->GetNumberOfV0s();
4559   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4560   AliV0 v0temp;
4561   for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4562     AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4563     if (!v0mi) continue;
4564     Int_t     itrack0   = v0mi->GetIndex(0);
4565     Int_t     itrack1   = v0mi->GetIndex(1);
4566     AliESDtrack *esd0   = event->GetTrack(itrack0);
4567     AliESDtrack *esd1   = event->GetTrack(itrack1);
4568     if (!esd0||!esd1) continue;
4569     AliITStrackMI tpc0(*esd0);
4570     AliITStrackMI tpc1(*esd1);
4571     Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4572     Double_t alpha =TMath::ATan2(y,x);   //I.B.
4573     if (v0mi->GetRr()>85){
4574       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4575         v0temp.SetParamN(tpc0);
4576         v0temp.SetParamP(tpc1);
4577         v0temp.Update(primvertex);
4578         if (kFALSE) cstream<<"Refit"<<
4579           "V0.="<<v0mi<<
4580           "V0refit.="<<&v0temp<<
4581           "Tr0.="<<&tpc0<<
4582           "Tr1.="<<&tpc1<<
4583           "\n";
4584         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4585         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4586           v0mi->SetParamN(tpc0);
4587           v0mi->SetParamP(tpc1);
4588           v0mi->Update(primvertex);
4589         }
4590       }
4591       continue;
4592     }
4593     if (v0mi->GetRr()>35){
4594       CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4595       CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4596       if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4597         v0temp.SetParamN(tpc0);
4598         v0temp.SetParamP(tpc1);
4599         v0temp.Update(primvertex);
4600         if (kFALSE) cstream<<"Refit"<<
4601           "V0.="<<v0mi<<
4602           "V0refit.="<<&v0temp<<
4603           "Tr0.="<<&tpc0<<
4604           "Tr1.="<<&tpc1<<
4605           "\n";
4606         //Bo:   if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4607         if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4608           v0mi->SetParamN(tpc0);
4609           v0mi->SetParamP(tpc1);
4610           v0mi->Update(primvertex);
4611         }
4612       }
4613       continue;
4614     }
4615     CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4616     CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4617     //    if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4618     if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4619       v0temp.SetParamN(tpc0);
4620       v0temp.SetParamP(tpc1);
4621       v0temp.Update(primvertex);
4622       if (kFALSE) cstream<<"Refit"<<
4623         "V0.="<<v0mi<<
4624         "V0refit.="<<&v0temp<<
4625         "Tr0.="<<&tpc0<<
4626         "Tr1.="<<&tpc1<<
4627         "\n";
4628       //Bo:      if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4629       if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4630         v0mi->SetParamN(tpc0);
4631         v0mi->SetParamP(tpc1);
4632         v0mi->Update(primvertex);
4633       }
4634     }
4635   }
4636 }
4637 //------------------------------------------------------------------------
4638 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4639   //--------------------------------------------------------------------
4640   // Fill a look-up table with mean material
4641   //--------------------------------------------------------------------
4642
4643   Int_t n=1000;
4644   Double_t mparam[7];
4645   Double_t point1[3],point2[3];
4646   Double_t phi,cosphi,sinphi,z;
4647   // 0-5 layers, 6 pipe, 7-8 shields
4648   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4649   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4650
4651   Int_t ifirst=0,ilast=0;
4652   if(material.Contains("Pipe")) {
4653     ifirst=6; ilast=6;
4654   } else if(material.Contains("Shields")) {
4655     ifirst=7; ilast=8;
4656   } else if(material.Contains("Layers")) {
4657     ifirst=0; ilast=5;
4658   } else {
4659     Error("BuildMaterialLUT","Wrong layer name\n");
4660   }
4661
4662   for(Int_t imat=ifirst; imat<=ilast; imat++) {
4663     Double_t param[5]={0.,0.,0.,0.,0.};
4664     for (Int_t i=0; i<n; i++) {
4665       phi = 2.*TMath::Pi()*gRandom->Rndm();
4666       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4667       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4668       point1[0] = rmin[imat]*cosphi;
4669       point1[1] = rmin[imat]*sinphi;
4670       point1[2] = z;
4671       point2[0] = rmax[imat]*cosphi;
4672       point2[1] = rmax[imat]*sinphi;
4673       point2[2] = z;
4674       AliTracker::MeanMaterialBudget(point1,point2,mparam);
4675       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4676     }
4677     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4678     if(imat<=5) {
4679       fxOverX0Layer[imat] = param[1];
4680       fxTimesRhoLayer[imat] = param[0]*param[4];
4681     } else if(imat==6) {
4682       fxOverX0Pipe = param[1];
4683       fxTimesRhoPipe = param[0]*param[4];
4684     } else if(imat==7) {
4685       fxOverX0Shield[0] = param[1];
4686       fxTimesRhoShield[0] = param[0]*param[4];
4687     } else if(imat==8) {
4688       fxOverX0Shield[1] = param[1];
4689       fxTimesRhoShield[1] = param[0]*param[4];
4690     }
4691   }
4692   /*
4693   printf("%s\n",material.Data());
4694   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4695   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4696   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4697   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4698   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4699   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4700   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4701   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4702   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4703   */
4704   return;
4705 }
4706 //------------------------------------------------------------------------
4707 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4708                                               TString direction) {
4709   //-------------------------------------------------------------------
4710   // Propagate beyond beam pipe and correct for material
4711   // (material budget in different ways according to fUseTGeo value)
4712   //-------------------------------------------------------------------
4713
4714   // Define budget mode:
4715   // 0: material from AliITSRecoParam (hard coded)
4716   // 1: material from TGeo in one step (on the fly)
4717   // 2: material from lut
4718   // 3: material from TGeo in one step (same for all hypotheses)
4719   Int_t mode;
4720   switch(fUseTGeo) {
4721   case 0:
4722     mode=0;
4723     break;
4724   case 1:
4725     mode=1;
4726     break;
4727   case 2:
4728     mode=2;
4729     break;
4730   case 3:
4731     if(fTrackingPhase.Contains("Clusters2Tracks"))
4732       { mode=3; } else { mode=1; }
4733     break;
4734   case 4:
4735     if(fTrackingPhase.Contains("Clusters2Tracks"))
4736       { mode=3; } else { mode=2; }
4737     break;
4738   default:
4739     mode=0;
4740     break;
4741   }
4742   if(fTrackingPhase.Contains("Default")) mode=0;
4743
4744   Int_t index=fCurrentEsdTrack;
4745
4746   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4747   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4748   Double_t xToGo;
4749   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4750
4751   Double_t xOverX0,x0,lengthTimesMeanDensity;
4752   Bool_t anglecorr=kTRUE;
4753
4754   switch(mode) {
4755   case 0:
4756     xOverX0 = AliITSRecoParam::GetdPipe();
4757     x0 = AliITSRecoParam::GetX0Be();
4758     lengthTimesMeanDensity = xOverX0*x0;
4759     lengthTimesMeanDensity *= dir;
4760     if (!t->Propagate(xToGo)) return 0;
4761     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4762     break;
4763   case 1:
4764     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4765     break;
4766   case 2:
4767     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4768     xOverX0 = fxOverX0Pipe;
4769     lengthTimesMeanDensity = fxTimesRhoPipe;
4770     lengthTimesMeanDensity *= dir;
4771     if (!t->Propagate(xToGo)) return 0;
4772     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4773     break;
4774   case 3:
4775     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4776     if(fxOverX0PipeTrks[index]<0) {
4777       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4778       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4779                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4780       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4781       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4782       return 1;
4783     }
4784     xOverX0 = fxOverX0PipeTrks[index];
4785     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4786     lengthTimesMeanDensity *= dir;
4787     if (!t->Propagate(xToGo)) return 0;
4788     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4789     break;
4790   }
4791
4792   return 1;
4793 }
4794 //------------------------------------------------------------------------
4795 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4796                                                 TString shield,
4797                                                 TString direction) {
4798   //-------------------------------------------------------------------
4799   // Propagate beyond SPD or SDD shield and correct for material
4800   // (material budget in different ways according to fUseTGeo value)
4801   //-------------------------------------------------------------------
4802
4803   // Define budget mode:
4804   // 0: material from AliITSRecoParam (hard coded)
4805   // 1: material from TGeo in steps of X cm (on the fly)
4806   //    X = AliITSRecoParam::GetStepSizeTGeo()
4807   // 2: material from lut
4808   // 3: material from TGeo in one step (same for all hypotheses)
4809   Int_t mode;
4810   switch(fUseTGeo) {
4811   case 0:
4812     mode=0;
4813     break;
4814   case 1:
4815     mode=1;
4816     break;
4817   case 2:
4818     mode=2;
4819     break;
4820   case 3:
4821     if(fTrackingPhase.Contains("Clusters2Tracks"))
4822       { mode=3; } else { mode=1; }
4823     break;
4824   case 4:
4825     if(fTrackingPhase.Contains("Clusters2Tracks"))
4826       { mode=3; } else { mode=2; }
4827     break;
4828   default:
4829     mode=0;
4830     break;
4831   }
4832   if(fTrackingPhase.Contains("Default")) mode=0;
4833
4834   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4835   Double_t rToGo;
4836   Int_t    shieldindex=0;
4837   if (shield.Contains("SDD")) { // SDDouter
4838     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4839     shieldindex=1;
4840   } else if (shield.Contains("SPD")) {        // SPDouter
4841     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4842     shieldindex=0;
4843   } else {
4844     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4845     return 0;
4846   }
4847   Double_t xToGo;
4848   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4849
4850   Int_t index=2*fCurrentEsdTrack+shieldindex;
4851
4852   Double_t xOverX0,x0,lengthTimesMeanDensity;
4853   Bool_t anglecorr=kTRUE;
4854   Int_t nsteps=1;
4855
4856   switch(mode) {
4857   case 0:
4858     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4859     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4860     lengthTimesMeanDensity = xOverX0*x0;
4861     lengthTimesMeanDensity *= dir;
4862     if (!t->Propagate(xToGo)) return 0;
4863     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4864     break;
4865   case 1:
4866     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4867     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4868     break;
4869   case 2:
4870     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4871     xOverX0 = fxOverX0Shield[shieldindex];
4872     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4873     lengthTimesMeanDensity *= dir;
4874     if (!t->Propagate(xToGo)) return 0;
4875     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4876     break;
4877   case 3:
4878     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4879     if(fxOverX0ShieldTrks[index]<0) {
4880       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4881       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4882                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4883       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4884       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4885       return 1;
4886     }
4887     xOverX0 = fxOverX0ShieldTrks[index];
4888     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4889     lengthTimesMeanDensity *= dir;
4890     if (!t->Propagate(xToGo)) return 0;
4891     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4892     break;
4893   }
4894
4895   return 1;
4896 }
4897 //------------------------------------------------------------------------
4898 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4899                                                Int_t layerindex,
4900                                                Double_t oldGlobXYZ[3],
4901                                                TString direction) {
4902   //-------------------------------------------------------------------
4903   // Propagate beyond layer and correct for material
4904   // (material budget in different ways according to fUseTGeo value)
4905   //-------------------------------------------------------------------
4906
4907   // Define budget mode:
4908   // 0: material from AliITSRecoParam (hard coded)
4909   // 1: material from TGeo in stepsof X cm (on the fly)
4910   //    X = AliITSRecoParam::GetStepSizeTGeo()
4911   // 2: material from lut
4912   // 3: material from TGeo in one step (same for all hypotheses)
4913   Int_t mode;
4914   switch(fUseTGeo) {
4915   case 0:
4916     mode=0;
4917     break;
4918   case 1:
4919     mode=1;
4920     break;
4921   case 2:
4922     mode=2;
4923     break;
4924   case 3:
4925     if(fTrackingPhase.Contains("Clusters2Tracks"))
4926       { mode=3; } else { mode=1; }
4927     break;
4928   case 4:
4929     if(fTrackingPhase.Contains("Clusters2Tracks"))
4930       { mode=3; } else { mode=2; }
4931     break;
4932   default:
4933     mode=0;
4934     break;
4935   }
4936   if(fTrackingPhase.Contains("Default")) mode=0;
4937
4938   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4939
4940   Double_t r=fgLayers[layerindex].GetR();
4941   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4942
4943   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4944   Double_t xToGo;
4945   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4946
4947   Int_t index=6*fCurrentEsdTrack+layerindex;
4948
4949
4950   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4951   Double_t mparam[7];
4952   Bool_t anglecorr=kTRUE;
4953   Int_t nsteps=1;
4954   Double_t rOld,xOld;
4955
4956   Bool_t addTime = kFALSE;
4957   switch(mode) {
4958   case 0:
4959     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4960     lengthTimesMeanDensity = xOverX0*x0;
4961     // Bring the track beyond the material
4962     if (!t->Propagate(xToGo)) return 0;
4963     lengthTimesMeanDensity *= dir;
4964     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4965     break;
4966   case 1:
4967     rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4968     if (!t->GetLocalXat(rOld,xOld)) return 0;
4969     if (!t->Propagate(xOld)) return 0; // back before material (no correction)
4970     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4971     if (!t->PropagateToTGeo(xToGo,nsteps,addTime)) return 0; // cross the material and apply correction
4972     break;
4973   case 2:
4974     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4975     xOverX0 = fxOverX0Layer[layerindex];
4976     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4977     // Bring the track beyond the material
4978     if (!t->Propagate(xToGo)) return 0;
4979     lengthTimesMeanDensity *= dir;
4980     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4981     break;
4982   case 3:
4983     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4984     // Bring the track beyond the material
4985     if (!t->Propagate(xToGo)) return 0;
4986     Double_t globXYZ[3];
4987     if (!t->GetXYZ(globXYZ)) return 0;
4988     if (fxOverX0LayerTrks[index]<0) {
4989       AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4990       if(mparam[1]>900000) return 0;
4991       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4992                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4993       xOverX0=mparam[1]/angle;
4994       lengthTimesMeanDensity=mparam[0]*mparam[4]/angle;
4995       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0);
4996       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity);
4997     }
4998     xOverX0 = fxOverX0LayerTrks[index];
4999     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
5000     lengthTimesMeanDensity *= dir;
5001     if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
5002     break;
5003   }
5004
5005   return 1;
5006 }
5007 //------------------------------------------------------------------------
5008 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
5009   //-----------------------------------------------------------------
5010   // Initialize LUT for storing material for each prolonged track
5011   //-----------------------------------------------------------------
5012   fxOverX0PipeTrks = new Float_t[ntracks];
5013   fxTimesRhoPipeTrks = new Float_t[ntracks];
5014   fxOverX0ShieldTrks = new Float_t[ntracks*2];
5015   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
5016   fxOverX0LayerTrks = new Float_t[ntracks*6];
5017   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
5018
5019   for(Int_t i=0; i<ntracks; i++) {
5020     fxOverX0PipeTrks[i] = -1.;
5021     fxTimesRhoPipeTrks[i] = -1.;
5022   }
5023   for(Int_t j=0; j<ntracks*2; j++) {
5024     fxOverX0ShieldTrks[j] = -1.;
5025     fxTimesRhoShieldTrks[j] = -1.;
5026   }
5027   for(Int_t k=0; k<ntracks*6; k++) {
5028     fxOverX0LayerTrks[k] = -1.;
5029     fxTimesRhoLayerTrks[k] = -1.;
5030   }
5031
5032   fNtracks = ntracks;
5033
5034   return;
5035 }
5036 //------------------------------------------------------------------------
5037 void AliITStrackerMI::DeleteTrksMaterialLUT() {
5038   //-----------------------------------------------------------------
5039   // Delete LUT for storing material for each prolonged track
5040   //-----------------------------------------------------------------
5041   if(fxOverX0PipeTrks) {
5042     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
5043   }
5044   if(fxOverX0ShieldTrks) {
5045     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
5046   }
5047
5048   if(fxOverX0LayerTrks) {
5049     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
5050   }
5051   if(fxTimesRhoPipeTrks) {
5052     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
5053   }
5054   if(fxTimesRhoShieldTrks) {
5055     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
5056   }
5057   if(fxTimesRhoLayerTrks) {
5058     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
5059   }
5060   return;
5061 }
5062 //------------------------------------------------------------------------
5063 Int_t AliITStrackerMI::CheckSkipLayer(AliITStrackMI *track,
5064                                       Int_t ilayer,Int_t idet) const {
5065   //-----------------------------------------------------------------
5066   // This method is used to decide whether to allow a prolongation
5067   // without clusters, because we want to skip the layer.
5068   // In this case the return value is > 0:
5069   // return 1: the user requested to skip a layer
5070   // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
5071   //-----------------------------------------------------------------
5072
5073   if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
5074
5075   if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
5076     // check if track will cross SPD outer layer
5077     Double_t phiAtSPD2,zAtSPD2;
5078     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
5079       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
5080     }
5081   }
5082
5083   return 0;
5084 }
5085 //------------------------------------------------------------------------
5086 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
5087                                      Int_t ilayer,Int_t idet,
5088                                      Double_t dz,Double_t dy,
5089                                      Bool_t noClusters) const {
5090   //-----------------------------------------------------------------
5091   // This method is used to decide whether to allow a prolongation
5092   // without clusters, because there is a dead zone in the road.
5093   // In this case the return value is > 0:
5094   // return 1: dead zone at z=0,+-7cm in SPD
5095   // return 2: all road is "bad" (dead or noisy) from the OCDB
5096   // return 3: something "bad" (dead or noisy) from the OCDB
5097   //-----------------------------------------------------------------
5098
5099   // check dead zones at z=0,+-7cm in the SPD
5100   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
5101     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5102                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5103                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
5104     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5105                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5106                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
5107     for (Int_t i=0; i<3; i++)
5108       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
5109         AliDebug(2,Form("crack SPD %d",ilayer));
5110         return 1;
5111       }
5112   }
5113
5114   // check bad zones from OCDB
5115   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
5116
5117   if (idet<0) return 0;
5118
5119   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
5120
5121   Int_t detType=-1;
5122   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
5123   if (ilayer==0 || ilayer==1) {        // ----------  SPD
5124     detType = 0;
5125   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
5126     detType = 1;
5127     detSizeFactorX *= 2.;
5128   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
5129     detType = 2;
5130   }
5131   AliITSsegmentation *segm = (AliITSsegmentation*)fDetTypeRec->GetSegmentationModel(detType);
5132   if (detType==2) segm->SetLayer(ilayer+1);
5133   Float_t detSizeX = detSizeFactorX*segm->Dx();
5134   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
5135
5136   // check if the road overlaps with bad chips
5137   Float_t xloc,zloc;
5138   LocalModuleCoord(ilayer,idet,track,xloc,zloc);
5139   Float_t zlocmin = zloc-dz;
5140   Float_t zlocmax = zloc+dz;
5141   Float_t xlocmin = xloc-dy;
5142   Float_t xlocmax = xloc+dy;
5143   Int_t chipsInRoad[100];
5144
5145   // check if road goes out of detector
5146   Bool_t touchNeighbourDet=kFALSE;
5147   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
5148   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
5149   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5150   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5151   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));
5152
5153   // check if this detector is bad
5154   if (det.IsBad()) {
5155     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
5156     if(!touchNeighbourDet) {
5157       return 2; // all detectors in road are bad
5158     } else {
5159       return 3; // at least one is bad
5160     }
5161   }
5162
5163   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
5164   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
5165   if (!nChipsInRoad) return 0;
5166
5167   Bool_t anyBad=kFALSE,anyGood=kFALSE;
5168   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
5169     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
5170     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
5171     if (det.IsChipBad(chipsInRoad[iCh])) {
5172       anyBad=kTRUE;
5173     } else {
5174       anyGood=kTRUE;
5175     }
5176   }
5177
5178   if (!anyGood) {
5179     if(!touchNeighbourDet) {
5180       AliDebug(2,"all bad in road");
5181       return 2;  // all chips in road are bad
5182     } else {
5183       return 3; // at least a bad chip in road
5184     }
5185   }
5186
5187   if (anyBad) {
5188     AliDebug(2,"at least a bad in road");
5189     return 3; // at least a bad chip in road
5190   }
5191
5192
5193   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
5194       || ilayer==4 || ilayer==5     // SSD
5195       || !noClusters) return 0;
5196
5197   // There are no clusters in road: check if there is at least
5198   // a bad SPD pixel or SDD anode
5199
5200   Int_t idetInITS=idet;
5201   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
5202
5203   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
5204     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
5205     return 3;
5206   }
5207   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
5208
5209   return 0;
5210 }
5211 //------------------------------------------------------------------------
5212 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
5213                                        AliITStrackMI *track,
5214                                        Float_t &xloc,Float_t &zloc) const {
5215   //-----------------------------------------------------------------
5216   // Gives position of track in local module ref. frame
5217   //-----------------------------------------------------------------
5218
5219   xloc=0.;
5220   zloc=0.;
5221
5222   if(idet<0) return kFALSE;
5223
5224   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
5225
5226   Int_t lad = Int_t(idet/ndet) + 1;
5227
5228   Int_t det = idet - (lad-1)*ndet + 1;
5229
5230   Double_t xyzGlob[3],xyzLoc[3];
5231
5232   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
5233   // take into account the misalignment: xyz at real detector plane
5234   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
5235
5236   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
5237
5238   xloc = (Float_t)xyzLoc[0];
5239   zloc = (Float_t)xyzLoc[2];
5240
5241   return kTRUE;
5242 }
5243 //------------------------------------------------------------------------
5244 Bool_t AliITStrackerMI::IsOKForPlaneEff(AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
5245 //
5246 // Method to be optimized further:
5247 // Aim: decide whether a track can be used for PlaneEff evaluation
5248 //      the decision is taken based on the track quality at the layer under study
5249 //      no information on the clusters on this layer has to be used
5250 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
5251 //      the cut is done on number of sigmas from the boundaries
5252 //
5253 //  Input: Actual track, layer [0,5] under study
5254 //  Output: none
5255 //  Return: kTRUE if this is a good track
5256 //
5257 // it will apply a pre-selection to obtain good quality tracks.
5258 // Here also  you will have the possibility to put a control on the
5259 // impact point of the track on the basic block, in order to exclude border regions
5260 // this will be done by calling a proper method of the AliITSPlaneEff class.
5261 //
5262 // input: AliITStrackMI* track, ilayer= layer number [0,5]
5263 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
5264 //
5265   Int_t index[AliITSgeomTGeo::kNLayers];
5266   Int_t k;
5267   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
5268   //
5269   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
5270     index[k]=clusters[k];
5271   }
5272
5273   if(!fPlaneEff)
5274     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
5275   AliITSlayer &layer=fgLayers[ilayer];
5276   Double_t r=layer.GetR();
5277   AliITStrackMI tmp(*track);
5278
5279 // require a minimal number of cluster in other layers and eventually clusters in closest layers
5280   Int_t ncl=0;
5281   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
5282     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
5283                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
5284     if (tmp.GetClIndex(lay)>0) ncl++;
5285   }
5286   Bool_t nextout = kFALSE;
5287   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
5288   else nextout = ((tmp.GetClIndex(ilayer+1)>0)? kTRUE : kFALSE );
5289   Bool_t nextin = kFALSE;
5290   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
5291   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
5292   if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
5293      return kFALSE;
5294   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
5295   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
5296   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
5297  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
5298
5299 // detector number
5300   Double_t phi,z;
5301   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
5302   Int_t idet=layer.FindDetectorIndex(phi,z);
5303   if(idet<0) { AliInfo(Form("cannot find detector"));
5304     return kFALSE;}
5305
5306   // here check if it has good Chi Square.
5307
5308   //propagate to the intersection with the detector plane
5309   const AliITSdetector &det=layer.GetDetector(idet);
5310   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
5311
5312   Float_t locx; //
5313   Float_t locz; //
5314   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
5315   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5316   if(key>fPlaneEff->Nblock()) return kFALSE;
5317   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
5318   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
5319   //***************
5320   // DEFINITION OF SEARCH ROAD FOR accepting a track
5321   //
5322   //For the time being they are hard-wired, later on from AliITSRecoParam
5323   // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
5324   // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
5325   Double_t nsigz=4;
5326   Double_t nsigx=4;
5327   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
5328   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
5329                                                 // done for RecPoints
5330
5331   // exclude tracks at boundary between detectors
5332   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
5333   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
5334   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
5335   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
5336   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5337
5338   if ( (locx-dx < blockXmn+boundaryWidth) ||
5339        (locx+dx > blockXmx-boundaryWidth) ||
5340        (locz-dz < blockZmn+boundaryWidth) ||
5341        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5342   return kTRUE;
5343 }
5344 //------------------------------------------------------------------------
5345 void AliITStrackerMI::UseTrackForPlaneEff(AliITStrackMI* track, Int_t ilayer) {
5346 //
5347 // This Method has to be optimized! For the time-being it uses the same criteria
5348 // as those used in the search of extra clusters for overlapping modules.
5349 //
5350 // Method Purpose: estabilish whether a track has produced a recpoint or not
5351 //                 in the layer under study (For Plane efficiency)
5352 //
5353 // inputs: AliITStrackMI* track  (pointer to a usable track)
5354 // outputs: none
5355 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5356 //               traversed by this very track. In details:
5357 //               - if a cluster can be associated to the track then call
5358 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5359 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5360 //
5361   if(!fPlaneEff)
5362     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5363   AliITSlayer &layer=fgLayers[ilayer];
5364   Double_t r=layer.GetR();
5365   AliITStrackMI tmp(*track);
5366
5367 // detector number
5368   Double_t phi,z;
5369   if (!tmp.GetPhiZat(r,phi,z)) return;
5370   Int_t idet=layer.FindDetectorIndex(phi,z);
5371
5372   if(idet<0) { AliInfo(Form("cannot find detector"));
5373     return;}
5374
5375
5376 //propagate to the intersection with the detector plane
5377   const AliITSdetector &det=layer.GetDetector(idet);
5378   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5379
5380
5381 //***************
5382 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5383 //
5384   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5385                     TMath::Sqrt(tmp.GetSigmaZ2() +
5386                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5387                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5388                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5389   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5390                     TMath::Sqrt(tmp.GetSigmaY2() +
5391                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5392                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5393                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5394
5395 // road in global (rphi,z) [i.e. in tracking ref. system]
5396   Double_t zmin = tmp.GetZ() - dz;
5397   Double_t zmax = tmp.GetZ() + dz;
5398   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5399   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5400
5401 // select clusters in road
5402   layer.SelectClusters(zmin,zmax,ymin,ymax);
5403
5404 // Define criteria for track-cluster association
5405   Double_t msz = tmp.GetSigmaZ2() +
5406   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5407   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5408   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5409   Double_t msy = tmp.GetSigmaY2() +
5410   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5411   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5412   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5413   if (tmp.GetConstrain()) {
5414     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5415     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5416   }  else {
5417     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5418     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5419   }
5420   msz = 1./msz; // 1/RoadZ^2
5421   msy = 1./msy; // 1/RoadY^2
5422 //
5423
5424   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5425   Int_t idetc=-1;
5426   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5427   //Double_t  tolerance=0.2;
5428   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5429     idetc = cl->GetDetectorIndex();
5430     if(idet!=idetc) continue;
5431     //Int_t ilay = cl->GetLayer();
5432
5433     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5434     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5435
5436     Double_t chi2=tmp.GetPredictedChi2(cl);
5437     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5438   }*/
5439   Float_t locx; //
5440   Float_t locz; //
5441   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5442 //
5443   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5444   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5445   if(key>fPlaneEff->Nblock()) return;
5446   Bool_t found=kFALSE;
5447   //if (ci>=0) {
5448   Double_t chi2;
5449   while ((cl=layer.GetNextCluster(clidx))!=0) {
5450     idetc = cl->GetDetectorIndex();
5451     if(idet!=idetc) continue;
5452     // here real control to see whether the cluster can be associated to the track.
5453     // cluster not associated to track
5454     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5455          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
5456     // calculate track-clusters chi2
5457     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5458                                                // in particular, the error associated to the cluster
5459     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5460     // chi2 cut
5461     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5462     found=kTRUE;
5463     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5464    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5465    // track->SetExtraModule(ilayer,idetExtra);
5466   }
5467   if(!fPlaneEff->UpDatePlaneEff(found,key))
5468        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5469   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5470     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
5471     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5472     Int_t cltype[2]={-999,-999};
5473
5474     tr[0]=locx;
5475     tr[1]=locz;
5476     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
5477     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
5478
5479     if (found){
5480       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5481       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5482       cltype[0]=layer.GetCluster(ci)->GetNy();
5483       cltype[1]=layer.GetCluster(ci)->GetNz();
5484
5485      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5486      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
5487      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5488      // It is computed properly by calling the method
5489      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5490      // T
5491      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5492       //if (tmp.PropagateTo(x,0.,0.)) {
5493         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5494         AliCluster c(*layer.GetCluster(ci));
5495         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5496         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5497         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
5498         clu[2]=TMath::Sqrt(c.GetSigmaY2());
5499         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5500       //}
5501     }
5502     fPlaneEff->FillHistos(key,found,tr,clu,cltype);
5503   }
5504 return;
5505 }
5506