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