HLT/ITS/tracking/AliITStrackerHLT.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: AliITStrackerHLT.cxx 32466 2009-05-20 07:51:56Z hristov $ */
  16
  17 //-------------------------------------------------------------------------
  18 //               Implementation of the ITS tracker class
  19 //    It reads AliITSRecPoint clusters and creates AliHLTITSTrack tracks
  20 //                   and fills with them the ESD
  21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  22 //          Current support and development:
  23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
  24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
  26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
  27 //-------------------------------------------------------------------------
  28
  29 //#include <TMatrixD.h>
  30 #include <TTree.h>
  31 #include <TDatabasePDG.h>
  32 #include <TString.h>
  33 #include <TRandom.h>
  34 #include <TTreeStream.h>
  35
  36
  37 #include "AliLog.h"
  38 #include "AliITSCalibrationSPD.h"
  39 #include "AliITSCalibrationSDD.h"
  40 #include "AliITSCalibrationSSD.h"
  41 #include "AliCDBEntry.h"
  42 #include "AliCDBManager.h"
  43 #include "AliAlignObj.h"
  44 #include "AliTrackPointArray.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliV0.h"
  49 #include "AliHelix.h"
  50 #include "AliITSChannelStatus.h"
  51 #include "AliITSRecPoint.h"
  52 #include "AliITSgeomTGeo.h"
  53 #include "AliITSReconstructor.h"
  54 #include "AliITSClusterParam.h"
  55 #include "AliITSsegmentation.h"
  56 #include "AliITSCalibration.h"
  57 #include "AliITSV0Finder.h"
  58 #include "AliITStrackerHLT.h"
  59 //#include "AliHLTTPCCATrackParam.h"
  60 //#include "AliHLTVertexer.h"
  61
  62
  63 ClassImp(AliITStrackerHLT)
  64
  65 Bool_t AliITStrackerHLT::TransportToX( AliExternalTrackParam *t, double x ) const
  66 {
  67   return t->PropagateTo( x, t->GetBz() );
  68 }
  69
  70 Bool_t AliITStrackerHLT::TransportToPhiX( AliExternalTrackParam *t, double phi, double x ) const
  71 {
  72   return t->Propagate( phi, x, t->GetBz() );
  73 }
  74
  75
  76 //------------------------------------------------------------------------
  77 Int_t AliITStrackerHLT::UpdateMI(AliHLTITSTrack* track, const AliITSRecPoint* cl,Double_t /*chi2*/,Int_t index) const
  78 {
  79   //
  80   // Update ITS track
  81   //
  82
  83   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
  84
  85   Int_t layer = (index & 0xf0000000) >> 28;
  86
  87   // Take into account the mis-alignment (bring track to cluster plane)
  88
  89   Double_t xTrOrig=track->GetX();
  90   if (!TransportToX( track, xTrOrig + cl->GetX() ) ) return 0;
  91
  92   Double_t err2Y, err2Z;
  93
  94   GetClusterErrors2( layer, cl, track, err2Y, err2Z );
  95
  96   Double_t p[2]={ cl->GetY(), cl->GetZ()};
  97   Double_t cov[3]={err2Y, 0., err2Z};
  98
  99   Int_t updated = 1;
 100   //if( layer!=2 && layer!=3 )
 101   updated = track->AliExternalTrackParam::Update(p,cov);
 102
 103   int n = track->GetNumberOfClusters();
 104   track->SetClusterIndex(n,index);
 105   track->SetNumberOfClusters(n+1);
 106
 107   return updated;
 108 }
 109
 110
 111 AliITStrackerHLT::AliITStrackerHLT()
 112   :AliTracker(),
 113    fRecoParam(0),
 114    fLayers(new AliHLTITSLayer[AliITSgeomTGeo::kNLayers]),
 115    fEsd(0),
 116    fUseTGeo(3),
 117    fxOverX0Pipe(-1.),
 118    fxTimesRhoPipe(-1.),
 119    fDebugStreamer(0),
 120    fITSChannelStatus(0),
 121    fTracks()
 122 {
 123   //Default constructor
 124   Int_t i;
 125   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
 126   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 127   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 128 }
 129 //------------------------------------------------------------------------
 130 AliITStrackerHLT::AliITStrackerHLT(const Char_t *geom)
 131 : AliTracker(),
 132   fRecoParam(0),
 133   fLayers(new AliHLTITSLayer[AliITSgeomTGeo::kNLayers]),
 134   fEsd(0),
 135   fUseTGeo(3),
 136   fxOverX0Pipe(-1.),
 137   fxTimesRhoPipe(-1.),
 138   fDebugStreamer(0),
 139   fITSChannelStatus(0),
 140   fTracks()
 141 {
 142   //--------------------------------------------------------------------
 143   //This is the AliITStrackerHLT constructor
 144   //--------------------------------------------------------------------
 145   if (geom) {
 146     AliWarning("\"geom\" is actually a dummy argument !");
 147   }
 148
 149   if(AliGeomManager::GetGeometry()==NULL){
 150     AliGeomManager::LoadGeometry();
 151   }
 152
 153   fRecoParam = AliITSReconstructor::GetRecoParam();
 154   if( !fRecoParam ){
 155     AliITSReconstructor *tmp = new AliITSReconstructor();
 156     tmp->Init();
 157     fRecoParam = AliITSRecoParam::GetLowFluxParam();
 158     tmp->AliReconstructor::SetRecoParam(fRecoParam);
 159   }
 160   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 161     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 162     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 163
 164     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 165     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 166     Double_t poff=TMath::ATan2(y,x);
 167     Double_t zoff=z;
 168     Double_t r=TMath::Sqrt(x*x + y*y);
 169
 170     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 171     r += TMath::Sqrt(x*x + y*y);
 172     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 173     r += TMath::Sqrt(x*x + y*y);
 174     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 175     r += TMath::Sqrt(x*x + y*y);
 176     r*=0.25;
 177
 178     new (fLayers+i-1) AliHLTITSLayer(r,poff,zoff,nlad,ndet);
 179
 180     for (Int_t j=1; j<nlad+1; j++) {
 181       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 182         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 183         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 184         m.Multiply(tm);
 185         Double_t txyz[3]={0.};
 186         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 187         m.LocalToMaster(txyz,xyz);
 188         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 189         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 190
 191         if (phi<0) phi+=TMath::TwoPi();
 192         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 193
 194         AliHLTITSDetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1);
 195         new(&det) AliHLTITSDetector(r,phi);
 196         // compute the real radius (with misalignment)
 197         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 198         mmisal.Multiply(tm);
 199         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 200         mmisal.LocalToMaster(txyz,xyz);
 201         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 202         det.SetRmisal(rmisal);
 203
 204       } // end loop on detectors
 205     } // end loop on ladders
 206   } // end loop on layers
 207
 208
 209   Double_t xyzVtx[]={ fRecoParam->GetXVdef(),
 210                       fRecoParam->GetYVdef(),
 211                       fRecoParam->GetZVdef()};
 212   Double_t ersVtx[]={ fRecoParam->GetSigmaXVdef(),
 213                       fRecoParam->GetSigmaYVdef(),
 214                       fRecoParam->GetSigmaZVdef()};
 215
 216   SetVertex(xyzVtx,ersVtx);
 217
 218   // store positions of centre of SPD modules (in z)
 219   Double_t tr[3];
 220   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 221   fSPDdetzcentre[0] = tr[2];
 222   AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 223   fSPDdetzcentre[1] = tr[2];
 224   AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 225   fSPDdetzcentre[2] = tr[2];
 226   AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 227   fSPDdetzcentre[3] = tr[2];
 228
 229   //fUseTGeo = fRecoParam->GetUseTGeoInTracker();
 230   //if(fRecoParam->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 231   //AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 232   //fUseTGeo = 3;
 233   //}
 234
 235   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 236   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 237
 238   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 239
 240 }
 241 //------------------------------------------------------------------------
 242 AliITStrackerHLT::AliITStrackerHLT(const AliITStrackerHLT &tracker)
 243 :AliTracker(tracker),
 244  fRecoParam( tracker.fRecoParam),
 245  fLayers(new AliHLTITSLayer[AliITSgeomTGeo::kNLayers]),
 246  fEsd(tracker.fEsd),
 247  fUseTGeo(tracker.fUseTGeo),
 248  fxOverX0Pipe(tracker.fxOverX0Pipe),
 249  fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 250  fDebugStreamer(tracker.fDebugStreamer),
 251  fITSChannelStatus(tracker.fITSChannelStatus),
 252  fTracks()
 253 {
 254   //Copy constructor
 255   Int_t i;
 256   for(i=0;i<4;i++) {
 257     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 258   }
 259   for(i=0;i<6;i++) {
 260     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 261     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 262   }
 263   for(i=0;i<2;i++) {
 264     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 265     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 266   }
 267 }
 268 //------------------------------------------------------------------------
 269 AliITStrackerHLT & AliITStrackerHLT::operator=(const AliITStrackerHLT &tracker){
 270   //Assignment operator
 271   this->~AliITStrackerHLT();
 272   new(this) AliITStrackerHLT(tracker);
 273   return *this;
 274 }
 275 //------------------------------------------------------------------------
 276 AliITStrackerHLT::~AliITStrackerHLT()
 277 {
 278   //
 279   //destructor
 280   //
 281   if (fDebugStreamer) {
 282     //fDebugStreamer->Close();
 283     delete fDebugStreamer;
 284   }
 285   if(fITSChannelStatus) delete fITSChannelStatus;
 286   delete [] fLayers;
 287 }
 288
 289
 290 //------------------------------------------------------------------------
 291 void AliITStrackerHLT::LoadClusters( std::vector<AliITSRecPoint> clusters)
 292 {
 293   //--------------------------------------------------------------------
 294   //This function loads ITS clusters
 295   //--------------------------------------------------------------------
 296
 297   //SignDeltas(clusters,GetZ());
 298   //std::cout<<"CA ITS: NClusters "<<clusters.size()<<std::endl;
 299   for( int i=0; i<AliITSgeomTGeo::GetNLayers(); i++ ){
 300      fLayers[i].ResetClusters(); //road defined by the cluster density
 301   }
 302
 303   for( unsigned int icl=0; icl<clusters.size(); icl++ ){
 304     //std::cout<<"CA ITS: icl "<<icl<<std::endl;
 305
 306     AliITSRecPoint &cl = clusters[icl];
 307     //std::cout<<"CA ITS: icl "<<icl<<": "
 308     //<<cl.GetX()<<" "<<cl.GetY()<<" "<<cl.GetZ()<<" "<<cl.GetDetectorIndex()<<" "<<cl.GetLayer()
 309     //<<" "<<cl.GetVolumeId()
 310     //<<std::endl;
 311
 312
 313     Int_t i=cl.GetLayer();
 314     //Int_t ndet=fLayers[i].GetNdetectors();
 315     //Int_t detector=cl.GetDetectorIndex();
 316     if (!cl.Misalign()) AliWarning("Can't misalign this cluster !"); //SG!!!
 317     fLayers[i].InsertCluster(new AliITSRecPoint(cl));
 318   }
 319
 320   for( int i=0; i<AliITSgeomTGeo::GetNLayers(); i++ ){
 321     if(0)for( int detector = 0; detector<fLayers[i].GetNdetectors(); detector++ ){ //SG!!!
 322       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 323       // zwindow cm from the dead zone
 324
 325       fRecoParam->GetAddVirtualClustersInDeadZone();
 326
 327       if (i<2 && fRecoParam->GetAddVirtualClustersInDeadZone()) {
 328         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*fRecoParam->GetXPassDeadZoneHits()) {
 329           Int_t lab[4]   = {0,0,0,detector};
 330           Int_t info[3]  = {0,0,i};
 331           Float_t q      = 0.; // this identifies virtual clusters
 332           Float_t hit[5] = {xdead,
 333                             0.,
 334                             fRecoParam->GetSigmaXDeadZoneHit2(),
 335                             fRecoParam->GetSigmaZDeadZoneHit2(),
 336                             q};
 337           Bool_t local   = kTRUE;
 338           Double_t zwindow = fRecoParam->GetZWindowDeadZone();
 339           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 340           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 341             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 342           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 343           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 344             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 345           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 346           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 347             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 348           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 349           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 350             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 351           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 352           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 353             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 354           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 355           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 356             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 357         }
 358       } // "virtual" clusters in SPD
 359     }
 360     fLayers[i].ResetRoad(); //road defined by the cluster density
 361     fLayers[i].SortClusters();
 362   }
 363 }
 364
 365
 366 //------------------------------------------------------------------------
 367 Int_t AliITStrackerHLT::LoadClusters(TTree *cTree) {
 368   //--------------------------------------------------------------------
 369   //This function loads ITS clusters
 370   //--------------------------------------------------------------------
 371   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 372   if (!branch) {
 373     Error("LoadClusters"," can't get the branch !\n");
 374     return 1;
 375   }
 376
 377   static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 378   branch->SetAddress(&clusters);
 379
 380   Int_t i=0,j=0,ndet=0;
 381   Int_t detector=0;
 382   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 383     ndet=fLayers[i].GetNdetectors();
 384     Int_t jmax = j + fLayers[i].GetNladders()*ndet;
 385     for (; j<jmax; j++) {
 386       if (!cTree->GetEvent(j)) continue;
 387       Int_t ncl=clusters->GetEntriesFast();
 388       SignDeltas(clusters,GetZ());
 389
 390       while (ncl--) {
 391         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 392         detector=c->GetDetectorIndex();
 393
 394         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 395
 396         fLayers[i].InsertCluster(new AliITSRecPoint(*c));
 397       }
 398       clusters->Delete();
 399       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 400       // zwindow cm from the dead zone
 401       if (i<2 && fRecoParam->GetAddVirtualClustersInDeadZone()) {
 402         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*fRecoParam->GetXPassDeadZoneHits()) {
 403           Int_t lab[4]   = {0,0,0,detector};
 404           Int_t info[3]  = {0,0,i};
 405           Float_t q      = 0.; // this identifies virtual clusters
 406           Float_t hit[5] = {xdead,
 407                             0.,
 408                             fRecoParam->GetSigmaXDeadZoneHit2(),
 409                             fRecoParam->GetSigmaZDeadZoneHit2(),
 410                             q};
 411           Bool_t local   = kTRUE;
 412           Double_t zwindow = fRecoParam->GetZWindowDeadZone();
 413           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 414           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 415             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 416           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 417           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 418             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 419           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 420           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 421             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 422           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 423           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 424             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 425           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 426           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 427             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 428           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 429           if (TMath::Abs(fLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 430             fLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 431         }
 432       } // "virtual" clusters in SPD
 433
 434     }
 435     //
 436     fLayers[i].ResetRoad(); //road defined by the cluster density
 437     fLayers[i].SortClusters();
 438   }
 439
 440   dummy.Clear();
 441
 442   return 0;
 443 }
 444 //------------------------------------------------------------------------
 445 void AliITStrackerHLT::UnloadClusters() {
 446   //--------------------------------------------------------------------
 447   //This function unloads ITS clusters
 448   //--------------------------------------------------------------------
 449   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayers[i].ResetClusters();
 450 }
 451
 452
 453
 454 //------------------------------------------------------------------------
 455 Int_t AliITStrackerHLT::CorrectForTPCtoITSDeadZoneMaterial(AliHLTITSTrack *t) {
 456   //--------------------------------------------------------------------
 457   // Correction for the material between the TPC and the ITS
 458   //--------------------------------------------------------------------
 459   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 460       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 461       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 462       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 463   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 464       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 465       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 466       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 467   } else {
 468     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
 469     return 0;
 470   }
 471
 472   return 1;
 473 }
 474
 475 #include "TStopwatch.h"
 476
 477 void AliITStrackerHLT::Reconstruct( std::vector<AliExternalTrackParam> tracksTPC )
 478 {
 479   //--------------------------------------------------------------------
 480   // This functions reconstructs ITS tracks
 481   // The clusters must be already loaded !
 482   //--------------------------------------------------------------------
 483
 484   Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 485   fTracks.clear();
 486
 487   for( unsigned int itr=0; itr<tracksTPC.size(); itr++ ){
 488
 489     AliHLTITSTrack tMI( tracksTPC[itr] );
 490     AliHLTITSTrack *t = &tMI;
 491     t->SetTPCtrackId( itr );
 492     t->SetMass(pimass);
 493     t->SetExpQ(0);
 494
 495     if (!CorrectForTPCtoITSDeadZoneMaterial(t))  continue;
 496
 497     //Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 498
 499     FollowProlongationTree(t);
 500     int nclu=0;
 501     for(Int_t i=0; i<6; i++) {
 502       if( t->GetClusterIndex(i)>=0 ) nclu++;
 503     }
 504     //cout<<"N assigned ITS clusters = "<<nclu<<std::endl;
 505     if( nclu>0 ){
 506       t->SetLabel(-1);//tpcLabel);
 507       t->SetFakeRatio(1.);
 508       CookLabel(t,0.); //For comparison only
 509       //cout<<"label = "<<t->GetLabel()<<" / "<<tpcLabel<<endl;
 510       TransportToX(t, 0 );
 511       //cout<<"\n fill track : parameters at "<<t->GetX()<<": "<< TMath::Sqrt(TMath::Abs(t->GetSigmaY2()))<<" "<< TMath::Sqrt(TMath::Abs(t->GetSigmaY2()))<<endl;
 512         //t->Print();
 513       fTracks.push_back( *t );
 514     }
 515   }
 516
 517   //AliHLTVertexer vertexer;
 518   //vertexer.SetESD( event );
 519   //vertexer.FindPrimaryVertex();
 520   //vertexer.FindV0s();
 521 }
 522
 523
 524
 525 //------------------------------------------------------------------------
 526 Int_t AliITStrackerHLT::Clusters2Tracks(AliESDEvent *event) {
 527   //--------------------------------------------------------------------
 528   // This functions reconstructs ITS tracks
 529   // The clusters must be already loaded !
 530   //--------------------------------------------------------------------
 531
 532   TStopwatch timer;
 533
 534   fEsd = event;         // store pointer to the esd
 535   std::vector<AliExternalTrackParam> tracksTPC;
 536
 537   for( int itr=0; itr<event->GetNumberOfTracks(); itr++ ){
 538
 539     AliESDtrack *esdTrack = event->GetTrack(itr);
 540
 541     if ((esdTrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 542     if (esdTrack->GetStatus()&AliESDtrack::kTPCout) continue;
 543     if (esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
 544     if (esdTrack->GetKinkIndex(0)>0) continue;   //kink daughter
 545
 546     AliHLTITSTrack t(*esdTrack);
 547     t.SetTPCtrackId( itr );
 548     tracksTPC.push_back( t );
 549   }
 550
 551   Reconstruct( tracksTPC );
 552
 553   for( unsigned int itr=0; itr<fTracks.size(); itr++ ){
 554     AliHLTITSTrack &t = fTracks[itr];
 555     UpdateESDtrack(event->GetTrack(t.TPCtrackId()), &t, AliESDtrack::kITSin);
 556   }
 557
 558   //AliHLTVertexer vertexer;
 559   //vertexer.SetESD( event );
 560   //vertexer.FindPrimaryVertex();
 561   //vertexer.FindV0s();
 562
 563   timer.Stop();
 564   static double totalTime = 0;
 565   static int nEvnts = 0;
 566   totalTime+=timer.CpuTime();
 567   nEvnts++;
 568   std::cout<<"\n\n ITS tracker time = "<<totalTime/nEvnts<<" [s/ev]  for "<<nEvnts<<" events\n\n "<<std::endl;
 569   return 0;
 570 }
 571
 572
 573 //------------------------------------------------------------------------
 574 Int_t AliITStrackerHLT::PropagateBack(AliESDEvent * /*event*/) {
 575   return 0;
 576 }
 577
 578 //------------------------------------------------------------------------
 579 Int_t AliITStrackerHLT::RefitInward(AliESDEvent * /*event*/ ) {
 580   return 0;
 581 }
 582 //------------------------------------------------------------------------
 583 AliCluster *AliITStrackerHLT::GetCluster(Int_t index) const {
 584   //--------------------------------------------------------------------
 585   //       Return pointer to a given cluster
 586   //--------------------------------------------------------------------
 587   Int_t l=(index & 0xf0000000) >> 28;
 588   Int_t c=(index & 0x0fffffff) >> 00;
 589   return fLayers[l].GetCluster(c);
 590 }
 591 //------------------------------------------------------------------------
 592 Bool_t AliITStrackerHLT::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 593   //--------------------------------------------------------------------
 594   // Get track space point with index i
 595   //--------------------------------------------------------------------
 596
 597   Int_t l=(index & 0xf0000000) >> 28;
 598   Int_t c=(index & 0x0fffffff) >> 00;
 599   AliITSRecPoint *cl = fLayers[l].GetCluster(c);
 600   Int_t idet = cl->GetDetectorIndex();
 601
 602   Float_t xyz[3];
 603   Float_t cov[6];
 604   cl->GetGlobalXYZ(xyz);
 605   cl->GetGlobalCov(cov);
 606   p.SetXYZ(xyz, cov);
 607   p.SetCharge(cl->GetQ());
 608   p.SetDriftTime(cl->GetDriftTime());
 609   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 610   switch (l) {
 611   case 0:
 612     iLayer = AliGeomManager::kSPD1;
 613     break;
 614   case 1:
 615     iLayer = AliGeomManager::kSPD2;
 616     break;
 617   case 2:
 618     iLayer = AliGeomManager::kSDD1;
 619     break;
 620   case 3:
 621     iLayer = AliGeomManager::kSDD2;
 622     break;
 623   case 4:
 624     iLayer = AliGeomManager::kSSD1;
 625     break;
 626   case 5:
 627     iLayer = AliGeomManager::kSSD2;
 628     break;
 629   default:
 630     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 631     break;
 632   };
 633   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 634   p.SetVolumeID((UShort_t)volid);
 635   return kTRUE;
 636 }
 637 //------------------------------------------------------------------------
 638 Bool_t AliITStrackerHLT::GetTrackPointTrackingError(Int_t index,
 639                         AliTrackPoint& p, const AliESDtrack *t) {
 640   //--------------------------------------------------------------------
 641   // Get track space point with index i
 642   // (assign error estimated during the tracking)
 643   //--------------------------------------------------------------------
 644
 645   Int_t l=(index & 0xf0000000) >> 28;
 646   Int_t c=(index & 0x0fffffff) >> 00;
 647   const AliITSRecPoint *cl = fLayers[l].GetCluster(c);
 648   Int_t idet = cl->GetDetectorIndex();
 649
 650   const AliHLTITSDetector &det=fLayers[l].GetDetector(idet);
 651
 652   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 653   Float_t detxy[2];
 654   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 655   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 656   Double_t alpha = t->GetAlpha();
 657   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 658   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
 659   phi += alpha-det.GetPhi();
 660   Float_t tgphi = TMath::Tan(phi);
 661
 662   Float_t tgl = t->GetTgl(); // tgl about const along track
 663   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 664
 665   Float_t errlocalx,errlocalz;
 666   Bool_t addMisalErr=kFALSE;
 667
 668   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
 669
 670   Float_t xyz[3];
 671   Float_t cov[6];
 672   cl->GetGlobalXYZ(xyz);
 673   //  cl->GetGlobalCov(cov);
 674   Float_t pos[3] = {0.,0.,0.};
 675   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
 676   tmpcl.GetGlobalCov(cov);
 677
 678   p.SetXYZ(xyz, cov);
 679   p.SetCharge(cl->GetQ());
 680   p.SetDriftTime(cl->GetDriftTime());
 681
 682   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 683   switch (l) {
 684   case 0:
 685     iLayer = AliGeomManager::kSPD1;
 686     break;
 687   case 1:
 688     iLayer = AliGeomManager::kSPD2;
 689     break;
 690   case 2:
 691     iLayer = AliGeomManager::kSDD1;
 692     break;
 693   case 3:
 694     iLayer = AliGeomManager::kSDD2;
 695     break;
 696   case 4:
 697     iLayer = AliGeomManager::kSSD1;
 698     break;
 699   case 5:
 700     iLayer = AliGeomManager::kSSD2;
 701     break;
 702   default:
 703     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 704     break;
 705   };
 706   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 707
 708   p.SetVolumeID((UShort_t)volid);
 709   return kTRUE;
 710 }
 711
 712
 713 //------------------------------------------------------------------------
 714 void AliITStrackerHLT::FollowProlongationTree(AliHLTITSTrack * track )
 715 {
 716     //cout<<endl;
 717   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
 718     //cout<<"\nLayer "<<ilayer<<endl;
 719
 720     AliHLTITSLayer &layer=fLayers[ilayer];
 721
 722     //cout<<" shield material.. "<<endl;
 723
 724     // material between SSD and SDD, SDD and SPD
 725     if (ilayer==3 && !CorrectForShieldMaterial(track,"SDD","inward")) continue;
 726     if (ilayer==1 && !CorrectForShieldMaterial(track,"SPD","inward")) continue;
 727
 728     int idet;
 729
 730     {
 731       // propagate to the layer radius
 732
 733       Double_t r = layer.GetR(), phi,z;
 734       Double_t xToGo;
 735       //cout<<" propagate to layer R= "<<r<<" .."<<endl;
 736       if( !track->GetLocalXat(r,xToGo) ) continue;
 737       if( !TransportToX(track, xToGo) ) continue;
 738
 739       // detector number
 740
 741       if (!track->GetPhiZat(r,phi,z)) continue;
 742       idet=layer.FindDetectorIndex(phi,z);
 743       //cout<<" detector number = "<<idet<<endl;
 744    }
 745
 746
 747     //cout<<" correct for the layer material .. "<<endl;
 748
 749     // correct for the layer material
 750     {
 751       Double_t trackGlobXYZ1[3];
 752       if (!track->GetXYZ(trackGlobXYZ1)) continue;
 753       CorrectForLayerMaterial(track,ilayer,trackGlobXYZ1,"inward");
 754     }
 755
 756     // track outside layer acceptance in z
 757
 758     if( idet<0 ) continue;
 759
 760     // propagate to the intersection with the detector plane
 761
 762     //cout<<" propagate to the intersection with the detector .. "<<endl;
 763
 764     const AliHLTITSDetector &det=layer.GetDetector( idet );
 765     if (!TransportToPhiX( track, det.GetPhi(), det.GetR() ) ) continue;
 766
 767     // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
 768
 769     // road in global (rphi,z) [i.e. in tracking ref. system]
 770
 771     Double_t zmin,zmax,ymin,ymax;
 772
 773     //cout<<" ComputeRoad .. "<<endl;
 774     if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
 775
 776     //cout<<" Road: y["<<ymin<<","<<ymax<<"], z["<<zmin<<","<<zmax<<"] "<<endl;
 777
 778     // select clusters in road
 779
 780     //cout<<" SelectClusters .. "<<endl;
 781     layer.SelectClusters(zmin,zmax,ymin,ymax);
 782
 783     // Define criteria for track-cluster association
 784
 785     Double_t msz = track->GetSigmaZ2() +
 786       fRecoParam->GetNSigmaZLayerForRoadZ()*
 787       fRecoParam->GetNSigmaZLayerForRoadZ()*
 788       fRecoParam->GetSigmaZ2(ilayer);
 789
 790     Double_t msy = track->GetSigmaY2() +
 791       fRecoParam->GetNSigmaYLayerForRoadY()*
 792       fRecoParam->GetNSigmaYLayerForRoadY()*
 793       fRecoParam->GetSigmaY2(ilayer);
 794
 795      msz *= fRecoParam->GetNSigma2RoadZNonC();
 796      msy *= fRecoParam->GetNSigma2RoadYNonC();
 797
 798      msz = 1./msz; // 1/RoadZ^2
 799      msy = 1./msy; // 1/RoadY^2
 800
 801      //
 802      // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
 803      //
 804
 805      const AliITSRecPoint *cl=0;
 806      Int_t clidx=-1;
 807      Double_t chi2trkcl=fRecoParam->GetMaxChi2(); // init with big value
 808      Bool_t deadzoneSPD=kFALSE;
 809
 810      // check if the road contains a dead zone
 811      Bool_t noClusters = !layer.GetNextCluster(clidx,kTRUE);
 812
 813      Double_t dz=0.5*(zmax-zmin);
 814      Double_t dy=0.5*(ymax-ymin);
 815
 816      Int_t dead = 0;//CheckDeadZone(track,ilayer,idet,dz,dy,noClusters);
 817
 818      // create a prolongation without clusters (check also if there are no clusters in the road)
 819
 820      if (dead==1) { // dead zone at z=0,+-7cm in SPD
 821        deadzoneSPD=kTRUE;
 822      }
 823
 824      clidx=-1;
 825
 826      //cout<<" loop over clusters in the road .. "<<endl;
 827
 828      // loop over clusters in the road
 829      const AliITSRecPoint *bestCluster=0;
 830      double bestChi2 = 1.e10;
 831      AliHLTITSTrack bestTrack( *track );
 832      int bestIdx = -1;
 833      while ((cl=layer.GetNextCluster(clidx))!=0) {
 834        //cout<<" cluster: "<<cl->GetX()<<" "<<cl->GetY()<<" "<<cl->GetZ()<<endl;
 835        AliHLTITSTrack t(*track);
 836        if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
 837
 838        Int_t idetc=cl->GetDetectorIndex();
 839
 840        //cout<<" cluster detector: "<<idetc<<endl;
 841
 842        if ( idet !=idetc ) { // new cluster's detector
 843          const AliHLTITSDetector &detc=layer.GetDetector(idetc);
 844          if (!TransportToPhiX( track, detc.GetPhi(),detc.GetR()) ) continue;
 845          t = *track;
 846          idet = idetc;
 847        }
 848
 849        // take into account misalignment (bring track to real detector plane)
 850
 851        if (!TransportToX( &t, t.GetX() + cl->GetX() ) ) continue;
 852        double chi2 = ( (t.GetZ()-cl->GetZ())*(t.GetZ()-cl->GetZ())*msz +
 853                        (t.GetY()-cl->GetY())*(t.GetY()-cl->GetY())*msy   );
 854        //cout<<" chi2="<<chi2<<endl;
 855        if ( chi2 < bestChi2 ){
 856          bestChi2 = chi2;
 857          bestCluster = cl;
 858          bestTrack = t;
 859          bestIdx = clidx;
 860          continue;
 861        }
 862      }
 863
 864      //cout<<" best chi2= "<<bestChi2<<endl;
 865
 866      if( bestCluster && bestChi2 <=1. ){
 867
 868        //cout<<" cluster found "<<endl;
 869        *track = bestTrack;
 870
 871
 872       // calculate track-clusters chi2
 873        chi2trkcl = GetPredictedChi2MI(track,bestCluster,ilayer);
 874        //cout<<" track-clusters chi2 = "<<chi2trkcl<<endl;
 875        //cout<<" max chi2 = "<<fRecoParam->GetMaxChi2s(ilayer)<<endl;
 876
 877        // chi2 cut
 878        if (chi2trkcl < fRecoParam->GetMaxChi2s(ilayer)) {
 879          if (bestCluster->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
 880          //cout<<"set index.."<<endl;
 881          //cout<<"set index ok"<<endl;
 882          if (bestCluster->GetQ()!=0) { // real cluster
 883            //cout<<" UpdateMI ... "<<endl;
 884            if (!UpdateMI(track,bestCluster,chi2trkcl,(ilayer<<28)+bestIdx)) {
 885              continue;
 886            }
 887          }
 888        }
 889      }
 890      //cout<<" goto next layer "<<endl;
 891
 892   }
 893 }
 894
 895
 896 //------------------------------------------------------------------------
 897 AliHLTITSLayer & AliITStrackerHLT::GetLayer(Int_t layer) const
 898 {
 899   //--------------------------------------------------------------------
 900   //
 901   //
 902   return fLayers[layer];
 903 }
 904
 905
 906
 907 //------------------------------------------------------------------------
 908 void AliITStrackerHLT::CookLabel(AliHLTITSTrack *track,Float_t wrong) const
 909 {
 910   // get MC label for the track
 911
 912   Int_t mcLabel = -1;
 913
 914   vector<int> labels;
 915   Int_t nClusters = track->GetNumberOfClusters();
 916
 917   for (Int_t i=0; i<nClusters; i++){
 918     Int_t cindex = track->GetClusterIndex(i);
 919     //Int_t l=(cindex & 0xf0000000) >> 28;
 920     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
 921     if ( cl->GetLabel(0) >= 0 ) labels.push_back(cl->GetLabel(0)) ;
 922     if ( cl->GetLabel(1) >= 0 ) labels.push_back(cl->GetLabel(1)) ;
 923     if ( cl->GetLabel(2) >= 0 ) labels.push_back(cl->GetLabel(2)) ;
 924   }
 925   std::sort( labels.begin(), labels.end() );
 926   labels.push_back( -1 ); // put -1 to the end
 927   int labelMax = -1, labelCur = -1, nLabelsMax = 0, nLabelsCurr = 0;
 928   for ( unsigned int iLab = 0; iLab < labels.size(); iLab++ ) {
 929     if ( labels[iLab] != labelCur ) {
 930       if ( labelCur >= 0 && nLabelsMax< nLabelsCurr ) {
 931         nLabelsMax = nLabelsCurr;
 932         labelMax = labelCur;
 933       }
 934       labelCur = labels[iLab];
 935       nLabelsCurr = 0;
 936     }
 937     nLabelsCurr++;
 938   }
 939
 940   if( labelMax>=0 && nLabelsMax <= wrong * nClusters ) labelMax = -labelMax;
 941
 942   mcLabel = labelMax;
 943
 944   track->SetLabel( mcLabel );
 945 }
 946
 947
 948
 949 void AliITStrackerHLT::GetClusterErrors2( Int_t layer, const AliITSRecPoint *cluster, AliHLTITSTrack* track, double &err2Y, double &err2Z ) const
 950 {
 951   Float_t erry,errz;
 952   Float_t theta = track->GetTgl();
 953   Float_t phi   = track->GetSnp();
 954   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
 955   Float_t expQ = TMath::Max((float)track->GetExpQ(),30.f);
 956   AliITSClusterParam::GetError(layer,cluster,theta,phi,expQ,erry,errz);
 957   err2Y = erry*erry;
 958   err2Z = errz*errz;
 959 }
 960
 961
 962 //------------------------------------------------------------------------
 963 Double_t AliITStrackerHLT::GetPredictedChi2MI(AliHLTITSTrack* track, const AliITSRecPoint *cluster,Int_t layer)
 964 {
 965   //
 966   // Compute predicted chi2
 967   //
 968
 969   AliHLTITSTrack t(*track);
 970   if (!t.Propagate( t.GetX()+cluster->GetX())) return 1000.;
 971
 972   Double_t  err2Y,err2Z;
 973   GetClusterErrors2( layer, cluster, &t, err2Y, err2Z );
 974
 975   Double_t chi2 = t.GetPredictedChi2(cluster->GetY(),cluster->GetZ(),err2Y,err2Z);
 976
 977   Float_t ny,nz;
 978   Float_t theta = track->GetTgl();
 979   Float_t phi   = track->GetSnp();
 980   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
 981   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
 982   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
 983   if (delta>1){
 984     chi2+=0.5*TMath::Min(delta/2,2.);
 985     chi2+=2.*cluster->GetDeltaProbability();
 986   }
 987   return chi2;
 988 }
 989
 990
 991
 992
 993 //------------------------------------------------------------------------
 994 void AliITStrackerHLT::SignDeltas(const TObjArray *clusterArray, Float_t vz)
 995 {
 996   //
 997   // Clusters from delta electrons?
 998   //
 999   Int_t entries = clusterArray->GetEntriesFast();
1000   if (entries<4) return;
1001   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
1002   Int_t layer = cluster->GetLayer();
1003   if (layer>1) return;
1004   Int_t index[10000];
1005   Int_t ncandidates=0;
1006   Float_t r = (layer>0)? 7:4;
1007   //
1008   for (Int_t i=0;i<entries;i++){
1009     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
1010     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
1011     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
1012     index[ncandidates] = i;  //candidate to belong to delta electron track
1013     ncandidates++;
1014     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
1015       cl0->SetDeltaProbability(1);
1016     }
1017   }
1018   //
1019   //
1020   //
1021   for (Int_t i=0;i<ncandidates;i++){
1022     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
1023     if (cl0->GetDeltaProbability()>0.8) continue;
1024     //
1025     Int_t ncl = 0;
1026     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
1027     sumy=sumz=sumy2=sumyz=sumw=0.0;
1028     for (Int_t j=0;j<ncandidates;j++){
1029       if (i==j) continue;
1030       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
1031       //
1032       Float_t dz = cl0->GetZ()-cl1->GetZ();
1033       Float_t dy = cl0->GetY()-cl1->GetY();
1034       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
1035         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
1036         y[ncl] = cl1->GetY();
1037         z[ncl] = cl1->GetZ();
1038         sumy+= y[ncl]*weight;
1039         sumz+= z[ncl]*weight;
1040         sumy2+=y[ncl]*y[ncl]*weight;
1041         sumyz+=y[ncl]*z[ncl]*weight;
1042         sumw+=weight;
1043         ncl++;
1044       }
1045     }
1046     if (ncl<4) continue;
1047     Float_t det = sumw*sumy2  - sumy*sumy;
1048     Float_t delta=1000;
1049     if (TMath::Abs(det)>0.01){
1050       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
1051       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
1052       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
1053     }
1054     else{
1055       Float_t z0  = sumyz/sumy;
1056       delta = TMath::Abs(cl0->GetZ()-z0);
1057     }
1058     if ( delta<0.05) {
1059       cl0->SetDeltaProbability(1-20.*delta);
1060     }
1061   }
1062 }
1063 //------------------------------------------------------------------------
1064 void AliITStrackerHLT::UpdateESDtrack(AliESDtrack *tESD, AliHLTITSTrack* track, ULong_t flags) const
1065 {
1066   //
1067   // Update ESD track
1068   //
1069   tESD->UpdateTrackParams(track,flags);
1070   AliHLTITSTrack * oldtrack = (AliHLTITSTrack*)(tESD->GetITStrack());
1071   if (oldtrack) delete oldtrack;
1072   tESD->SetITStrack(new AliHLTITSTrack(*track));
1073 }
1074
1075
1076
1077
1078 //------------------------------------------------------------------------
1079 void AliITStrackerHLT::BuildMaterialLUT(TString material) {
1080   //--------------------------------------------------------------------
1081   // Fill a look-up table with mean material
1082   //--------------------------------------------------------------------
1083
1084   Int_t n=1000;
1085   Double_t mparam[7];
1086   Double_t point1[3],point2[3];
1087   Double_t phi,cosphi,sinphi,z;
1088   // 0-5 layers, 6 pipe, 7-8 shields
1089   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
1090   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
1091
1092   Int_t ifirst=0,ilast=0;
1093   if(material.Contains("Pipe")) {
1094     ifirst=6; ilast=6;
1095   } else if(material.Contains("Shields")) {
1096     ifirst=7; ilast=8;
1097   } else if(material.Contains("Layers")) {
1098     ifirst=0; ilast=5;
1099   } else {
1100     Error("BuildMaterialLUT","Wrong layer name\n");
1101   }
1102
1103   for(Int_t imat=ifirst; imat<=ilast; imat++) {
1104     Double_t param[5]={0.,0.,0.,0.,0.};
1105     for (Int_t i=0; i<n; i++) {
1106       phi = 2.*TMath::Pi()*gRandom->Rndm();
1107       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
1108       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
1109       point1[0] = rmin[imat]*cosphi;
1110       point1[1] = rmin[imat]*sinphi;
1111       point1[2] = z;
1112       point2[0] = rmax[imat]*cosphi;
1113       point2[1] = rmax[imat]*sinphi;
1114       point2[2] = z;
1115       AliTracker::MeanMaterialBudget(point1,point2,mparam);
1116       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
1117     }
1118     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
1119     if(imat<=5) {
1120       fxOverX0Layer[imat] = param[1];
1121       fxTimesRhoLayer[imat] = param[0]*param[4];
1122     } else if(imat==6) {
1123       fxOverX0Pipe = param[1];
1124       fxTimesRhoPipe = param[0]*param[4];
1125     } else if(imat==7) {
1126       fxOverX0Shield[0] = param[1];
1127       fxTimesRhoShield[0] = param[0]*param[4];
1128     } else if(imat==8) {
1129       fxOverX0Shield[1] = param[1];
1130       fxTimesRhoShield[1] = param[0]*param[4];
1131     }
1132   }
1133   /*
1134   printf("%s\n",material.Data());
1135   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
1136   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
1137   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
1138   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
1139   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
1140   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
1141   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
1142   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
1143   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
1144   */
1145   return;
1146 }
1147
1148 //------------------------------------------------------------------------
1149 Int_t AliITStrackerHLT::CorrectForPipeMaterial(AliHLTITSTrack *t,
1150                                               TString direction) {
1151   //-------------------------------------------------------------------
1152   // Propagate beyond beam pipe and correct for material
1153   // (material budget in different ways according to fUseTGeo value)
1154   // Add time if going outward (PropagateTo or PropagateToTGeo)
1155   //-------------------------------------------------------------------
1156
1157   // Define budget mode:
1158   // 0: material from AliITSRecoParam (hard coded)
1159   // 1: material from TGeo in one step (on the fly)
1160   // 2: material from lut
1161   // 3: material from TGeo in one step (same for all hypotheses)
1162   Int_t mode;
1163   switch(fUseTGeo) {
1164   case 0:
1165     mode=0;
1166     break;
1167   case 1:
1168     mode=1;
1169     break;
1170   case 2:
1171     mode=2;
1172     break;
1173   case 3:
1174     mode=3;
1175     break;
1176   case 4:
1177     mode=3;
1178     break;
1179   default:
1180     mode=0;
1181     break;
1182   }
1183
1184   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1185   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
1186   Double_t xToGo;
1187   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1188
1189   Double_t xOverX0,x0,lengthTimesMeanDensity;
1190
1191   switch(mode) {
1192   case 0:
1193     xOverX0 = AliITSRecoParam::GetdPipe();
1194     x0 = AliITSRecoParam::GetX0Be();
1195     lengthTimesMeanDensity = xOverX0*x0;
1196     lengthTimesMeanDensity *= dir;
1197     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1198     break;
1199   case 1:
1200     if (!t->PropagateToTGeo(xToGo,1)) return 0;
1201     break;
1202   case 2:
1203     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
1204     xOverX0 = fxOverX0Pipe;
1205     lengthTimesMeanDensity = fxTimesRhoPipe;
1206     lengthTimesMeanDensity *= dir;
1207     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1208     break;
1209   case 3:
1210     double xOverX0PipeTrks, xTimesRhoPipeTrks;
1211     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1212     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1213                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
1214     xOverX0PipeTrks = TMath::Abs(xOverX0)/angle;
1215     xTimesRhoPipeTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1216     xOverX0 = xOverX0PipeTrks;
1217     lengthTimesMeanDensity = xTimesRhoPipeTrks;
1218     lengthTimesMeanDensity *= dir;
1219     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1220     break;
1221   }
1222
1223   return 1;
1224 }
1225 //------------------------------------------------------------------------
1226 Int_t AliITStrackerHLT::CorrectForShieldMaterial(AliHLTITSTrack *t,
1227                                                 TString shield,
1228                                                 TString direction) {
1229   //-------------------------------------------------------------------
1230   // Propagate beyond SPD or SDD shield and correct for material
1231   // (material budget in different ways according to fUseTGeo value)
1232   // Add time if going outward (PropagateTo or PropagateToTGeo)
1233   //-------------------------------------------------------------------
1234
1235   // Define budget mode:
1236   // 0: material from AliITSRecoParam (hard coded)
1237   // 1: material from TGeo in steps of X cm (on the fly)
1238   //    X = AliITSRecoParam::GetStepSizeTGeo()
1239   // 2: material from lut
1240   // 3: material from TGeo in one step (same for all hypotheses)
1241   Int_t mode;
1242   switch(fUseTGeo) {
1243   case 0:
1244     mode=0;
1245     break;
1246   case 1:
1247     mode=1;
1248     break;
1249   case 2:
1250     mode=2;
1251     break;
1252   case 3:
1253     mode=3;
1254     break;
1255   case 4:
1256     mode=3;
1257     break;
1258   default:
1259     mode=0;
1260     break;
1261   }
1262
1263
1264   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1265   Double_t rToGo;
1266   Int_t    shieldindex=0;
1267   if (shield.Contains("SDD")) { // SDDouter
1268     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
1269     shieldindex=1;
1270   } else if (shield.Contains("SPD")) {        // SPDouter
1271     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
1272     shieldindex=0;
1273   } else {
1274     Error("CorrectForShieldMaterial"," Wrong shield name\n");
1275     return 0;
1276   }
1277   Double_t xToGo;
1278   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1279
1280   Double_t xOverX0,x0,lengthTimesMeanDensity;
1281   Int_t nsteps=1;
1282
1283   switch(mode) {
1284   case 0:
1285     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
1286     x0 = AliITSRecoParam::GetX0shield(shieldindex);
1287     lengthTimesMeanDensity = xOverX0*x0;
1288     lengthTimesMeanDensity *= dir;
1289     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1290     break;
1291   case 1:
1292     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
1293     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1294     break;
1295   case 2:
1296     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
1297     xOverX0 = fxOverX0Shield[shieldindex];
1298     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
1299     lengthTimesMeanDensity *= dir;
1300     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1301     break;
1302   case 3:
1303     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1304     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1305                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
1306     double xOverX0ShieldTrks = TMath::Abs(xOverX0)/angle;
1307     double xTimesRhoShieldTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1308     xOverX0 = xOverX0ShieldTrks;
1309     lengthTimesMeanDensity = xTimesRhoShieldTrks;
1310     lengthTimesMeanDensity *= dir;
1311     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1312     break;
1313   }
1314
1315   return 1;
1316 }
1317 //------------------------------------------------------------------------
1318 Int_t AliITStrackerHLT::CorrectForLayerMaterial(AliHLTITSTrack *t,
1319                                                Int_t layerindex,
1320                                                Double_t oldGlobXYZ[3],
1321                                                TString direction) {
1322   //-------------------------------------------------------------------
1323   // Propagate beyond layer and correct for material
1324   // (material budget in different ways according to fUseTGeo value)
1325   // Add time if going outward (PropagateTo or PropagateToTGeo)
1326   //-------------------------------------------------------------------
1327
1328   // Define budget mode:
1329   // 0: material from AliITSRecoParam (hard coded)
1330   // 1: material from TGeo in stepsof X cm (on the fly)
1331   //    X = AliITSRecoParam::GetStepSizeTGeo()
1332   // 2: material from lut
1333   // 3: material from TGeo in one step (same for all hypotheses)
1334   Int_t mode;
1335   switch(fUseTGeo) {
1336   case 0:
1337     mode=0;
1338     break;
1339   case 1:
1340     mode=1;
1341     break;
1342   case 2:
1343     mode=2;
1344     break;
1345   case 3:
1346     mode=3;
1347     break;
1348   case 4:
1349     mode=3;
1350     break;
1351   default:
1352     mode=0;
1353     break;
1354   }
1355
1356
1357   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1358
1359   Double_t r=fLayers[layerindex].GetR();
1360   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
1361
1362   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
1363   Double_t xToGo;
1364   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1365
1366
1367   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
1368   Int_t nsteps=1;
1369
1370   // back before material (no correction)
1371   Double_t rOld,xOld;
1372   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1373   if (!t->GetLocalXat(rOld,xOld)) return 0;
1374   if (!TransportToX( t, xOld)) return 0;
1375
1376   switch(mode) {
1377   case 0:
1378     xOverX0 = fLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
1379     lengthTimesMeanDensity = xOverX0*x0;
1380     lengthTimesMeanDensity *= dir;
1381     // Bring the track beyond the material
1382     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1383     break;
1384   case 1:
1385     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
1386     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1387     break;
1388   case 2:
1389     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
1390     xOverX0 = fxOverX0Layer[layerindex];
1391     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1392     lengthTimesMeanDensity *= dir;
1393     // Bring the track beyond the material
1394     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1395     break;
1396   case 3:
1397     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
1398     if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
1399     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1400                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
1401     double xOverX0LayerTrks = TMath::Abs(xOverX0)/angle;
1402     double xTimesRhoLayerTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1403     xOverX0 = xOverX0LayerTrks;
1404     lengthTimesMeanDensity = xTimesRhoLayerTrks;
1405     lengthTimesMeanDensity *= dir;
1406     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1407     break;
1408   }
1409
1410
1411   return 1;
1412 }
1413
1414
1415 //------------------------------------------------------------------------
1416 Bool_t AliITStrackerHLT::LocalModuleCoord(Int_t ilayer,Int_t idet,
1417                                        const AliHLTITSTrack *track,
1418                                        Float_t &xloc,Float_t &zloc) const {
1419   //-----------------------------------------------------------------
1420   // Gives position of track in local module ref. frame
1421   //-----------------------------------------------------------------
1422
1423   xloc=0.;
1424   zloc=0.;
1425
1426   if(idet<0) return kFALSE;
1427
1428   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
1429
1430   Int_t lad = Int_t(idet/ndet) + 1;
1431
1432   Int_t det = idet - (lad-1)*ndet + 1;
1433
1434   Double_t xyzGlob[3],xyzLoc[3];
1435
1436   AliHLTITSDetector &detector = fLayers[ilayer].GetDetector(idet);
1437   // take into account the misalignment: xyz at real detector plane
1438   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
1439
1440   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
1441
1442   xloc = (Float_t)xyzLoc[0];
1443   zloc = (Float_t)xyzLoc[2];
1444
1445   return kTRUE;
1446 }
1447
1448 //------------------------------------------------------------------------
1449 Bool_t AliITStrackerHLT::ComputeRoad(AliHLTITSTrack* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1450   //--------------------------------------------------------------------
1451   // This function computes the rectangular road for this track
1452   //--------------------------------------------------------------------
1453
1454
1455   AliHLTITSDetector &det = fLayers[ilayer].GetDetector(idet);
1456   // take into account the misalignment: propagate track to misaligned detector plane
1457   if (!TransportToPhiX( track, det.GetPhi(),det.GetRmisal() ) ) return kFALSE;
1458
1459   Double_t dz=fRecoParam->GetNSigmaRoadZ()*
1460                     TMath::Sqrt(track->GetSigmaZ2() +
1461                                 fRecoParam->GetNSigmaZLayerForRoadZ()*
1462                                 fRecoParam->GetNSigmaZLayerForRoadZ()*
1463                                 fRecoParam->GetSigmaZ2(ilayer));
1464   Double_t dy=fRecoParam->GetNSigmaRoadY()*
1465                     TMath::Sqrt(track->GetSigmaY2() +
1466                                 fRecoParam->GetNSigmaYLayerForRoadY()*
1467                                 fRecoParam->GetNSigmaYLayerForRoadY()*
1468                                 fRecoParam->GetSigmaY2(ilayer));
1469
1470   // track at boundary between detectors, enlarge road
1471   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1472   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1473        (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1474        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1475        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1476     Float_t tgl = TMath::Abs(track->GetTgl());
1477     if (tgl > 1.) tgl=1.;
1478     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1479     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1480     Float_t snp = TMath::Abs(track->GetSnp());
1481     if (snp > fRecoParam->GetMaxSnp()) return kFALSE;
1482     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1483   } // boundary
1484
1485   // add to the road a term (up to 2-3 mm) to deal with misalignments
1486   dy = TMath::Sqrt(dy*dy + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1487   dz = TMath::Sqrt(dz*dz + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1488
1489   Double_t r = fLayers[ilayer].GetR();
1490   zmin = track->GetZ() - dz;
1491   zmax = track->GetZ() + dz;
1492   ymin = track->GetY() + r*det.GetPhi() - dy;
1493   ymax = track->GetY() + r*det.GetPhi() + dy;
1494
1495   // bring track back to idead detector plane
1496   if (!TransportToPhiX( track, det.GetPhi(),det.GetR())) return kFALSE;
1497
1498   return kTRUE;
1499 }
1500