]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITSUTrackerGlo.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added implementation for phi at the TPC ouer radius for AOD
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUTrackerGlo.cxx
CommitLineData
32d38de2 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>
76390254 28#include <TFile.h>
32d38de2 29
30#include "AliITSUTrackerGlo.h"
31#include "AliESDEvent.h"
32#include "AliESDtrack.h"
33#include "AliITSURecoDet.h"
34#include "AliITSURecoSens.h"
35#include "AliITSUReconstructor.h"
36#include "AliITSReconstructor.h"
37#include "AliITSUSeed.h"
38#include "AliITSUAux.h"
c61e50c3 39#include "AliITSUClusterPix.h"
32d38de2 40using namespace AliITSUAux;
41using namespace TMath;
42
76390254 43#define _FILL_CONTROL_HISTOS_
32d38de2 44
76390254 45#ifdef _FILL_CONTROL_HISTOS_
46//
47TObjArray* cHistoArr;
48enum {kHResY=0,kHResYP=10,kHResZ=20,kHResZP=30,kHChi2Cl=40};
49enum {kHistosPhase=100};
50//
51#endif
e2d2481c 52
c8d1f258 53//----------------- tmp stuff -----------------
c8d1f258 54
32d38de2 55ClassImp(AliITSUTrackerGlo)
e2d2481c 56
57const Double_t AliITSUTrackerGlo::fgkToler = 1e-6;// tolerance for layer on-surface check
58
59
32d38de2 60//_________________________________________________________________________
61AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
62: fReconstructor(rec)
63 ,fITS(0)
0091e9f0 64 ,fCurrESDtrack(0)
0ddbf657 65 ,fCurrESDtrMClb(kDummyLabel)
32d38de2 66 ,fCurrMass(kPionMass)
0ddbf657 67 ,fCountProlongationTrials(0)
68 ,fCountITSin(0)
69 ,fCountITSout(0)
70 ,fCountITSrefit(0)
716ccba7 71 ,fHypStore(100)
38997c3c 72 ,fNBranchesAdded(0)
73 ,fNCandidatesAdded(0)
716ccba7 74 ,fCurrHyp(0)
c61e50c3 75 ,fSeedsPool("AliITSUSeed",0)
b8b59e05 76 ,fFreeSeedsID(0)
77 ,fNFreeSeeds(0)
78 ,fLastSeedID(0)
42c3d4bd 79 ,fDefTrackConds(0)
80 ,fCurrTrackCond(0)
c51c3124 81 ,fTrackPhase(-1)
8b16dbae 82 ,fClInfo(0)
32d38de2 83{
84 // Default constructor
85 if (rec) Init(rec);
86}
87
88//_________________________________________________________________________
89AliITSUTrackerGlo::~AliITSUTrackerGlo()
90{
91 // Default destructor
92 //
8b16dbae 93 delete[] fClInfo;
32d38de2 94 //
76390254 95#ifdef _FILL_CONTROL_HISTOS_
96 if (cHistoArr) {
97 TFile* ctrOut = TFile::Open("itsuControlHistos.root","recreate");
98 ctrOut->cd();
99 AliInfo("Storing control histos");
100 cHistoArr->Print();
101 // ctrOut->WriteObject(cHistoArr,"controlH","kSingleKey");
102 cHistoArr->Write();
103 ctrOut->Close();
104 delete ctrOut;
105 cHistoArr = 0;
106 }
107#endif
108 //
32d38de2 109}
110
111//_________________________________________________________________________
112void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
113{
114 // init with external reconstructor
c61e50c3 115 //
3d4dc3e2 116 fITS = rec->GetITSInterface();
8b16dbae 117 int nLr = fITS->GetNLayersActive();
118 fClInfo = new Int_t[nLr<<1];
c61e50c3 119 //
120 fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
b8b59e05 121 fFreeSeedsID.Set(1000);
c61e50c3 122 //
32d38de2 123}
124
42c3d4bd 125//_________________________________________________________________________
126void AliITSUTrackerGlo::CreateDefaultTrackCond()
127{
128 // creates default tracking conditions to be used when recoparam does not provide them
129 int nLr = fITS->GetNLayersActive();
130 fClInfo = new Int_t[nLr<<1];
131 //
132 AliITSUTrackCond* cond = new AliITSUTrackCond();
133 //
134 cond->SetNLayers(fITS->GetNLayersActive());
135 cond->AddNewCondition(nLr);
136 cond->AddGroupPattern( 0xffff ); // require all layers hit
137 //
138 fDefTrackConds.AddLast(cond);
139 //
140 AliInfo("Created conditions: ");
141 for (int i=0;i<fDefTrackConds.GetEntriesFast();i++) fDefTrackConds[i]->Print();
142 //
143}
144
145
32d38de2 146//_________________________________________________________________________
147Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
148{
149 //
76390254 150 SetTrackingPhase(kClus2Tracks);
151 //
152#ifdef _FILL_CONTROL_HISTOS_
153 //
154 if (!cHistoArr) { // create control histos
155 const int kNResDef=7;
156 const double kResDef[kNResDef]={0.05,0.05,0.3, 0.05,1,0.5,1.5};
157 cHistoArr = new TObjArray();
158 cHistoArr->SetOwner(kTRUE);
159 const double ptMax=10;
160 const double plMax=10;
161 const double chiMax=100;
162 const int nptbins=50;
163 const int nresbins=400;
164 const int nplbins=50;
165 const int nchbins=200;
166 int nblr = fITS->GetNLayersActive();
167 TString ttl;
168 for (int stp=0;stp<kNTrackingPhases;stp++) {
169 for (int ilr=0;ilr<nblr;ilr++) {
170 int hoffs = stp*kHistosPhase + ilr;
171 double mxdf = ilr>=kNResDef ? kResDef[kNResDef-1] : kResDef[ilr];
172 ttl = Form("S%d_residY%d",stp,ilr);
173 TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
174 cHistoArr->AddAtAndExpand(hdy,hoffs + kHResY);
175 hdy->SetDirectory(0);
176 //
177 ttl = Form("S%d_residYPull%d",stp,ilr);
178 TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
179 cHistoArr->AddAtAndExpand(hdyp,hoffs + kHResYP);
180 hdyp->SetDirectory(0);
181 //
182 ttl = Form("S%d_residZ%d",stp,ilr);
183 TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
184 cHistoArr->AddAtAndExpand(hdz,hoffs + kHResZ);
185 hdz->SetDirectory(0);
186 //
187 ttl = Form("S%d_residZPull%d",stp,ilr);
188 TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
189 hdzp->SetDirectory(0);
190 cHistoArr->AddAtAndExpand(hdzp,hoffs + kHResZP);
191 //
192 ttl = Form("S%d_chi2Cl%d",stp,ilr);
193 TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
194 hchi->SetDirectory(0);
195 cHistoArr->AddAtAndExpand(hchi,hoffs + kHChi2Cl);
196 }
197 }
198 }
199 //
200#endif
201
42c3d4bd 202 TObjArray *trackConds = 0;
32d38de2 203 //
0ddbf657 204 fCountProlongationTrials = 0;
205 fCountITSin = 0;
206 fCountITSout = 0;
207 fCountITSrefit = 0;
208 //
b8b59e05 209 ResetSeedsPool();
716ccba7 210 int nTrESD = esdEv->GetNumberOfTracks();
c51c3124 211 AliInfo(Form("Will try to find prolongations for %d tracks",nTrESD));
42c3d4bd 212 int nTrackCond = AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();
213 if (nTrackCond<1) {
214 if (!fDefTrackConds.GetEntriesFast()) {
215 AliInfo("No tracking conditions found in recoparams, creating default one requesting all layers hit");
216 CreateDefaultTrackCond();
217 }
218 trackConds = &fDefTrackConds;
219 nTrackCond = trackConds->GetEntriesFast();
220 }
221 else trackConds = AliITSUReconstructor::GetRecoParam()->GetTrackingConditions();
222 //
15b02d69 223 static Bool_t first = kTRUE;
224 if (first) {
225 AliITSUReconstructor::GetRecoParam()->Print();
226 first = kFALSE;
227 }
228 fHypStore.Delete();
229 if (fHypStore.GetSize()<nTrESD) fHypStore.Expand(nTrESD+100);
230 //
231 fITS->ProcessClusters();
232 //
42c3d4bd 233 for (int icnd=0;icnd<nTrackCond;icnd++) {
234 fCurrTrackCond = (AliITSUTrackCond*)trackConds->UncheckedAt(icnd);
235 // select ESD tracks to propagate
236 for (int itr=0;itr<nTrESD;itr++) {
0ddbf657 237 fCurrESDtrack = esdEv->GetTrack(itr);
238 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
239 //
240 if (!NeedToProlong(fCurrESDtrack)) continue; // are we interested in this track?
38997c3c 241 AliDebug(+1,Form("Processing track %d | M=%.3f Pt=%.3f | MCLabel: %d",itr,fCurrESDtrack->GetMass(kTRUE),fCurrESDtrack->Pt(),fCurrESDtrMClb));//RS
0ddbf657 242 FindTrack(fCurrESDtrack, itr);
42c3d4bd 243 }
244 //
38997c3c 245 if (AliDebugLevelClass()>+2) {
15b02d69 246 AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
247 fHypStore.Print();
248 }
42c3d4bd 249 FinalizeHypotheses();
250 }
716ccba7 251 //
0ddbf657 252 AliInfo(Form("%d ITSin for %d tried TPC seeds out of %d ESD tracks\n",fCountITSin,fCountProlongationTrials,nTrESD));
32d38de2 253 return 0;
254}
255
256//_________________________________________________________________________
c51c3124 257Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent *esdEv)
32d38de2 258{
259 //
c51c3124 260 // Do outward fits in ITS
261 //
e7d83d38 262 SetTrackingPhase(kPropBack);
c51c3124 263 int nTrESD = esdEv->GetNumberOfTracks();
15b02d69 264 AliDebug(1,Form("Will propagate back %d tracks",nTrESD));
c51c3124 265 //
266 double bz0 = GetBz();
267 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
268 AliITSUTrackHyp dummyTr,*currTr=0;
269 const double kWatchStep=10.; // for larger steps watch arc vs segment difference
270 Double_t times[AliPID::kSPECIES];
271 //
272 //
273 for (int itr=0;itr<nTrESD;itr++) {
0ddbf657 274 fCurrESDtrack = esdEv->GetTrack(itr);
275 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
c51c3124 276 // Start time integral and add distance from current position to vertex
0ddbf657 277 if (fCurrESDtrack->IsOn(AliESDtrack::kITSout)) continue; // already done
c51c3124 278 //
0ddbf657 279 fCurrESDtrack->GetXYZ(xyzTrk);
c51c3124 280 Double_t dst = 0.; // set initial track lenght, tof
281 {
282 double dxs = xyzTrk[0] - xyzVtx[0];
283 double dys = xyzTrk[1] - xyzVtx[1];
284 double dzs = xyzTrk[2] - xyzVtx[2];
285 // RS: for large segment steps d use approximation of cicrular arc s by
286 // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
287 // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
288 // Hence s^2/d^2 = (1+1/6 p^2)^2
289 dst = dxs*dxs + dys*dys;
290 if (dst > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
0ddbf657 291 double crv = Abs(fCurrESDtrack->GetC(bz0));
c51c3124 292 double fcarc = 1.+crv*crv*dst/6.;
293 dst *= fcarc*fcarc;
294 }
295 dst += dzs*dzs;
296 dst = Sqrt(dst);
297 }
298 //
0ddbf657 299 fCurrESDtrack->SetStatus(AliESDtrack::kTIME);
c51c3124 300 //
0ddbf657 301 if (!fCurrESDtrack->IsOn(AliESDtrack::kITSin)) { // case of tracks w/o ITS prolongation: just set the time integration
302 dummyTr.AliExternalTrackParam::operator=(*fCurrESDtrack);
c51c3124 303 dummyTr.StartTimeIntegral();
e2d2481c 304 dummyTr.AddTimeStep(dst);
c51c3124 305 dummyTr.GetIntegratedTimes(times);
0ddbf657 306 fCurrESDtrack->SetIntegratedTimes(times);
307 fCurrESDtrack->SetIntegratedLength(dummyTr.GetIntegratedLength());
c51c3124 308 continue;
309 }
310 //
311 currTr = GetTrackHyp(itr);
312 currTr->StartTimeIntegral();
313 currTr->AddTimeStep(dst);
15b02d69 314 // printf("Before resetCov: "); currTr->AliExternalTrackParam::Print();
e2d2481c 315 currTr->ResetCovariance(10000);
68a0f687 316 if (RefitTrack(currTr,fITS->GetRMax())) { // propagate to exit from the ITS/TPC screen
317 UpdateESDTrack(currTr,AliESDtrack::kITSout);
0ddbf657 318 fCountITSout++;
68a0f687 319 }
e7d83d38 320 else {
0ddbf657 321 AliDebug(2,Form("Refit Failed for track %d | ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
322 //currTr->AliExternalTrackParam::Print();
323 //currTr->GetWinner()->Print();
e7d83d38 324 }
c51c3124 325 //
326 }
32d38de2 327 //
0ddbf657 328 AliInfo(Form("%d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
329 fCountITSout,fCountITSin,fCountProlongationTrials,nTrESD));
330 //
32d38de2 331 return 0;
332}
333
334//_________________________________________________________________________
e7d83d38 335Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent *esdEv)
32d38de2 336{
337 //
e7d83d38 338 // refit the tracks inward, using current cov.matrix
339 //
340 SetTrackingPhase(kRefitInw);
341 Int_t nTrESD = esdEv->GetNumberOfTracks();
15b02d69 342 AliDebug(1,Form("Will refit inward %d tracks",nTrESD));
e7d83d38 343 AliITSUTrackHyp *currTr=0;
344 //
345 for (int itr=0;itr<nTrESD;itr++) {
0ddbf657 346 fCurrESDtrack = esdEv->GetTrack(itr);
347 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
e7d83d38 348 // Start time integral and add distance from current position to vertex
0ddbf657 349 UInt_t trStat = fCurrESDtrack->GetStatus();
e7d83d38 350 if ( !(trStat & AliESDtrack::kITSout) ) continue;
351 if ( trStat & AliESDtrack::kITSrefit ) continue; // already done
352 if ( (trStat & AliESDtrack::kTPCout) && !(trStat & AliESDtrack::kTPCrefit) ) continue;
353 //
354 currTr = GetTrackHyp(itr);
0ddbf657 355 currTr->AliExternalTrackParam::operator=(*fCurrESDtrack); // fetch current ESDtrack kinematics
e7d83d38 356 if (RefitTrack(currTr,fITS->GetRMin())) { // propagate up to inside radius of the beam pipe
357 UpdateESDTrack(currTr,AliESDtrack::kITSrefit);
0ddbf657 358 fCountITSrefit++;
e7d83d38 359 }
360 else {
0ddbf657 361 AliDebug(2,Form("Refit Failed for track %d |ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 362 }
363 }
32d38de2 364 //
0ddbf657 365 AliInfo(Form("%d ITSrefit in %d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
366 fCountITSrefit,fCountITSout,fCountITSin,fCountProlongationTrials,nTrESD));
367 //
32d38de2 368 return 0;
369}
370
371//_________________________________________________________________________
372Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
373{
374 // read from tree (if pointer provided) or directly from the ITS reco interface
375 //
376 return fReconstructor->LoadClusters(treeRP);
377}
378
379//_________________________________________________________________________
380void AliITSUTrackerGlo::UnloadClusters()
381{
382 //
383 // To be implemented
384 //
385
386 Info("UnloadClusters","To be implemented");
387}
388//_________________________________________________________________________
389AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
390{
391 //
392 // To be implemented
393 //
394
395 Info("GetCluster","To be implemented");
396 return 0x0;
397}
398
399//_________________________________________________________________________
400Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
401{
402 // do we need to match this track to ITS?
403 //
404 static double bz = GetBz();
405 if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
406 esdTr->IsOn(AliESDtrack::kTPCout) ||
407 esdTr->IsOn(AliESDtrack::kITSin) ||
408 esdTr->GetKinkIndex(0)>0) return kFALSE;
409 //
c61e50c3 410 if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
32d38de2 411 //
412 float dtz[2];
413 esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz);
414 // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
c61e50c3 415 if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation())
416 && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
417 Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
32d38de2 418 //
419 return kTRUE;
420}
421
422//_________________________________________________________________________
716ccba7 423void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
32d38de2 424{
425 // find prolongaion candidates finding for single seed
426 //
3e4e3c23 427 AliITSUSeed seedUC; // copy of the seed from the upper layer
428 AliITSUSeed seedT; // transient seed between the seedUC and new prolongation hypothesis
429 //
3e4e3c23 430 if (!InitHypothesis(esdTr,esdID)) return; // initialize prolongations hypotheses tree
32d38de2 431 //
173b3073 432 AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
f8832015 433 //
c61e50c3 434 TObjArray clArr; // container for transfer of clusters matching to seed
32d38de2 435 //
3e4e3c23 436 int nLrActive = fITS->GetNLayersActive();
437 for (int ila=nLrActive;ila--;) {
32d38de2 438 int ilaUp = ila+1; // prolong seeds from layer above
3e4e3c23 439 //
440 // for the outermost layer the seed is created from the ESD track
441 int nSeedsUp = (ilaUp==nLrActive) ? 1 : fCurrHyp->GetNSeeds(ilaUp);
38997c3c 442 int maxNBranches = fCurrTrackCond->GetMaxBranches(ila);
443 int maxNCandidates = fCurrTrackCond->GetMaxCandidates(ila);
3e4e3c23 444 //
32d38de2 445 for (int isd=0;isd<nSeedsUp;isd++) {
3e4e3c23 446 AliITSUSeed* seedU;
447 if (ilaUp==nLrActive) {
448 seedU = 0;
449 seedUC.InitFromESDTrack(esdTr);
450 }
451 else {
452 seedU = fCurrHyp->GetSeed(ilaUp,isd); // seed on prev.active layer to prolong
453 seedUC = *seedU; // its copy will be prolonged
b8b59e05 454 seedUC.SetParent(seedU);
3e4e3c23 455 }
943e1898 456 seedUC.ResetFMatrix(); // reset the matrix for propagation to next layer
32d38de2 457 // go till next active layer
0ddbf657 458 AliDebug(2,Form("working on Lr:%d Seed:%d of %d for esdID=%d (MClb:%d) | pT=%.3f",ila,isd,nSeedsUp,esdID,fCurrESDtrMClb,seedUC.Pt()));
f8832015 459 if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) {
32d38de2 460 //
0ddbf657 461 AliDebug(2,Form("Transport failed | esdID=%d (MClb:%d)",esdID,fCurrESDtrMClb));
32d38de2 462 // Check if the seed satisfies to track definition
b8b59e05 463 if (NeedToKill(&seedUC,kTransportFailed) && seedU) KillSeed(seedU,kTRUE);
32d38de2 464 continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
465 }
466 AliITSURecoLayer* lrA = fITS->GetLayerActive(ila);
f8832015 467 if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
b8b59e05 468 if (NeedToKill(&seedUC,kRWCheckFailed) && seedU) KillSeed(seedU,kTRUE);
32d38de2 469 continue;
470 }
173b3073 471 int nsens = lrA->FindSensors(&fTrImpData[kTrPhi0], hitSens); // find detectors which may be hit by the track
0ddbf657 472 AliDebug(2,Form("Will check %d sensors on lr:%d | esdID=%d (MClb:%d)",nsens,ila,esdID,fCurrESDtrMClb));
32d38de2 473 //
38997c3c 474 seedUC.SetLr(ila);
475 //
943e1898 476 double bz = GetBz();
32d38de2 477 for (int isn=nsens;isn--;) {
c61e50c3 478 AliITSURecoSens* sens = hitSens[isn];
38997c3c 479 int ncl = sens->GetNClusters();
480 if (!ncl) continue;
481 seedT = seedUC;
f8832015 482 //
943e1898 483 // We need to propagate the seed to sensor on lrA staying the frame of the sensor from prev.layer,
484 // since the transport matrix should be defined in this frame.
485 double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
486 if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) continue;
8b16dbae 487 if (!PropagateSeed(&seedT,xs,fCurrMass)) continue;
488 // if (!seedT.PropagateToX(xs,bz)) continue;
716ccba7 489 // if (!seedT.Rotate(sens->GetPhiTF())) continue;
490 if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
943e1898 491 //
f8832015 492 int clID0 = sens->GetFirstClusterId();
38997c3c 493 for (int icl=ncl;icl--;) {
f8832015 494 int res = CheckCluster(&seedT,ila,clID0+icl);
495 //
496 if (res==kStopSearchOnSensor) break; // stop looking on this sensor
497 if (res==kClusterNotMatching) continue; // cluster does not match
498 // cluster is matching and it was added to the hypotheses tree
c61e50c3 499 }
32d38de2 500 }
38997c3c 501 // cluster search is done. Do we need to have a version of this seed skipping current layer
502 if (!NeedToKill(&seedUC,kMissingCluster)) {
b2d3bc3a 503 AliITSUSeed* seedSkp = NewSeedFromPool(&seedUC);
716ccba7 504 double penalty = -AliITSUReconstructor::GetRecoParam()->GetMissPenalty(ila);
505 // to do: make penalty to account for probability to miss the cluster for good reason
506 seedSkp->SetChi2Cl(penalty);
507 AddProlongationHypothesis(seedSkp,ila);
508 }
38997c3c 509 // transfer the new branches of the seed to the hypothesis container
510 if (fNBranchesAdded) ValidateAllowedBranches(maxNBranches);
511 //
32d38de2 512 }
38997c3c 513 if (fNCandidatesAdded) ValidateAllowedCandidates(ila,maxNCandidates);
514 // ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
515 if (AliDebugLevelClass()>2) { //RS
15b02d69 516 printf(">>> All hypotheses on lr %d: \n",ila);
517 for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();}
518 }
519 //
520 /*
716ccba7 521 for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
716ccba7 522 if (ila!=0) continue;
523 double vecL[5] = {0};
524 double matL[15] = {0};
525 AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
526 while(sp->GetParent()) {
527 sp->Smooth(vecL,matL);
528 if (sp->GetLayerID()>=fITS->GetNLayersActive()-1) break;
529 sp = (AliITSUSeed*)sp->GetParent();
530 }
531 */
32d38de2 532 }
533 //
716ccba7 534 SaveCurrentTrackHypotheses();
535 //
32d38de2 536}
537
538//_________________________________________________________________________
3e4e3c23 539Bool_t AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esdID)
32d38de2 540{
c61e50c3 541 // init prolongaion candidates finding for single seed
716ccba7 542 fCurrHyp = GetTrackHyp(esdID);
543 if (fCurrHyp) return kTRUE;
544 //
0ddbf657 545 fCountProlongationTrials++;
546 //
32d38de2 547 fCurrMass = esdTr->GetMass();
548 if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
c61e50c3 549 //
3e4e3c23 550 fCurrHyp = new AliITSUTrackHyp(fITS->GetNLayersActive());
551 fCurrHyp->SetESDTrack(esdTr);
716ccba7 552 fCurrHyp->SetUniqueID(esdID);
c51c3124 553 fCurrHyp->SetMass(fCurrMass);
716ccba7 554 SetTrackHyp(fCurrHyp,esdID);
3e4e3c23 555 //
32d38de2 556 return kTRUE;
557 // TO DO
558}
559
560//_________________________________________________________________________
561Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
562{
0ddbf657 563 // transport seed from layerFrom to the entrance (along the direction of the track) of layerTo
32d38de2 564 //
8b16dbae 565 //
566 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
567 //
32d38de2 568 int dir = lTo > lFrom ? 1:-1;
569 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
570 Bool_t checkFirst = kTRUE;
571 while(lFrom!=lTo) {
32d38de2 572 if (lrFr) {
b2d3bc3a 573 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) {
0ddbf657 574 //printf("FailHere0 Dir=%d\n",dir);
575 //seed->Print("etp");
b2d3bc3a 576 return kFALSE; // go till the end of current layer
577 }
e2d2481c 578 checkFirst = kFALSE;
32d38de2 579 }
580 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
581 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
582 //
583 // go the entrance of the layer, assuming no materials in between
8b16dbae 584 double xToGo = lrTo->GetR(-dir);
0ddbf657 585 // double xts = xToGo;
b2d3bc3a 586 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
0ddbf657 587 //printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
588 //seed->Print("etp");
589 //
590
b2d3bc3a 591 return kFALSE;
592 }
593 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
0ddbf657 594 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
595 //seed->Print("etp");
b2d3bc3a 596 return kFALSE;
597 }
32d38de2 598 lrFr = lrTo;
599 }
600 return kTRUE;
601 //
602}
603
3e4e3c23 604//_________________________________________________________________________
605Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo)
606{
607 // transport track from layerFrom to the entrance of layerTo
608 //
8b16dbae 609 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
610 //
3e4e3c23 611 int dir = lTo > lFrom ? 1:-1;
612 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
613 Bool_t checkFirst = kTRUE;
614 while(lFrom!=lTo) {
3e4e3c23 615 if (lrFr) {
e2d2481c 616 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
617 checkFirst = kFALSE;
3e4e3c23 618 }
619 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
620 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
621 //
0ddbf657 622 // go the entrance of the layer, assuming no materials in between
623 double xToGo = lrTo->GetR(-dir);
624 // double xts = xToGo;
625 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
626 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
627 // seed->Print("etp");
628 return kFALSE;
629 }
630 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
631 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
632 //seed->Print("etp");
633 return kFALSE;
634 }
3e4e3c23 635 lrFr = lrTo;
636 }
637 return kTRUE;
638 //
639}
640
e2d2481c 641//_________________________________________________________________________
642Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
643{
644 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
645 // If check is requested, do this only provided the track has not exited the layer already
646 double xToGo = lr->GetR(dir);
647 if (check) { // do we need to track till the surface of the current layer ?
648 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
15b02d69 649 AliDebug(3,Form(" dir:%d Cur: %e Tgt: %e",dir,Sqrt(curR2),xToGo));
0ddbf657 650 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
651 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 652 }
653 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
654 // go via layer to its boundary, applying material correction.
655 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 656 //
e2d2481c 657 return kTRUE;
658 //
659}
660
661//_________________________________________________________________________
662Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
663{
664 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
665 // If check is requested, do this only provided the track has not exited the layer already
666 double xToGo = lr->GetR(dir);
667 if (check) { // do we need to track till the surface of the current layer ?
668 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 669 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
670 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 671 }
672 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
673 // go via layer to its boundary, applying material correction.
674 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 675 //
e2d2481c 676 return kTRUE;
677 //
678}
679
680
681//_________________________________________________________________________
682Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
683{
684 // go to the entrance of lr in direction dir, w/o applying material corrections.
685 // If check is requested, do this only provided the track did not reach the layer already
686 double xToGo = lr->GetR(-dir);
687 if (check) { // do we need to track till the surface of the current layer ?
688 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 689 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
690 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 691 }
692 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
693 // go via layer to its boundary, applying material correction.
694 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
695 return kTRUE;
696 //
697}
698
699//_________________________________________________________________________
700Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
701{
702 // go to the entrance of lr in direction dir, w/o applying material corrections.
703 // If check is requested, do this only provided the track did not reach the layer already
704 double xToGo = lr->GetR(-dir);
705 if (check) { // do we need to track till the surface of the current layer ?
706 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 707 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
708 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 709 }
710 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
711 // go via layer to its boundary, applying material correction.
712 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
713 return kTRUE;
714 //
715}
716
32d38de2 717//_________________________________________________________________________
c61e50c3 718Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
32d38de2 719{
720 // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
721 // as well as some aux info
722 double bz = GetBz();
723 AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
c61e50c3 724 seed->GetXYZ(&fTrImpData[kTrXIn]); // lab position at the entrance from above
943e1898 725 static AliExternalTrackParam sc; // seed copy for manipalitions
726 sc = *seed;
32d38de2 727 //
c61e50c3 728 fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
943e1898 729 if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
c61e50c3 730 double dr = lrA->GetDR(); // approximate X dist at the inner radius
943e1898 731 if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
32d38de2 732 // special case: track does not reach inner radius, might be tangential
943e1898 733 double r = sc.GetD(0,0,bz);
32d38de2 734 double x;
943e1898 735 if (!sc.GetXatLabR(r,x,bz,-1)) {
736 sc.Print();
737 AliFatal(Form("This should not happen: r=%f",r));
32d38de2 738 }
943e1898 739 dr = Abs(sc.GetX() - x);
740 if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) {
741 sc.Print();
742 AliFatal(Form("This should not happen: x=%f",x));
32d38de2 743 }
744 }
745 //
c61e50c3 746 fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
943e1898 747 double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
748 double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
c61e50c3 749 sgy = Sqrt(sgy)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY();
750 sgz = Sqrt(sgz)*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ();
751 fTrImpData[kTrPhi0] = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
173b3073 752 fTrImpData[kTrZ0] = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
c61e50c3 753 fTrImpData[kTrDPhi] = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]) + sgy/lrA->GetR();
173b3073 754 fTrImpData[kTrDZ] = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn]) + sgz;
32d38de2 755 //
756 return kTRUE;
757}
758
b8b59e05 759//________________________________________
760void AliITSUTrackerGlo::ResetSeedsPool()
761{
762 // mark all seeds in the pool as unused
0ddbf657 763 AliDebug(-1,Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
b8b59e05 764 fNFreeSeeds = 0;
765 fSeedsPool.Clear(); // seeds don't allocate memory
766}
767
768
769//________________________________________
770void AliITSUTrackerGlo::MarkSeedFree(AliITSUSeed *sd)
32d38de2 771{
b8b59e05 772 // account that this seed is "deleted"
773 int id = sd->GetPoolID();
774 if (id<0) {
775 AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
776 return;
777 }
778 // AliInfo(Form("%d %p",id, seed));
779 fSeedsPool.RemoveAt(id);
780 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
781 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
32d38de2 782}
f8832015 783
784//_________________________________________________________________________
785Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
786{
787 // Check if the cluster (in tracking frame!) is matching to track.
788 // The track must be already propagated to sensor tracking frame.
789 // Returns: kStopSearchOnSensor if the search on given sensor should be stopped,
790 // kClusterMatching if the cluster is matching
791 // kClusterMatching otherwise
792 //
f8832015 793 const double kTolerX = 5e-4;
794 AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
795 //
c8d1f258 796 Bool_t goodCl = kFALSE;
0ddbf657 797 int currLabel = Abs(fCurrESDtrMClb);
716ccba7 798 //
8b16dbae 799 if (cl->GetLabel(0)>=0) {
800 for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}
801 }
716ccba7 802 else goodCl = kTRUE;
c8d1f258 803 //
f8832015 804 if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
c8d1f258 805 if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
0ddbf657 806 if (goodCl && AliDebugLevelClass()>2) {
807 AliDebug(2,Form("Lost good cl: Failed propagation. |ESDtrack#%d (MClb:%d)",fCurrESDtrack->GetID(),fCurrESDtrMClb));
15b02d69 808 cl->Print();
809 }
c8d1f258 810 return kStopSearchOnSensor; // propagation failed, seedT is intact
811 }
f8832015 812 }
3dd9c283 813 double dy = cl->GetY()-track->GetY();
c8d1f258 814 double dz = cl->GetZ()-track->GetZ();
0091e9f0 815 //
76390254 816#ifdef _FILL_CONTROL_HISTOS_
817 int hcOffs = fTrackPhase*kHistosPhase + lr;
818 double htrPt=-1;
819 if (goodCl && cHistoArr && track->GetChi2Penalty()<1e-5 && !track->ContainsFake()/* && track->Charge()>0*/) {
820 htrPt = track->Pt();
821 ((TH2*)cHistoArr->At(kHResY+hcOffs))->Fill(htrPt,dy);
822 ((TH2*)cHistoArr->At(kHResZ+hcOffs))->Fill(htrPt,dz);
823 double errY = track->GetSigmaY2();
824 double errZ = track->GetSigmaZ2();
825 if (errY>0) ((TH2*)cHistoArr->At(kHResYP+hcOffs))->Fill(htrPt,dy/Sqrt(errY));
826 if (errZ>0) ((TH2*)cHistoArr->At(kHResZP+hcOffs))->Fill(htrPt,dz/Sqrt(errZ));
827 }
828#endif
829 //
f8832015 830 double dy2 = dy*dy;
831 double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
832 AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadY(); // RS TOOPTIMIZE
833 if (dy2>tol2) { // the clusters are sorted in Z(col) then in Y(row).
15b02d69 834 if (goodCl && AliDebugLevelClass()>2) {
0ddbf657 835 AliDebug(2,Form("Loose good cl: dy2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",dy2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
15b02d69 836 cl->Print();
837 }
f8832015 838 if (dy>0) return kStopSearchOnSensor; // No chance that other cluster of this sensor will match (all Y's will be even larger)
839 else return kClusterNotMatching; // Other clusters may match
840 }
c8d1f258 841 double dz2 = dz*dz;
f8832015 842 tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
843 AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ()*AliITSUReconstructor::GetRecoParam()->GetNSigmaRoadZ(); // RS TOOPTIMIZE
c8d1f258 844 if (dz2>tol2) {
15b02d69 845 if (goodCl && AliDebugLevelClass()>2) {
0ddbf657 846 AliDebug(2,Form("Lost good cl: dz2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",dz2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
15b02d69 847 cl->Print();
848 }
c8d1f258 849 return kClusterNotMatching; // Other clusters may match
850 }
f8832015 851 //
852 // check chi2
853 Double_t p[2]={cl->GetY(), cl->GetZ()};
854 Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
855 double chi2 = track->GetPredictedChi2(p,cov);
76390254 856 //
857#ifdef _FILL_CONTROL_HISTOS_
858 if (htrPt>0) {
859 ((TH2*)cHistoArr->At(kHChi2Cl+hcOffs))->Fill(htrPt,chi2);
860 }
861#endif
862 //
c8d1f258 863 if (chi2>AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr)) {
15b02d69 864 if (goodCl && AliDebugLevelClass()>2) {
0ddbf657 865 AliDebug(2,Form("Lost good cl: Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
866 chi2,AliITSUReconstructor::GetRecoParam()->GetMaxTr2ClChi2(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb));
3dd9c283 867 track->Print("etp");
868 cl->Print("");
869 }
c8d1f258 870 return kClusterNotMatching;
871 }
f8832015 872 //
873 track = NewSeedFromPool(track); // input track will be reused, use its clone for updates
943e1898 874 if (!track->Update()) {
38997c3c 875 if (goodCl && AliDebugLevelClass()>2) {
876 AliDebug(2,"Lost good cluster: failed to update");
877 track->Print("etp");
878 cl->Print("");
879 printf("Loose good cl: Failed update |");
880 }
b8b59e05 881 MarkSeedFree(track);
c8d1f258 882 return kClusterNotMatching;
883 }
f8832015 884 track->SetChi2Cl(chi2);
885 track->SetLrClusterID(lr,clID);
886 cl->IncreaseClusterUsage();
887 //
716ccba7 888 track->SetFake(!goodCl);
889 //
0ddbf657 890 AliDebug(2,Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)| ",
716ccba7 891 goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
c8d1f258 892 //
38997c3c 893 AddSeedBranch(track);
f8832015 894 //
895 return kClusterMatching;
896}
c8d1f258 897
898//_________________________________________________________________________
899Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
900{
901 // check if the seed should not be discarded
902 const UShort_t kMask = 0xffff;
903 if (flag==kMissingCluster) {
904 int lastChecked = seed->GetLayerID();
905 UShort_t patt = seed->GetHitsPattern();
906 if (lastChecked) patt |= ~(kMask<<lastChecked); // not all layers were checked, complete unchecked once by potential hits
42c3d4bd 907 Bool_t seedOK = fCurrTrackCond->CheckPattern(patt);
c8d1f258 908 return !seedOK;
909 }
910 return kTRUE;
911}
44785f3e 912
913//______________________________________________________________________________
914Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
915{
916 // propagate seed to given x applying material correction if requested
917 const Double_t kEpsilon = 1e-5;
918 Double_t xpos = seed->GetX();
919 Int_t dir = (xpos<xToGo) ? 1:-1;
920 Double_t xyz0[3],xyz1[3],param[7];
921 //
8b16dbae 922 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
f497470c 923 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
44785f3e 924 while ( (xToGo-xpos)*dir > kEpsilon){
925 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
926 Double_t x = xpos+step;
927 Double_t bz=GetBz(); // getting the local Bz
0ddbf657 928 if (!seed->PropagateToX(x,bz)) return kFALSE;
8b16dbae 929 double ds = 0;
930 if (matCorr || updTime) {
44785f3e 931 xyz0[0]=xyz1[0]; // global pos at the beginning of step
932 xyz0[1]=xyz1[1];
933 xyz0[2]=xyz1[2];
934 seed->GetXYZ(xyz1); // // global pos at the end of step
8b16dbae 935 if (matCorr) {
936 MeanMaterialBudget(xyz0,xyz1,param);
937 Double_t xrho=param[0]*param[4], xx0=param[1];
938 if (dir>0) xrho = -xrho; // outward should be negative
b2d3bc3a 939 if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) {AliInfo("Fail2"); seed->Print("etp"); return kFALSE;}
8b16dbae 940 ds = param[4];
941 }
942 else { // matCorr is not requested but time integral is
943 double d0 = xyz1[0]-xyz0[0];
944 double d1 = xyz1[1]-xyz0[1];
945 double d2 = xyz1[2]-xyz0[2];
946 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
947 }
44785f3e 948 }
8b16dbae 949 if (updTime) seed->AddTimeStep(ds);
44785f3e 950 xpos = seed->GetX();
951 }
952 return kTRUE;
953}
716ccba7 954
3e4e3c23 955//______________________________________________________________________________
956Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
957{
958 // propagate seed to given x applying material correction if requested
959 const Double_t kEpsilon = 1e-5;
960 Double_t xpos = seed->GetX();
961 Int_t dir = (xpos<xToGo) ? 1:-1;
962 Double_t xyz0[3],xyz1[3],param[7];
963 //
964 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
965 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
966 while ( (xToGo-xpos)*dir > kEpsilon){
967 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
968 Double_t x = xpos+step;
969 Double_t bz=GetBz(); // getting the local Bz
970 if (!seed->PropagateTo(x,bz)) return kFALSE;
971 double ds = 0;
972 if (matCorr || updTime) {
973 xyz0[0]=xyz1[0]; // global pos at the beginning of step
974 xyz0[1]=xyz1[1];
975 xyz0[2]=xyz1[2];
976 seed->GetXYZ(xyz1); // // global pos at the end of step
977 //
978 if (matCorr) {
979 MeanMaterialBudget(xyz0,xyz1,param);
980 Double_t xrho=param[0]*param[4], xx0=param[1];
981 if (dir>0) xrho = -xrho; // outward should be negative
982 if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
983 ds = param[4];
984 }
985 else { // matCorr is not requested but time integral is
986 double d0 = xyz1[0]-xyz0[0];
987 double d1 = xyz1[1]-xyz0[1];
988 double d2 = xyz1[2]-xyz0[2];
989 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
990 }
991 }
992 if (updTime) seed->AddTimeStep(ds);
993 //
994 xpos = seed->GetX();
995 }
996 return kTRUE;
997}
998
716ccba7 999//______________________________________________________________________________
1000void AliITSUTrackerGlo::SaveCurrentTrackHypotheses()
1001{
1002 // RS: shall we clean up killed seeds?
38997c3c 1003 ((TObjArray*)fCurrHyp->GetLayerSeeds(0))->Sort();
716ccba7 1004 fCurrHyp = 0;
1005 // TODO
1006
1007}
3e4e3c23 1008
1009//______________________________________________________________________________
1010void AliITSUTrackerGlo::FinalizeHypotheses()
1011{
1012 // select winner for each hypothesis, remove cl. sharing conflicts
15b02d69 1013 AliDebug(2,"TODO");
3e4e3c23 1014 //
1015 int nh = fHypStore.GetEntriesFast();
1016 for (int ih=0;ih<nh;ih++) {
1017 AliITSUTrackHyp* hyp = (AliITSUTrackHyp*) fHypStore.UncheckedAt(ih);
69e0f089 1018 if (!hyp || !hyp->DefineWinner()) continue; // TODO
0ddbf657 1019 if (hyp->GetESDTrack()->IsOn(AliESDtrack::kITSin)) continue;
1020 fCountITSin++;
68a0f687 1021 CookMCLabel(hyp);
1022 UpdateESDTrack(hyp,AliESDtrack::kITSin);
3e4e3c23 1023 }
1024
1025}
c51c3124 1026
1027//______________________________________________________________________________
68a0f687 1028void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp,Int_t flag)
c51c3124 1029{
1030 // update ESD track with current best hypothesis
1031 AliESDtrack* esdTr = hyp->GetESDTrack();
1032 if (!esdTr) return;
1033 AliITSUSeed* win = hyp->GetWinner();
1034 if (!win) return;
68a0f687 1035 switch (flag) {
1036 case AliESDtrack::kITSin:
1037 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1038 // TODO: set cluster info
1039 break;
1040 //
1041 case AliESDtrack::kITSout:
1042 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1043 // TODO: avoid setting friend
1044 break;
1045 //
1046 case AliESDtrack::kITSrefit:
1047 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1048 // TODO: avoid setting cluster info
1049 break;
1050 default:
1051 AliFatal(Form("Unknown flag %d",flag));
1052 }
c51c3124 1053 //
1054 // transfer module indices
1055 // TODO
1056}
1057
1058//______________________________________________________________________________
68a0f687 1059Bool_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest, Bool_t stopAtLastCl)
c51c3124 1060{
1061 // refit track till radius rDest
8b16dbae 1062 AliITSUTrackHyp tmpTr;
1063 //
c51c3124 1064 double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
1065 int dir,lrStart,lrStop;
1066 //
8b16dbae 1067 dir = rCurr<rDest ? 1 : -1;
1068 lrStart = fITS->FindFirstLayerID(rCurr,dir);
1069 lrStop = fITS->FindLastLayerID(rDest,dir);
15b02d69 1070 if (AliDebugLevelClass()>2) {
1071 printf("Refit %d: Lr:%d (%f) -> Lr:%d (%f)\n",dir,lrStart,rCurr, lrStop,rDest);
1072 printf("Before refit: "); trc->AliExternalTrackParam::Print();
1073 }
8b16dbae 1074 if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
1075 //
68a0f687 1076 int ncl = trc->FetchClusterInfo(fClInfo);
8b16dbae 1077 fCurrMass = trc->GetMass();
1078 tmpTr.AliKalmanTrack::operator=(*trc);
1079 tmpTr.SetChi2(0);
68a0f687 1080 int iclLr[2],nclLr,clCount=0;
c51c3124 1081 //
1082 for (int ilr=lrStart;ilr!=lrStop;ilr+=dir) {
1083 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
1084 if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
8b16dbae 1085 int ilrA2,ilrA = lr->GetActiveID();
1086 // passive layer or active w/o hits will be traversed on the way to next cluster
1087 if (!lr->IsActive() || fClInfo[ilrA2=(ilrA<<1)]<0) continue;
1088 //
e2d2481c 1089 if (ilr!=lrStart && !TransportToLayer(&tmpTr,lrStart,ilr)) {
0ddbf657 1090 AliDebug(2,Form("Failed to transport %d -> %d | ESDtrack#%d (MClb:%d)\n",lrStart,ilr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1091 //tmpTr.AliExternalTrackParam::Print();
1092 //trc->GetWinner()->Print("etp");
e2d2481c 1093 return kFALSE; // go to the entrance to the layer
1094 }
8b16dbae 1095 lrStart = ilr;
1096 //
1097 // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
1098 nclLr=0;
1099 if (dir>0) { // clusters are stored in increasing radius order
1100 iclLr[nclLr++]=fClInfo[ilrA2++];
1101 if (fClInfo[ilrA2]>=0) iclLr[nclLr++]=fClInfo[ilrA2];
1102 }
1103 else {
1104 if ( fClInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=fClInfo[ilrA2+1];
1105 iclLr[nclLr++]=fClInfo[ilrA2];
1106 }
1107 //
1108 for (int icl=0;icl<nclLr;icl++) {
1109 AliITSUClusterPix* clus = (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
1110 AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
e7d83d38 1111 if (!tmpTr.Rotate(sens->GetPhiTF())) {
0ddbf657 1112 AliDebug(2,Form("Failed on rotate to %f | ESDtrack#%d (MClb:%d)",sens->GetPhiTF(),fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 1113 return kFALSE;
1114 }
1115 //printf("Refit cl:%d on lr %d Need to go %.4f -> %.4f\n",icl,ilrA,tmpTr.GetX(),sens->GetXTF()+clus->GetX());
1116 if (!PropagateSeed(&tmpTr,sens->GetXTF()+clus->GetX(),fCurrMass)) {
0ddbf657 1117 AliDebug(2,Form("Failed on propagate to %f (dir=%d) | ESDtrack#%d (MClb:%d)",sens->GetXTF()+clus->GetX(),dir,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1118 //tmpTr.AliExternalTrackParam::Print("");
1119 //trc->GetWinner()->Print("etp");
e7d83d38 1120 return kFALSE;
1121 }
76390254 1122 //
1123#ifdef _FILL_CONTROL_HISTOS_
1124 int hcOffs = fTrackPhase*kHistosPhase + ilrA;
1125 double htrPt=-1;
1126 if (cHistoArr && trc->GetLabel()>=0/* && trc->Charge()>0*/) {
1127 htrPt = tmpTr.Pt();
1128 double dy = clus->GetY()-tmpTr.GetY();
1129 double dz = clus->GetZ()-tmpTr.GetZ();
1130 ((TH2*)cHistoArr->At(kHResY+hcOffs))->Fill(htrPt,dy);
1131 ((TH2*)cHistoArr->At(kHResZ+hcOffs))->Fill(htrPt,dz);
1132 double errY = tmpTr.GetSigmaY2();
1133 double errZ = tmpTr.GetSigmaZ2();
1134 if (errY>0) ((TH2*)cHistoArr->At(kHResYP+hcOffs))->Fill(htrPt,dy/Sqrt(errY));
1135 if (errZ>0) ((TH2*)cHistoArr->At(kHResZP+hcOffs))->Fill(htrPt,dz/Sqrt(errZ));
1136 }
1137#endif
1138 //
1139 double chi2;
1140 if ( (chi2=tmpTr.Update(clus))<0 ) {
0ddbf657 1141 AliDebug(2,Form("Failed on Update | ESDtrack#%d (MClb:%d)",fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 1142 return kFALSE;
1143 }
76390254 1144#ifdef _FILL_CONTROL_HISTOS_
1145 if (htrPt>0) {
1146 ((TH2*)cHistoArr->At(kHChi2Cl+hcOffs))->Fill(htrPt,chi2);
1147 }
1148#endif
e7d83d38 1149 //printf("AfterRefit: "); tmpTr.AliExternalTrackParam::Print();
68a0f687 1150 if (stopAtLastCl && ++clCount==ncl) return kTRUE; // it was requested to not propagate after last update
8b16dbae 1151 }
c51c3124 1152 //
1153 }
8b16dbae 1154 // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
1155 // Still, try to go as close as possible to rDest.
1156 //
1157 if (lrStart!=lrStop) {
e7d83d38 1158 //printf("Going to last layer %d -> %d\n",lrStart,lrStop);
1159 if (!TransportToLayer(&tmpTr,lrStart,lrStop)) {
0ddbf657 1160 AliDebug(1,Form("Failed on TransportToLayer %d->%d | ESDtrack#%d (MClb:%d)",lrStart,lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1161 //tmpTr.AliExternalTrackParam::Print();
1162 //trc->GetWinner()->Print("etp");
e7d83d38 1163 return kTRUE;
1164 }
1165 if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) {
0ddbf657 1166 AliDebug(3,Form("Failed on GoToExitFromLayer %d | ESDtrack#%d (MClb:%d)",lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 1167 return kTRUE; // go till the exit from layer
1168 }
e2d2481c 1169 //
e7d83d38 1170 //printf("On exit from last layer\n");
15b02d69 1171 //tmpTr.AliExternalTrackParam::Print();
e2d2481c 1172 // go to the destination radius
1173 if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),dir)) return kTRUE;
1174 if (!PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) return kTRUE;
8b16dbae 1175 }
1176 trc->AliKalmanTrack::operator=(tmpTr);
15b02d69 1177 if (AliDebugLevelClass()>2) {
1178 printf("After refit (now at lr %d): ",lrStart); trc->AliExternalTrackParam::Print();
1179 }
8b16dbae 1180 return kTRUE;
c51c3124 1181}
dd2117a2 1182
1183//__________________________________________________________________
1184void AliITSUTrackerGlo::CookMCLabel(AliITSUTrackHyp* hyp)
1185{
1186 // build MC label
1187 //
1188 const int kMaxLbPerCl = 3;
1189 int lbID[kMaxLayers*6],lbStat[kMaxLayers*6];
1190 Int_t lr,nLab=0,nCl=0;
1191 AliITSUSeed *seed = hyp->GetWinner();
1192 while(seed) {
1193 int clID = seed->GetLrCluster(lr);
1194 if (clID>=0) {
1195 AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1196 nCl++;
1197 for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
1198 int trLb = cl->GetLabel(imc);
704d9398 1199 if (trLb<0) break;
dd2117a2 1200 // search this mc track in already accounted ones
1201 int iLab;
1202 for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
1203 if (iLab<nLab) lbStat[iLab]++;
1204 else {
1205 lbID[nLab] = trLb;
1206 lbStat[nLab++] = 1;
1207 }
1208 } // loop over given cluster's labels
1209 }
1210 seed = (AliITSUSeed*)seed->GetParent();
1211 } // loop over clusters
1212 //
41343470 1213 if (nCl && nLab) {
dd2117a2 1214 int maxLab=0,nTPCok=0;
1215 AliESDtrack* esdTr = hyp->GetESDTrack();
1216 int tpcLab = esdTr ? Abs(esdTr->GetTPCLabel()) : -kDummyLabel;
1217 for (int ilb=nLab;ilb--;) {
1218 int st = lbStat[ilb];
1219 if (lbStat[maxLab]<st) maxLab=ilb;
1220 if (tpcLab==lbID[ilb]) nTPCok += st;
1221 }
1222 hyp->SetFakeRatio(1.-float(nTPCok)/nCl);
1223 hyp->SetLabel( nTPCok==nCl ? tpcLab : -tpcLab);
1224 hyp->SetITSLabel( lbStat[maxLab]==nCl ? lbStat[maxLab] : -lbStat[maxLab]); // winner label
1225 return;
1226 }
1227 //
1228 hyp->SetFakeRatio(-1.);
1229 hyp->SetLabel( kDummyLabel );
1230 hyp->SetITSLabel( kDummyLabel );
1231}
38997c3c 1232
1233//__________________________________________________________________
1234Bool_t AliITSUTrackerGlo::AddSeedBranch(AliITSUSeed* seed)
1235{
1236 // add new prolongation branch for the seed to temporary array. The seeds are sorted according to Compare f-n
1237 if (fNCandidatesAdded+fNBranchesAdded>=AliITSUTrackCond::kMaxCandidates ) {
1238 AliError(Form("Number of candidates for current seed reached limit of AliITSUTrackCond::kMaxCandidates=%d, increase it",AliITSUTrackCond::kMaxCandidates));
1239 return kFALSE;
1240 }
1241 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded]; // note: fNCandidatesAdded is incremented after adding all branches of current seed
1242 int slot=fNBranchesAdded++;
1243 for (int slotF=slot;slotF--;) { // slotF is always slot-1
1244 AliITSUSeed* si = branches[slotF];
1245 if (si->Compare(seed)<0) break; // found the last seed with better quality
1246 // otherwise, shift the worse seed to the next slot
1247 branches[slot] = si;
1248 slot = slotF; // slot should be slotF+1
1249 }
1250 // if needed, move worse seeds
1251 branches[slot] = seed;
1252 return kTRUE;
1253 //
1254}
1255
1256//__________________________________________________________________
1257void AliITSUTrackerGlo::ValidateAllowedBranches(Int_t acceptMax)
1258{
1259 // keep allowed number of branches for current seed and disable extras
1260 int nb = Min(fNBranchesAdded,acceptMax);
1261 // if (nb<fNBranchesAdded) printf("ValidateAllowedBranches: %d of %d (%d) on lr %d\n",nb,fNBranchesAdded,fNCandidatesAdded,ilr);
1262 // disable unused branches
1263 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded];
1264 for (int ib=nb;ib<fNBranchesAdded;ib++) {
1265 /*
1266 if (!branches[ib]->ContainsFake()) {
1267 printf("Suppress good branch as %d of %d |",ib,fNBranchesAdded); branches[ib]->Print();
1268 printf("Survivors : \n");
1269 for (int j=0;j<nb;j++) branches[j]->Print();
1270 }
1271 */
1272 MarkSeedFree(branches[ib]);
1273 }
1274 fNCandidatesAdded += nb; // update total candidates counter
1275 fNBranchesAdded = 0; // reset branches counter
1276 //
1277}
1278
1279//__________________________________________________________________
1280void AliITSUTrackerGlo::ValidateAllowedCandidates(Int_t ilr, Int_t acceptMax)
1281{
1282 // transfer allowed number of branches to hypothesis container
1283 //
1284 // sort candidates in increasing order of chi2
1285 Int_t index[AliITSUTrackCond::kMaxCandidates];
1286 Float_t chi2[AliITSUTrackCond::kMaxCandidates];
1287 for (int i=fNCandidatesAdded;i--;) chi2[i] = fLayerCandidates[i]->GetChi2GloNrm();
1288 Sort(fNCandidatesAdded,chi2,index,kFALSE);
1289 //
1290 int nb = Min(fNCandidatesAdded,acceptMax);
1291 // if (nb<fNCandidatesAdded) printf("ValidateAllowedCandidates: %d of %d on lr %d\n",nb,fNCandidatesAdded,ilr);
1292 //
1293 for (int ib=0;ib<nb;ib++) AddProlongationHypothesis(fLayerCandidates[index[ib]],ilr);
1294 // disable unused candidates
1295 for (int ib=nb;ib<fNCandidatesAdded;ib++) {
1296 /*
1297 if (!fLayerCandidates[index[ib]]->ContainsFake()) {
1298 printf("Suppress good candidate as %d of %d |",index[ib],fNCandidatesAdded); fLayerCandidates[index[ib]]->Print();
1299 }
1300 */
1301 MarkSeedFree(fLayerCandidates[index[ib]]);
1302 }
1303 fNCandidatesAdded = 0; // reset candidates counter
1304 //
1305}
1306