ITS/UPGRADE/AliITSUTrackerGlo.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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
  16 /*
  17   TO CHECK
  18   GetBz usage
  19
  20  */
  21
  22 //-------------------------------------------------------------------------
  23 //  Implementation of the ITS Upgrade tracker mother class.
  24 //-------------------------------------------------------------------------
  25 #include <TTree.h>
  26 #include <Riostream.h>
  27 #include <TMath.h>
  28
  29 #include "AliITSUTrackerGlo.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDtrack.h"
  32 #include "AliITSURecoDet.h"
  33 #include "AliITSURecoSens.h"
  34 #include "AliITSUReconstructor.h"
  35 #include "AliITSReconstructor.h"
  36 #include "AliITSUSeed.h"
  37 #include "AliITSUAux.h"
  38 #include "AliITSUClusterPix.h"
  39 using namespace AliITSUAux;
  40 using namespace TMath;
  41
  42
  43 ClassImp(AliITSUTrackerGlo)
  44 //_________________________________________________________________________
  45 AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
  46 :  fReconstructor(rec)
  47   ,fITS(0)
  48   ,fCurrMass(kPionMass)
  49   ,fSeedsLr(0)
  50   ,fSeedsPool("AliITSUSeed",0)
  51 {
  52   // Default constructor
  53   if (rec) Init(rec);
  54 }
  55
  56 //_________________________________________________________________________
  57 AliITSUTrackerGlo::~AliITSUTrackerGlo()
  58 {
  59  // Default destructor
  60  //
  61   delete fITS;
  62   delete[] fSeedsLr;
  63   //
  64 }
  65
  66 //_________________________________________________________________________
  67 void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
  68 {
  69   // init with external reconstructor
  70   //
  71   fITS = new AliITSURecoDet(rec->GetGeom(),"ITSURecoInterface");
  72   for (int ilr=fITS->GetNLayersActive();ilr--;) {
  73     fITS->GetLayerActive(ilr)->SetClusters(rec->GetClusters(ilr));
  74   }
  75   //
  76   fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
  77   int n = fITS->GetNLayersActive()+1;
  78   fSeedsLr = new TObjArray[n];
  79   //
  80
  81 }
  82
  83 //_________________________________________________________________________
  84 Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
  85 {
  86   //
  87   //
  88   AliITSUReconstructor::GetRecoParam()->Print();
  89
  90   fITS->ProcessClusters();
  91   // select ESD tracks to propagate
  92   int nTrESD = esdEv->GetNumberOfTracks();
  93   for (int itr=0;itr<nTrESD;itr++) {
  94     printf("Processing track %d\n",itr);
  95     AliESDtrack *esdTr = esdEv->GetTrack(itr);
  96     FindTrack(esdTr);
  97   }
  98
  99   return 0;
 100 }
 101
 102 //_________________________________________________________________________
 103 Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent * /*event*/)
 104 {
 105   //
 106   // To be implemented
 107   //
 108
 109  Info("PropagateBack","To be implemented");
 110   return 0;
 111 }
 112
 113 //_________________________________________________________________________
 114 Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent * /*event*/)
 115 {
 116   //
 117   // To be implemented
 118   //
 119
 120   Info("RefitInward","To be implemented");
 121   return 0;
 122 }
 123
 124 //_________________________________________________________________________
 125 Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
 126 {
 127   // read from tree (if pointer provided) or directly from the ITS reco interface
 128   //
 129   return fReconstructor->LoadClusters(treeRP);
 130 }
 131
 132 //_________________________________________________________________________
 133 void AliITSUTrackerGlo::UnloadClusters()
 134 {
 135   //
 136   // To be implemented
 137   //
 138
 139   Info("UnloadClusters","To be implemented");
 140 }
 141 //_________________________________________________________________________
 142 AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
 143 {
 144   //
 145   // To be implemented
 146   //
 147
 148   Info("GetCluster","To be implemented");
 149   return 0x0;
 150 }
 151
 152 //_________________________________________________________________________
 153 Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
 154 {
 155   // do we need to match this track to ITS?
 156   //
 157   static double bz = GetBz();
 158   if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
 159       esdTr->IsOn(AliESDtrack::kTPCout) ||
 160       esdTr->IsOn(AliESDtrack::kITSin)  ||
 161       esdTr->GetKinkIndex(0)>0) return kFALSE;
 162   //
 163   if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
 164   //
 165   float dtz[2];
 166   esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz);
 167   // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
 168   if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation())
 169        && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
 170            Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
 171   //
 172   return kTRUE;
 173 }
 174
 175 //_________________________________________________________________________
 176 void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr)
 177 {
 178   // find prolongaion candidates finding for single seed
 179   //
 180   if (!NeedToProlong(esdTr)) return;  // are we interested in this track?
 181   if (!InitSeed(esdTr))      return;  // initialize prolongations hypotheses tree
 182   //
 183   AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
 184   AliITSUSeed seedUC;  // copy of the seed from the upper layer
 185   AliITSUSeed seedT;   // transient seed between the seedUC and new prolongation hypothesis
 186   //
 187   TObjArray clArr; // container for transfer of clusters matching to seed
 188   //
 189   for (int ila=fITS->GetNLayersActive();ila--;) {
 190     int ilaUp = ila+1;                         // prolong seeds from layer above
 191     int nSeedsUp = GetNSeeds(ilaUp);
 192     printf("NSeeds %d\n",nSeedsUp);
 193     for (int isd=0;isd<nSeedsUp;isd++) {
 194       AliITSUSeed* seedU = GetSeed(ilaUp,isd);  // seed on prev.active layer to prolong
 195       seedUC = *seedU;
 196       seedUC.SetParent(seedU);
 197       // go till next active layer
 198       printf("Lr:%d Seed:%d\n",ila,isd);
 199       if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) {
 200         //
 201         printf("Tr failed\n");
 202         // Check if the seed satisfies to track definition
 203         if (NeedToKill(&seedUC,kTransportFailed)) KillSeed(ilaUp,isd);
 204         continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
 205       }
 206       AliITSURecoLayer* lrA = fITS->GetLayerActive(ila);
 207       if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
 208         if (NeedToKill(&seedUC,kRWCheckFailed)) KillSeed(ilaUp,isd);
 209         continue;
 210       }
 211       int nsens = lrA->FindSensors(&fTrImpData[kTrPhi0], hitSens);  // find detectors which may be hit by the track
 212       printf("Lr:%d Ns:%d\n",ila, nsens);
 213       //
 214       for (int isn=nsens;isn--;) {
 215         seedT = seedUC;
 216         AliITSURecoSens* sens = hitSens[isn];
 217         //
 218         if (!seedT.Propagate(sens->GetPhiTF(),sens->GetXTF(),GetBz())) continue; // propagation failed, seedT is intact
 219         int clID0 = sens->GetFirstClusterId();
 220         printf("Sn:%d Ncl:%d\n",isn,sens->GetNClusters());
 221         for (int icl=sens->GetNClusters();icl--;) {
 222           int res = CheckCluster(&seedT,ila,clID0+icl);
 223           //
 224           if (res==kStopSearchOnSensor) break;     // stop looking on this sensor
 225           if (res==kClusterNotMatching) continue;  // cluster does not match
 226           // cluster is matching and it was added to the hypotheses tree
 227         }
 228       }
 229       // cluster search is done. Do we need ta have a version of this seed skipping current layer
 230       seedT.SetLr(ila);
 231       if (!NeedToKill(&seedT,kMissingCluster)) AddProlongationHypothesis(NewSeedFromPool(&seedT) ,ila);
 232     }
 233   }
 234   //
 235   ResetSeedTree();
 236 }
 237
 238 //_________________________________________________________________________
 239 Bool_t AliITSUTrackerGlo::InitSeed(AliESDtrack *esdTr)
 240 {
 241   // init prolongaion candidates finding for single seed
 242   fCurrMass = esdTr->GetMass();
 243   if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
 244   //
 245   AliITSUSeed* seed = NewSeedFromPool();
 246   seed->SetLr(fITS->GetNLayersActive());   // fake layer
 247   seed->AliExternalTrackParam::operator=(*esdTr);
 248   seed->SetParent(esdTr);
 249   AddProlongationHypothesis(seed,fITS->GetNLayersActive());
 250   return kTRUE;
 251   // TO DO
 252 }
 253
 254 //_________________________________________________________________________
 255 void AliITSUTrackerGlo::ResetSeedTree()
 256 {
 257   // reset current hypotheses tree
 258   for (int i=fITS->GetNLayersActive()+1;i--;) fSeedsLr[fITS->GetNLayersActive()].Clear();
 259 }
 260
 261 //_________________________________________________________________________
 262 Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
 263 {
 264   // transport seed from layerFrom to the entrance of layerTo
 265   //
 266   const double kToler = 1e-6; // tolerance for layer on-surface check
 267   //
 268   int dir = lTo > lFrom ? 1:-1;
 269   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 270   Bool_t checkFirst = kTRUE;
 271   while(lFrom!=lTo) {
 272     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 273     if (lrFr) {
 274       Bool_t doLayer = kTRUE;
 275       double xToGo = dir>0 ? lrFr->GetRMax() : lrFr->GetRMin();
 276       if (checkFirst) { // do we need to track till the surface of the current layer ?
 277         checkFirst = kFALSE;
 278         if      (dir>0) { if (curR2-xToGo*xToGo>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
 279         else if (dir<0) { if (xToGo*xToGo-curR2>kToler) doLayer = kFALSE; } // on the surface or outside of the layer
 280       }
 281       if (doLayer) {
 282         if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 283         // go via layer to its boundary, applying material correction.
 284         if (!PropagateTrackTo(seed,xToGo,fCurrMass, lrFr->GetMaxStep(), kFALSE, -1, 0, kTRUE)) return kFALSE;
 285       }
 286     }
 287     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 288     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 289     //
 290     // go the entrance of the layer, assuming no materials in between
 291     double xToGo = dir>0 ? lrTo->GetRMin() : lrTo->GetRMax();
 292     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 293     if (!seed->PropagateTo(xToGo, GetBz())) return kFALSE; // RS: do we need BxByBz?
 294     lrFr = lrTo;
 295   }
 296   return kTRUE;
 297   //
 298 }
 299
 300 //_________________________________________________________________________
 301 Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
 302 {
 303   // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
 304   // as well as some aux info
 305   double bz = GetBz();
 306   AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
 307   seed->GetXYZ(&fTrImpData[kTrXIn]);    // lab position at the entrance from above
 308   //
 309   fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
 310   if (!seed->Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
 311   double dr  = lrA->GetDR();                              // approximate X dist at the inner radius
 312   if (!seed->GetXYZAt(seed->GetX()-dr, bz, fTrImpData + kTrXOut)) {
 313     // special case: track does not reach inner radius, might be tangential
 314     double r = seed->GetD(0,0,bz);
 315     double x;
 316     if (!seed->GetXatLabR(r,x,bz,-1)) {
 317       AliError(Form("This should not happen: r=%f",r));
 318       seed->Print();
 319       return kFALSE;
 320     }
 321     dr = Abs(seed->GetX() - x);
 322     if (!seed->GetXYZAt(x, bz, fTrImpData + kTrXOut)) {
 323       AliError(Form("This should not happen: x=%f",x));
 324       seed->Print();
 325       return kFALSE;
 326     }
 327   }
 328   //
 329   fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
 330   double sgy = seed->GetSigmaY2() + dr*dr*seed->GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
 331   double sgz = seed->GetSigmaZ2() + dr*dr*seed->GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
 332   sgy = Sqrt(sgy)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY();
 333   sgz = Sqrt(sgz)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ();
 334   fTrImpData[kTrPhi0] = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
 335   fTrImpData[kTrZ0]   = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
 336   fTrImpData[kTrDPhi] = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]) + sgy/lrA->GetR();
 337   fTrImpData[kTrDZ]   = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn])   + sgz;
 338   //
 339   return kTRUE;
 340 }
 341
 342 //_________________________________________________________________________
 343 AliITSUSeed* AliITSUTrackerGlo::NewSeedFromPool(const AliITSUSeed* src)
 344 {
 345   // create new seed, optionally copying from the source
 346   return src ?
 347     new(fSeedsPool[fSeedsPool.GetEntriesFast()]) AliITSUSeed(*src) :
 348     new(fSeedsPool[fSeedsPool.GetEntriesFast()]) AliITSUSeed();
 349 }
 350
 351 //_________________________________________________________________________
 352 Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
 353 {
 354   // Check if the cluster (in tracking frame!) is matching to track.
 355   // The track must be already propagated to sensor tracking frame.
 356   // Returns:  kStopSearchOnSensor if the search on given sensor should be stopped,
 357   //           kClusterMatching    if the cluster is matching
 358   //           kClusterMatching    otherwise
 359   //
 360   // The seed is already propagated to cluster
 361   const double kTolerX = 5e-4;
 362   AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
 363   //
 364   if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
 365     if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) return kStopSearchOnSensor; // propagation failed, seedT is intact
 366   }
 367   double dy = cl->GetY() - track->GetY();
 368   double dy2 = dy*dy;
 369   double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
 370     AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY(); // RS TOOPTIMIZE
 371   if (dy2>tol2) {                          // the clusters are sorted in Z(col) then in Y(row).
 372     if (dy>0) return kStopSearchOnSensor;  // No chance that other cluster of this sensor will match (all Y's will be even larger)
 373     else      return kClusterNotMatching;   // Other clusters may match
 374   }
 375   double dz2 = cl->GetZ()-track->GetZ();
 376   dz2 *= dz2;
 377   tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
 378     AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ(); // RS TOOPTIMIZE
 379   if (dz2>tol2) return kClusterNotMatching; // Other clusters may match
 380   //
 381   // check chi2
 382   Double_t p[2]={cl->GetY(), cl->GetZ()};
 383   Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
 384   double chi2 = track->GetPredictedChi2(p,cov);
 385   if (chi2>AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr)) return kClusterNotMatching;
 386   //
 387   track = NewSeedFromPool(track);  // input track will be reused, use its clone for updates
 388   if (!track->Update(p,cov)) return kClusterNotMatching;
 389   track->SetChi2Cl(chi2);
 390   track->SetLrClusterID(lr,clID);
 391   cl->IncreaseClusterUsage();
 392   //
 393   AddProlongationHypothesis(track,lr);
 394   //
 395   return kClusterMatching;
 396 }