ITS/AliITStrackerMI.cxx

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