]> 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 cluster sharing histogramming
[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"
cb50e082 40#include "AliITSUGeomTGeo.h"
bdb92766 41#include "AliCodeTimer.h"
32d38de2 42using namespace AliITSUAux;
43using namespace TMath;
44
c8d1f258 45//----------------- tmp stuff -----------------
c8d1f258 46
32d38de2 47ClassImp(AliITSUTrackerGlo)
e2d2481c 48
49const Double_t AliITSUTrackerGlo::fgkToler = 1e-6;// tolerance for layer on-surface check
50
51
32d38de2 52//_________________________________________________________________________
53AliITSUTrackerGlo::AliITSUTrackerGlo(AliITSUReconstructor* rec)
54: fReconstructor(rec)
55 ,fITS(0)
0091e9f0 56 ,fCurrESDtrack(0)
0ddbf657 57 ,fCurrESDtrMClb(kDummyLabel)
32d38de2 58 ,fCurrMass(kPionMass)
0ddbf657 59 ,fCountProlongationTrials(0)
60 ,fCountITSin(0)
61 ,fCountITSout(0)
62 ,fCountITSrefit(0)
716ccba7 63 ,fHypStore(100)
6cd80116 64 ,fLayerMaxCandidates(1000)
65 ,fLayerCandidates(0)
38997c3c 66 ,fNBranchesAdded(0)
67 ,fNCandidatesAdded(0)
716ccba7 68 ,fCurrHyp(0)
6cd80116 69 ,fWorkHyp(0)
c61e50c3 70 ,fSeedsPool("AliITSUSeed",0)
b8b59e05 71 ,fFreeSeedsID(0)
72 ,fNFreeSeeds(0)
73 ,fLastSeedID(0)
6cd80116 74 ,fNLrActive(0)
42c3d4bd 75 ,fDefTrackConds(0)
76 ,fCurrTrackCond(0)
c03e4f8a 77 ,fCurrActLrID(-1)
78 ,fCurrLayer(0)
c51c3124 79 ,fTrackPhase(-1)
8b16dbae 80 ,fClInfo(0)
53870004 81#ifdef _FILL_CONTROL_HISTOS_
82 ,fCHistoArr(0)
83#endif
32d38de2 84{
85 // Default constructor
86 if (rec) Init(rec);
87}
88
89//_________________________________________________________________________
90AliITSUTrackerGlo::~AliITSUTrackerGlo()
91{
92 // Default destructor
93 //
8b16dbae 94 delete[] fClInfo;
32d38de2 95 //
76390254 96#ifdef _FILL_CONTROL_HISTOS_
53870004 97 if (fCHistoArr) {
76390254 98 TFile* ctrOut = TFile::Open("itsuControlHistos.root","recreate");
99 ctrOut->cd();
100 AliInfo("Storing control histos");
53870004 101 fCHistoArr->Print();
102 // ctrOut->WriteObject(fCHistoArr,"controlH","kSingleKey");
103 fCHistoArr->Write();
76390254 104 ctrOut->Close();
105 delete ctrOut;
53870004 106 fCHistoArr = 0;
76390254 107 }
108#endif
109 //
32d38de2 110}
111
112//_________________________________________________________________________
113void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
114{
115 // init with external reconstructor
c61e50c3 116 //
3d4dc3e2 117 fITS = rec->GetITSInterface();
6cd80116 118 fNLrActive = fITS->GetNLayersActive();
119 fClInfo = new Int_t[fNLrActive*2];
120 fWorkHyp = new AliITSUTrackHyp(fNLrActive);
c61e50c3 121 //
6cd80116 122 if (fLayerMaxCandidates<1) fLayerMaxCandidates = 1000;
123 fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
c61e50c3 124 fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
b8b59e05 125 fFreeSeedsID.Set(1000);
c61e50c3 126 //
32d38de2 127}
128
42c3d4bd 129//_________________________________________________________________________
130void AliITSUTrackerGlo::CreateDefaultTrackCond()
131{
132 // creates default tracking conditions to be used when recoparam does not provide them
6cd80116 133
42c3d4bd 134 AliITSUTrackCond* cond = new AliITSUTrackCond();
135 //
6cd80116 136 cond->SetNLayers(fNLrActive);
137 cond->AddNewCondition(fNLrActive);
42c3d4bd 138 cond->AddGroupPattern( 0xffff ); // require all layers hit
ee54014a 139 cond->Init();
42c3d4bd 140 //
141 fDefTrackConds.AddLast(cond);
142 //
143 AliInfo("Created conditions: ");
144 for (int i=0;i<fDefTrackConds.GetEntriesFast();i++) fDefTrackConds[i]->Print();
145 //
146}
147
148
32d38de2 149//_________________________________________________________________________
150Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
151{
152 //
bdb92766 153 AliCodeTimerAuto("",0);
76390254 154 SetTrackingPhase(kClus2Tracks);
155 //
156#ifdef _FILL_CONTROL_HISTOS_
53870004 157 if (!fCHistoArr) BookControlHistos();
76390254 158#endif
c03e4f8a 159 static int evID = 0;
bdb92766 160 static TArrayI esdTrInd;
161 static TArrayF esdTrPt;
162 //
42c3d4bd 163 TObjArray *trackConds = 0;
32d38de2 164 //
0ddbf657 165 fCountProlongationTrials = 0;
166 fCountITSin = 0;
167 fCountITSout = 0;
168 fCountITSrefit = 0;
169 //
b8b59e05 170 ResetSeedsPool();
716ccba7 171 int nTrESD = esdEv->GetNumberOfTracks();
c51c3124 172 AliInfo(Form("Will try to find prolongations for %d tracks",nTrESD));
42c3d4bd 173 int nTrackCond = AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();
174 if (nTrackCond<1) {
175 if (!fDefTrackConds.GetEntriesFast()) {
176 AliInfo("No tracking conditions found in recoparams, creating default one requesting all layers hit");
177 CreateDefaultTrackCond();
178 }
179 trackConds = &fDefTrackConds;
180 nTrackCond = trackConds->GetEntriesFast();
181 }
182 else trackConds = AliITSUReconstructor::GetRecoParam()->GetTrackingConditions();
183 //
15b02d69 184 static Bool_t first = kTRUE;
185 if (first) {
186 AliITSUReconstructor::GetRecoParam()->Print();
187 first = kFALSE;
188 }
189 fHypStore.Delete();
190 if (fHypStore.GetSize()<nTrESD) fHypStore.Expand(nTrESD+100);
191 //
192 fITS->ProcessClusters();
193 //
bdb92766 194#ifdef _FILL_CONTROL_HISTOS_
53870004 195 FlagSplitClusters(); // tmp RS
bdb92766 196#endif
197 //
198 // the tracks will be reconstructed in decreasing pt order, sort them
199 if (esdTrInd.GetSize()<nTrESD) {
200 esdTrInd.Set(nTrESD+200);
201 esdTrPt.Set(nTrESD+200);
202 }
203 int* esdTrackIndex = esdTrInd.GetArray();
204 float* trPt = esdTrPt.GetArray();
205 for (int itr=nTrESD;itr--;) trPt[itr] = esdEv->GetTrack(itr)->Pt();
206 Sort(nTrESD,trPt,esdTrackIndex,kTRUE);
53870004 207 //
42c3d4bd 208 for (int icnd=0;icnd<nTrackCond;icnd++) {
209 fCurrTrackCond = (AliITSUTrackCond*)trackConds->UncheckedAt(icnd);
ee54014a 210 if (!fCurrTrackCond->IsInitDone()) fCurrTrackCond->Init();
42c3d4bd 211 // select ESD tracks to propagate
212 for (int itr=0;itr<nTrESD;itr++) {
bdb92766 213 int trID = esdTrackIndex[itr];
214 fCurrESDtrack = esdEv->GetTrack(trID);
0ddbf657 215 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
216 //
217 if (!NeedToProlong(fCurrESDtrack)) continue; // are we interested in this track?
c03e4f8a 218 /*
219 // if speceific tracks should be checked, create a global array int wtc[] = {evITS*10000+trID, ...};
220 Bool_t dbg = kFALSE;
221 int nwtc = sizeof(wtc)/sizeof(int);
222
223 for (int k=0;k<nwtc;k++) {
bdb92766 224 if (wtc[k]==evID*10000+trID) {
c03e4f8a 225 dbg = kTRUE;
226 printf("\n\n\n\n\n\n\n");
227 printf("At esdTr: %d MC: %d\n",wtc[k],fCurrESDtrMClb);
228 printf("\n\n\n\n\n\n\n");
229 break;
230 }
231 }
232 AliLog::SetClassDebugLevel("AliITSUTrackerGlo",dbg ? 10:0);
233 */
234
bdb92766 235 AliDebug(1,Form("Processing track %d(esd%d) | M=%.3f Pt=%.3f | MCLabel: %d",trID,itr,fCurrESDtrack->GetMass(kTRUE),fCurrESDtrack->Pt(),fCurrESDtrMClb));//RS
236 FindTrack(fCurrESDtrack, trID);
42c3d4bd 237 }
238 //
38997c3c 239 if (AliDebugLevelClass()>+2) {
15b02d69 240 AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
241 fHypStore.Print();
242 }
42c3d4bd 243 FinalizeHypotheses();
c03e4f8a 244 //
245 //AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0); // in case wtc array was defined, uncomment this
42c3d4bd 246 }
716ccba7 247 //
0ddbf657 248 AliInfo(Form("%d ITSin for %d tried TPC seeds out of %d ESD tracks\n",fCountITSin,fCountProlongationTrials,nTrESD));
c03e4f8a 249 evID++;
32d38de2 250 return 0;
251}
252
253//_________________________________________________________________________
c51c3124 254Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent *esdEv)
32d38de2 255{
256 //
c51c3124 257 // Do outward fits in ITS
bdb92766 258 AliCodeTimerAuto("",0);
c51c3124 259 //
e7d83d38 260 SetTrackingPhase(kPropBack);
c51c3124 261 int nTrESD = esdEv->GetNumberOfTracks();
15b02d69 262 AliDebug(1,Form("Will propagate back %d tracks",nTrESD));
c51c3124 263 //
264 double bz0 = GetBz();
265 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
266 AliITSUTrackHyp dummyTr,*currTr=0;
267 const double kWatchStep=10.; // for larger steps watch arc vs segment difference
268 Double_t times[AliPID::kSPECIES];
269 //
270 //
271 for (int itr=0;itr<nTrESD;itr++) {
0ddbf657 272 fCurrESDtrack = esdEv->GetTrack(itr);
273 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
c51c3124 274 // Start time integral and add distance from current position to vertex
0ddbf657 275 if (fCurrESDtrack->IsOn(AliESDtrack::kITSout)) continue; // already done
c51c3124 276 //
0ddbf657 277 fCurrESDtrack->GetXYZ(xyzTrk);
c51c3124 278 Double_t dst = 0.; // set initial track lenght, tof
279 {
280 double dxs = xyzTrk[0] - xyzVtx[0];
281 double dys = xyzTrk[1] - xyzVtx[1];
282 double dzs = xyzTrk[2] - xyzVtx[2];
283 // RS: for large segment steps d use approximation of cicrular arc s by
284 // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
285 // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
286 // Hence s^2/d^2 = (1+1/6 p^2)^2
287 dst = dxs*dxs + dys*dys;
288 if (dst > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
0ddbf657 289 double crv = Abs(fCurrESDtrack->GetC(bz0));
c51c3124 290 double fcarc = 1.+crv*crv*dst/6.;
291 dst *= fcarc*fcarc;
292 }
293 dst += dzs*dzs;
294 dst = Sqrt(dst);
295 }
296 //
0ddbf657 297 fCurrESDtrack->SetStatus(AliESDtrack::kTIME);
c51c3124 298 //
0ddbf657 299 if (!fCurrESDtrack->IsOn(AliESDtrack::kITSin)) { // case of tracks w/o ITS prolongation: just set the time integration
300 dummyTr.AliExternalTrackParam::operator=(*fCurrESDtrack);
c51c3124 301 dummyTr.StartTimeIntegral();
e2d2481c 302 dummyTr.AddTimeStep(dst);
c51c3124 303 dummyTr.GetIntegratedTimes(times);
0ddbf657 304 fCurrESDtrack->SetIntegratedTimes(times);
305 fCurrESDtrack->SetIntegratedLength(dummyTr.GetIntegratedLength());
c51c3124 306 continue;
307 }
308 //
309 currTr = GetTrackHyp(itr);
310 currTr->StartTimeIntegral();
311 currTr->AddTimeStep(dst);
15b02d69 312 // printf("Before resetCov: "); currTr->AliExternalTrackParam::Print();
e2d2481c 313 currTr->ResetCovariance(10000);
68a0f687 314 if (RefitTrack(currTr,fITS->GetRMax())) { // propagate to exit from the ITS/TPC screen
315 UpdateESDTrack(currTr,AliESDtrack::kITSout);
0ddbf657 316 fCountITSout++;
68a0f687 317 }
e7d83d38 318 else {
0ddbf657 319 AliDebug(2,Form("Refit Failed for track %d | ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
320 //currTr->AliExternalTrackParam::Print();
321 //currTr->GetWinner()->Print();
e7d83d38 322 }
c51c3124 323 //
324 }
32d38de2 325 //
0ddbf657 326 AliInfo(Form("%d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
327 fCountITSout,fCountITSin,fCountProlongationTrials,nTrESD));
328 //
32d38de2 329 return 0;
330}
331
332//_________________________________________________________________________
e7d83d38 333Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent *esdEv)
32d38de2 334{
335 //
e7d83d38 336 // refit the tracks inward, using current cov.matrix
bdb92766 337 AliCodeTimerAuto("",0);
e7d83d38 338 //
339 SetTrackingPhase(kRefitInw);
340 Int_t nTrESD = esdEv->GetNumberOfTracks();
b515ad7e 341 // AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
342
15b02d69 343 AliDebug(1,Form("Will refit inward %d tracks",nTrESD));
e7d83d38 344 AliITSUTrackHyp *currTr=0;
345 //
346 for (int itr=0;itr<nTrESD;itr++) {
0ddbf657 347 fCurrESDtrack = esdEv->GetTrack(itr);
348 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
e7d83d38 349 // Start time integral and add distance from current position to vertex
0ddbf657 350 UInt_t trStat = fCurrESDtrack->GetStatus();
e7d83d38 351 if ( !(trStat & AliESDtrack::kITSout) ) continue;
352 if ( trStat & AliESDtrack::kITSrefit ) continue; // already done
353 if ( (trStat & AliESDtrack::kTPCout) && !(trStat & AliESDtrack::kTPCrefit) ) continue;
354 //
355 currTr = GetTrackHyp(itr);
0ddbf657 356 currTr->AliExternalTrackParam::operator=(*fCurrESDtrack); // fetch current ESDtrack kinematics
e7d83d38 357 if (RefitTrack(currTr,fITS->GetRMin())) { // propagate up to inside radius of the beam pipe
358 UpdateESDTrack(currTr,AliESDtrack::kITSrefit);
0ddbf657 359 fCountITSrefit++;
e7d83d38 360 }
361 else {
0ddbf657 362 AliDebug(2,Form("Refit Failed for track %d |ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 363 }
364 }
32d38de2 365 //
0ddbf657 366 AliInfo(Form("%d ITSrefit in %d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
367 fCountITSrefit,fCountITSout,fCountITSin,fCountProlongationTrials,nTrESD));
368 //
b515ad7e 369 // AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0);
32d38de2 370 return 0;
371}
372
373//_________________________________________________________________________
374Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
375{
376 // read from tree (if pointer provided) or directly from the ITS reco interface
377 //
378 return fReconstructor->LoadClusters(treeRP);
379}
380
381//_________________________________________________________________________
382void AliITSUTrackerGlo::UnloadClusters()
383{
384 //
385 // To be implemented
386 //
387
388 Info("UnloadClusters","To be implemented");
389}
390//_________________________________________________________________________
391AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
392{
393 //
394 // To be implemented
395 //
396
397 Info("GetCluster","To be implemented");
398 return 0x0;
399}
400
401//_________________________________________________________________________
402Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr)
403{
404 // do we need to match this track to ITS?
405 //
406 static double bz = GetBz();
407 if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
408 esdTr->IsOn(AliESDtrack::kTPCout) ||
409 esdTr->IsOn(AliESDtrack::kITSin) ||
410 esdTr->GetKinkIndex(0)>0) return kFALSE;
411 //
c61e50c3 412 if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
32d38de2 413 //
414 float dtz[2];
415 esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz);
416 // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
c61e50c3 417 if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation())
418 && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
419 Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
32d38de2 420 //
cb50e082 421 // if (esdTr->Pt()<3) return kFALSE;//RS
32d38de2 422 return kTRUE;
423}
424
425//_________________________________________________________________________
716ccba7 426void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
32d38de2 427{
428 // find prolongaion candidates finding for single seed
429 //
3e4e3c23 430 AliITSUSeed seedUC; // copy of the seed from the upper layer
431 AliITSUSeed seedT; // transient seed between the seedUC and new prolongation hypothesis
432 //
6cd80116 433 AliITSUTrackHyp* hypTr = InitHypothesis(esdTr,esdID); // initialize prolongations hypotheses tree
434 fCurrHyp = fWorkHyp;
435 fCurrHyp->InitFrom(hypTr);
32d38de2 436 //
173b3073 437 AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
f8832015 438 //
c61e50c3 439 TObjArray clArr; // container for transfer of clusters matching to seed
32d38de2 440 //
6cd80116 441 for (int ila=fNLrActive;ila--;) {
c03e4f8a 442 fCurrActLrID = ila;
443 fCurrLayer = fITS->GetLayerActive(ila);
444 Bool_t noClSharing = fCurrTrackCond->GetClSharing(ila)==0;
32d38de2 445 int ilaUp = ila+1; // prolong seeds from layer above
3e4e3c23 446 //
447 // for the outermost layer the seed is created from the ESD track
6cd80116 448 int nSeedsUp = (ilaUp==fNLrActive) ? 1 : fCurrHyp->GetNSeeds(ilaUp);
38997c3c 449 int maxNBranches = fCurrTrackCond->GetMaxBranches(ila);
450 int maxNCandidates = fCurrTrackCond->GetMaxCandidates(ila);
3e4e3c23 451 //
32d38de2 452 for (int isd=0;isd<nSeedsUp;isd++) {
3e4e3c23 453 AliITSUSeed* seedU;
6cd80116 454 if (ilaUp==fNLrActive) {
3e4e3c23 455 seedU = 0;
456 seedUC.InitFromESDTrack(esdTr);
457 }
458 else {
459 seedU = fCurrHyp->GetSeed(ilaUp,isd); // seed on prev.active layer to prolong
c03e4f8a 460 if (seedU->IsKilled()) continue;
3e4e3c23 461 seedUC = *seedU; // its copy will be prolonged
b8b59e05 462 seedUC.SetParent(seedU);
3e4e3c23 463 }
943e1898 464 seedUC.ResetFMatrix(); // reset the matrix for propagation to next layer
32d38de2 465 // go till next active layer
0ddbf657 466 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 467 if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) {
32d38de2 468 //
0ddbf657 469 AliDebug(2,Form("Transport failed | esdID=%d (MClb:%d)",esdID,fCurrESDtrMClb));
32d38de2 470 // Check if the seed satisfies to track definition
b8b59e05 471 if (NeedToKill(&seedUC,kTransportFailed) && seedU) KillSeed(seedU,kTRUE);
32d38de2 472 continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
473 }
f8832015 474 if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
b8b59e05 475 if (NeedToKill(&seedUC,kRWCheckFailed) && seedU) KillSeed(seedU,kTRUE);
32d38de2 476 continue;
477 }
5785a9d8 478 /*
479 //RS toremove
480 int mcquest = -1;
c03e4f8a 481 if (!seedUC.ContainsFake() && AliDebugLevelClass()>2) {
5785a9d8 482 mcquest = fCurrESDtrMClb;
483 seedUC.Print("etp");
484 printf("FSParams: "); for (int ip=0;ip<kNTrImpData;ip++) printf("%+e ",fTrImpData[ip]); printf("\n");
485 }
486 //
c03e4f8a 487 int nsens = fCurrLayer->FindSensors(&fTrImpData[kTrPhi0], hitSens, mcquest); // find detectors which may be hit by the track
5785a9d8 488 */
c03e4f8a 489 int nsens = fCurrLayer->FindSensors(&fTrImpData[kTrPhi0], hitSens); // find detectors which may be hit by the track
0ddbf657 490 AliDebug(2,Form("Will check %d sensors on lr:%d | esdID=%d (MClb:%d)",nsens,ila,esdID,fCurrESDtrMClb));
32d38de2 491 //
38997c3c 492 seedUC.SetLr(ila);
493 //
943e1898 494 double bz = GetBz();
32d38de2 495 for (int isn=nsens;isn--;) {
c61e50c3 496 AliITSURecoSens* sens = hitSens[isn];
38997c3c 497 int ncl = sens->GetNClusters();
498 if (!ncl) continue;
499 seedT = seedUC;
f8832015 500 //
c03e4f8a 501 // We need to propagate the seed to sensor on fCurrLayer staying in the frame of the sensor from prev.layer,
943e1898 502 // since the transport matrix should be defined in this frame.
503 double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
53870004 504 if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) {
505 if (AliDebugLevelClass()>2) {
506 printf("Failed on GetTrackingXAtXAlpha: X=%.4f alp=%.3f\n",sens->GetXTF(),sens->GetPhiTF());
507 seedT.Print("etp");
508 }
509 continue;
510 }
511 if (xs<seedT.GetX()) {
512 if (!PropagateSeed(&seedT,xs,fCurrMass)) continue;
513 }
514 else { // some low precision tracks may hit the sensor plane outside of the layer radius
515 if (AliDebugLevelClass()>2) {
516 if (!seedT.ContainsFake()) {
517 printf("WRONG PROP on L%d, sens%d of %d: %.4f > %.4f\n",ila,isn,nsens,xs,seedT.GetX());
518 sens->Print();
519 seedT.Print("etp");
520 }
521 }
522 if (!seedT.PropagateParamOnlyTo(xs,bz)) continue;
523 }
8b16dbae 524 // if (!seedT.PropagateToX(xs,bz)) continue;
716ccba7 525 // if (!seedT.Rotate(sens->GetPhiTF())) continue;
526 if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
943e1898 527 //
f8832015 528 int clID0 = sens->GetFirstClusterId();
cb50e082 529 for (int icl=ncl;icl--;) { // don't change the order, clusters are sorted
c03e4f8a 530 //AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
531 int clID = clID0 + icl;
532 if (noClSharing && fCurrLayer->GetCluster(clID)->IsClusterUsed()) continue;
f8832015 533 int res = CheckCluster(&seedT,ila,clID0+icl);
c03e4f8a 534 //AliLog::SetClassDebugLevel("AliITSUTrackerGlo", 0);
f8832015 535 //
536 if (res==kStopSearchOnSensor) break; // stop looking on this sensor
537 if (res==kClusterNotMatching) continue; // cluster does not match
538 // cluster is matching and it was added to the hypotheses tree
c61e50c3 539 }
32d38de2 540 }
38997c3c 541 // cluster search is done. Do we need to have a version of this seed skipping current layer
542 if (!NeedToKill(&seedUC,kMissingCluster)) {
b2d3bc3a 543 AliITSUSeed* seedSkp = NewSeedFromPool(&seedUC);
ee54014a 544 double penalty = -fCurrTrackCond->GetMissPenalty(ila);
716ccba7 545 // to do: make penalty to account for probability to miss the cluster for good reason
546 seedSkp->SetChi2Cl(penalty);
c03e4f8a 547 AddSeedBranch(seedSkp);
548 // AddProlongationHypothesis(seedSkp,ila);
716ccba7 549 }
38997c3c 550 // transfer the new branches of the seed to the hypothesis container
551 if (fNBranchesAdded) ValidateAllowedBranches(maxNBranches);
552 //
32d38de2 553 }
38997c3c 554 if (fNCandidatesAdded) ValidateAllowedCandidates(ila,maxNCandidates);
555 // ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
556 if (AliDebugLevelClass()>2) { //RS
15b02d69 557 printf(">>> All hypotheses on lr %d: \n",ila);
558 for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();}
559 }
560 //
561 /*
716ccba7 562 for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
716ccba7 563 if (ila!=0) continue;
564 double vecL[5] = {0};
565 double matL[15] = {0};
566 AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
567 while(sp->GetParent()) {
568 sp->Smooth(vecL,matL);
6cd80116 569 if (sp->GetLayerID()>=fNLrActive-1) break;
716ccba7 570 sp = (AliITSUSeed*)sp->GetParent();
571 }
572 */
32d38de2 573 }
574 //
6cd80116 575 SaveReducedHypothesesTree(hypTr);
576 if (AliDebugLevelClass()>1) {
577 printf("\nSaved hypotheses for esdTrack %d (MCLab:%d)\n",esdID,fCurrESDtrMClb);
578 hypTr->Print("l");
579 }
580 fCurrHyp = 0;
716ccba7 581 //
32d38de2 582}
583
584//_________________________________________________________________________
6cd80116 585AliITSUTrackHyp* AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esdID)
32d38de2 586{
c61e50c3 587 // init prolongaion candidates finding for single seed
6cd80116 588 AliITSUTrackHyp* hyp = GetTrackHyp(esdID);
589 if (hyp) return hyp;
716ccba7 590 //
0ddbf657 591 fCountProlongationTrials++;
592 //
32d38de2 593 fCurrMass = esdTr->GetMass();
594 if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
c61e50c3 595 //
6cd80116 596 hyp = new AliITSUTrackHyp(fNLrActive);
597 hyp->SetESDTrack(esdTr);
598 hyp->SetUniqueID(esdID);
599 hyp->SetMass(fCurrMass);
600 SetTrackHyp(hyp,esdID);
3e4e3c23 601 //
6cd80116 602 return hyp;
32d38de2 603 // TO DO
604}
605
606//_________________________________________________________________________
607Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo)
608{
0ddbf657 609 // transport seed from layerFrom to the entrance (along the direction of the track) of layerTo
32d38de2 610 //
8b16dbae 611 //
612 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
613 //
32d38de2 614 int dir = lTo > lFrom ? 1:-1;
615 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
616 Bool_t checkFirst = kTRUE;
617 while(lFrom!=lTo) {
32d38de2 618 if (lrFr) {
b2d3bc3a 619 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) {
0ddbf657 620 //printf("FailHere0 Dir=%d\n",dir);
621 //seed->Print("etp");
b2d3bc3a 622 return kFALSE; // go till the end of current layer
623 }
e2d2481c 624 checkFirst = kFALSE;
32d38de2 625 }
626 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
627 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
628 //
629 // go the entrance of the layer, assuming no materials in between
8b16dbae 630 double xToGo = lrTo->GetR(-dir);
0ddbf657 631 // double xts = xToGo;
b2d3bc3a 632 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
0ddbf657 633 //printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
634 //seed->Print("etp");
635 //
636
b2d3bc3a 637 return kFALSE;
638 }
b515ad7e 639 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
b2d3bc3a 640 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
0ddbf657 641 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
642 //seed->Print("etp");
b2d3bc3a 643 return kFALSE;
644 }
32d38de2 645 lrFr = lrTo;
646 }
647 return kTRUE;
648 //
649}
650
3e4e3c23 651//_________________________________________________________________________
652Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo)
653{
654 // transport track from layerFrom to the entrance of layerTo
655 //
8b16dbae 656 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
657 //
3e4e3c23 658 int dir = lTo > lFrom ? 1:-1;
659 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
660 Bool_t checkFirst = kTRUE;
661 while(lFrom!=lTo) {
3e4e3c23 662 if (lrFr) {
e2d2481c 663 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
664 checkFirst = kFALSE;
3e4e3c23 665 }
666 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
667 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
668 //
0ddbf657 669 // go the entrance of the layer, assuming no materials in between
670 double xToGo = lrTo->GetR(-dir);
671 // double xts = xToGo;
672 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
673 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
674 // seed->Print("etp");
675 return kFALSE;
676 }
b515ad7e 677 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
678 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
679 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
680 //seed->Print("etp");
681 return kFALSE;
682 }
683 lrFr = lrTo;
684 }
685 return kTRUE;
686 //
687}
688
689//_________________________________________________________________________
690Bool_t AliITSUTrackerGlo::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
691{
692 // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
693 //
694 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
695 //
696 int dir = lTo > lFrom ? 1:-1;
697 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
698 Bool_t checkFirst = kTRUE;
699 while(lFrom!=lTo) {
700 if (lrFr) {
701 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
702 checkFirst = kFALSE;
703 }
704 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
705 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
706 //
707 // go the entrance of the layer, assuming no materials in between
708 double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
709 //
710 // double xts = xToGo;
711 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
712 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
713 // seed->Print("etp");
714 return kFALSE;
715 }
716 if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
717 //
718 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
0ddbf657 719 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
720 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
721 //seed->Print("etp");
722 return kFALSE;
723 }
3e4e3c23 724 lrFr = lrTo;
725 }
726 return kTRUE;
727 //
728}
729
e2d2481c 730//_________________________________________________________________________
731Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
732{
733 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
734 // If check is requested, do this only provided the track has not exited the layer already
735 double xToGo = lr->GetR(dir);
736 if (check) { // do we need to track till the surface of the current layer ?
737 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
15b02d69 738 AliDebug(3,Form(" dir:%d Cur: %e Tgt: %e",dir,Sqrt(curR2),xToGo));
0ddbf657 739 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
740 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 741 }
742 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
743 // go via layer to its boundary, applying material correction.
744 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 745 //
e2d2481c 746 return kTRUE;
747 //
748}
749
750//_________________________________________________________________________
751Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
752{
753 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
754 // If check is requested, do this only provided the track has not exited the layer already
755 double xToGo = lr->GetR(dir);
756 if (check) { // do we need to track till the surface of the current layer ?
757 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 758 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
759 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 760 }
b515ad7e 761 AliDebug(2,Form(" dir=%d : from R=%.4f -> R=%.4f",dir,Sqrt(seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY()), xToGo));
e2d2481c 762 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
763 // go via layer to its boundary, applying material correction.
764 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 765 //
e2d2481c 766 return kTRUE;
767 //
768}
769
770
771//_________________________________________________________________________
772Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
773{
774 // go to the entrance of lr in direction dir, w/o applying material corrections.
775 // If check is requested, do this only provided the track did not reach the layer already
776 double xToGo = lr->GetR(-dir);
777 if (check) { // do we need to track till the surface of the current layer ?
778 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 779 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
780 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 781 }
782 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
783 // go via layer to its boundary, applying material correction.
784 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
785 return kTRUE;
786 //
787}
788
789//_________________________________________________________________________
790Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
791{
792 // go to the entrance of lr in direction dir, w/o applying material corrections.
793 // If check is requested, do this only provided the track did not reach the layer already
794 double xToGo = lr->GetR(-dir);
795 if (check) { // do we need to track till the surface of the current layer ?
796 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 797 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
798 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 799 }
800 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
801 // go via layer to its boundary, applying material correction.
802 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
803 return kTRUE;
804 //
805}
806
32d38de2 807//_________________________________________________________________________
c61e50c3 808Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
32d38de2 809{
810 // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
811 // as well as some aux info
812 double bz = GetBz();
813 AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
c61e50c3 814 seed->GetXYZ(&fTrImpData[kTrXIn]); // lab position at the entrance from above
53870004 815 static AliExternalTrackParam sc; // seed copy for manipulations
943e1898 816 sc = *seed;
32d38de2 817 //
c61e50c3 818 fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
943e1898 819 if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
c61e50c3 820 double dr = lrA->GetDR(); // approximate X dist at the inner radius
943e1898 821 if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
32d38de2 822 // special case: track does not reach inner radius, might be tangential
943e1898 823 double r = sc.GetD(0,0,bz);
32d38de2 824 double x;
bdb92766 825 if (!sc.GetXatLabR(r,x,bz,-1)) return kFALSE;
943e1898 826 dr = Abs(sc.GetX() - x);
bdb92766 827 if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) return kFALSE;
32d38de2 828 }
829 //
c61e50c3 830 fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
943e1898 831 double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
832 double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
ee54014a 833 sgy = Sqrt(sgy)*fCurrTrackCond->GetNSigmaRoadY(ilrA);
834 sgz = Sqrt(sgz)*fCurrTrackCond->GetNSigmaRoadZ(ilrA);
5785a9d8 835 double dphi0 = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]);
836 double phi0 = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
837 if ( dphi0>(0.5*Pi()) ) {
838 // special case of angles around pi
839 dphi0 = Abs(phi0);
840 phi0 += Pi();
841 }
842
843 fTrImpData[kTrPhi0] = phi0;
173b3073 844 fTrImpData[kTrZ0] = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
5785a9d8 845 fTrImpData[kTrDPhi] = dphi0 + sgy/lrA->GetR();
173b3073 846 fTrImpData[kTrDZ] = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn]) + sgz;
32d38de2 847 //
848 return kTRUE;
849}
850
b8b59e05 851//________________________________________
852void AliITSUTrackerGlo::ResetSeedsPool()
853{
854 // mark all seeds in the pool as unused
0ddbf657 855 AliDebug(-1,Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
b8b59e05 856 fNFreeSeeds = 0;
857 fSeedsPool.Clear(); // seeds don't allocate memory
858}
859
860
861//________________________________________
862void AliITSUTrackerGlo::MarkSeedFree(AliITSUSeed *sd)
32d38de2 863{
b8b59e05 864 // account that this seed is "deleted"
865 int id = sd->GetPoolID();
866 if (id<0) {
867 AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
868 return;
869 }
870 // AliInfo(Form("%d %p",id, seed));
871 fSeedsPool.RemoveAt(id);
872 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
873 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
32d38de2 874}
f8832015 875
876//_________________________________________________________________________
877Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
878{
879 // Check if the cluster (in tracking frame!) is matching to track.
880 // The track must be already propagated to sensor tracking frame.
881 // Returns: kStopSearchOnSensor if the search on given sensor should be stopped,
882 // kClusterMatching if the cluster is matching
883 // kClusterMatching otherwise
884 //
f8832015 885 const double kTolerX = 5e-4;
c03e4f8a 886 // AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
887 AliCluster *cl = fCurrLayer->GetCluster(clID);
f8832015 888 //
c8d1f258 889 Bool_t goodCl = kFALSE;
0ddbf657 890 int currLabel = Abs(fCurrESDtrMClb);
716ccba7 891 //
8b16dbae 892 if (cl->GetLabel(0)>=0) {
893 for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}
894 }
53870004 895 Bool_t goodSeed = !track->ContainsFake();
c8d1f258 896 //
f8832015 897 if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
c8d1f258 898 if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
53870004 899 if (goodCl&&goodSeed && AliDebugLevelClass()>2 ) {
900 AliDebug(2,Form("Lost good cl on L:%d failed propagation. |ESDtrack#%d (MClb:%d)",lr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
901 track->Print("etp");
15b02d69 902 cl->Print();
903 }
c8d1f258 904 return kStopSearchOnSensor; // propagation failed, seedT is intact
905 }
f8832015 906 }
3dd9c283 907 double dy = cl->GetY()-track->GetY();
c8d1f258 908 double dz = cl->GetZ()-track->GetZ();
0091e9f0 909 //
76390254 910#ifdef _FILL_CONTROL_HISTOS_
bdb92766 911 int hcOffs = (1+fTrackPhase)*kHistosPhase + lr;
76390254 912 double htrPt=-1;
cb50e082 913 if (goodCl && (((AliITSUClusterPix*)cl)->GetNPix()>1 || !((AliITSUClusterPix*)cl)->IsSplit()) && goodSeed && fCHistoArr /* && track->GetChi2Penalty()<1e-5*/) {
76390254 914 htrPt = track->Pt();
53870004 915 ((TH2*)fCHistoArr->At(kHResY+hcOffs))->Fill(htrPt,dy);
916 ((TH2*)fCHistoArr->At(kHResZ+hcOffs))->Fill(htrPt,dz);
76390254 917 double errY = track->GetSigmaY2();
918 double errZ = track->GetSigmaZ2();
53870004 919 if (errY>0) ((TH2*)fCHistoArr->At(kHResYP+hcOffs))->Fill(htrPt,dy/Sqrt(errY));
920 if (errZ>0) ((TH2*)fCHistoArr->At(kHResZP+hcOffs))->Fill(htrPt,dz/Sqrt(errZ));
76390254 921 }
922#endif
923 //
f8832015 924 double dy2 = dy*dy;
925 double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
ee54014a 926 fCurrTrackCond->GetNSigmaRoadY(lr)*fCurrTrackCond->GetNSigmaRoadY(lr); // RS TOOPTIMIZE
f8832015 927 if (dy2>tol2) { // the clusters are sorted in Z(col) then in Y(row).
53870004 928 if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
929 AliDebug(2,Form("Lost good cl: dy2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",dy2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
930 track->Print("etp");
15b02d69 931 cl->Print();
53870004 932 PrintSeedClusters(track);
15b02d69 933 }
cb50e082 934 // clusters are sorted in increasing Y and the loop goes from last (highers Y) to 1st.
935 // if the track has too large y for this cluster (dy<0) it will be so for all other clusters of the sensor
936 if (dy<0) return kStopSearchOnSensor; // No chance that other cluster of this sensor will match (all Y's will be even larger)
f8832015 937 else return kClusterNotMatching; // Other clusters may match
938 }
c8d1f258 939 double dz2 = dz*dz;
f8832015 940 tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
ee54014a 941 fCurrTrackCond->GetNSigmaRoadZ(lr)*fCurrTrackCond->GetNSigmaRoadZ(lr); // RS TOOPTIMIZE
c8d1f258 942 if (dz2>tol2) {
53870004 943 if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
944 AliDebug(2,Form("Lost good cl on L:%d : dz2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",lr,dz2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
945 track->Print("etp");
15b02d69 946 cl->Print();
53870004 947 PrintSeedClusters(track);
15b02d69 948 }
c8d1f258 949 return kClusterNotMatching; // Other clusters may match
950 }
f8832015 951 //
952 // check chi2
953 Double_t p[2]={cl->GetY(), cl->GetZ()};
954 Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
955 double chi2 = track->GetPredictedChi2(p,cov);
76390254 956 //
957#ifdef _FILL_CONTROL_HISTOS_
958 if (htrPt>0) {
53870004 959 ((TH2*)fCHistoArr->At(kHChi2Cl+hcOffs))->Fill(htrPt,chi2);
76390254 960 }
961#endif
962 //
ee54014a 963 if (chi2>fCurrTrackCond->GetMaxTr2ClChi2(lr)) {
53870004 964 if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
965 AliDebug(2,Form("Lost good cl on L:%d : Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
ee54014a 966 lr,chi2,fCurrTrackCond->GetMaxTr2ClChi2(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb));
3dd9c283 967 track->Print("etp");
968 cl->Print("");
53870004 969 PrintSeedClusters(track);
3dd9c283 970 }
c8d1f258 971 return kClusterNotMatching;
972 }
f8832015 973 //
974 track = NewSeedFromPool(track); // input track will be reused, use its clone for updates
943e1898 975 if (!track->Update()) {
53870004 976 if (goodCl&&goodSeed && AliDebugLevelClass()>2) {
977 AliDebug(2,Form("Lost good cluster on L:%d : failed to update",lr));
38997c3c 978 track->Print("etp");
979 cl->Print("");
53870004 980 PrintSeedClusters(track);
38997c3c 981 }
b8b59e05 982 MarkSeedFree(track);
c8d1f258 983 return kClusterNotMatching;
984 }
f8832015 985 track->SetChi2Cl(chi2);
986 track->SetLrClusterID(lr,clID);
c03e4f8a 987 // cl->IncreaseClusterUsage(); // do this only for winners
f8832015 988 //
716ccba7 989 track->SetFake(!goodCl);
990 //
0ddbf657 991 AliDebug(2,Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)| ",
716ccba7 992 goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
c8d1f258 993 //
38997c3c 994 AddSeedBranch(track);
c03e4f8a 995#ifdef _FILL_CONTROL_HISTOS_
996 if (htrPt>0) {
997 ((TH2*)fCHistoArr->At(kHChi2Nrm+hcOffs))->Fill(htrPt,track->GetChi2GloNrm());
998 }
999#endif
f8832015 1000 //
1001 return kClusterMatching;
1002}
c8d1f258 1003
1004//_________________________________________________________________________
1005Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
1006{
1007 // check if the seed should not be discarded
1008 const UShort_t kMask = 0xffff;
1009 if (flag==kMissingCluster) {
1010 int lastChecked = seed->GetLayerID();
1011 UShort_t patt = seed->GetHitsPattern();
c03e4f8a 1012 if (lastChecked) patt |= ~(kMask<<lastChecked); // not all layers were checked, complete unchecked ones by potential hits
42c3d4bd 1013 Bool_t seedOK = fCurrTrackCond->CheckPattern(patt);
c8d1f258 1014 return !seedOK;
1015 }
1016 return kTRUE;
1017}
44785f3e 1018
1019//______________________________________________________________________________
1020Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
1021{
1022 // propagate seed to given x applying material correction if requested
1023 const Double_t kEpsilon = 1e-5;
1024 Double_t xpos = seed->GetX();
1025 Int_t dir = (xpos<xToGo) ? 1:-1;
1026 Double_t xyz0[3],xyz1[3],param[7];
1027 //
8b16dbae 1028 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
f497470c 1029 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
44785f3e 1030 while ( (xToGo-xpos)*dir > kEpsilon){
1031 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1032 Double_t x = xpos+step;
1033 Double_t bz=GetBz(); // getting the local Bz
0ddbf657 1034 if (!seed->PropagateToX(x,bz)) return kFALSE;
8b16dbae 1035 double ds = 0;
1036 if (matCorr || updTime) {
44785f3e 1037 xyz0[0]=xyz1[0]; // global pos at the beginning of step
1038 xyz0[1]=xyz1[1];
1039 xyz0[2]=xyz1[2];
1040 seed->GetXYZ(xyz1); // // global pos at the end of step
8b16dbae 1041 if (matCorr) {
1042 MeanMaterialBudget(xyz0,xyz1,param);
1043 Double_t xrho=param[0]*param[4], xx0=param[1];
1044 if (dir>0) xrho = -xrho; // outward should be negative
b2d3bc3a 1045 if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) {AliInfo("Fail2"); seed->Print("etp"); return kFALSE;}
8b16dbae 1046 ds = param[4];
1047 }
1048 else { // matCorr is not requested but time integral is
1049 double d0 = xyz1[0]-xyz0[0];
1050 double d1 = xyz1[1]-xyz0[1];
1051 double d2 = xyz1[2]-xyz0[2];
1052 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1053 }
44785f3e 1054 }
8b16dbae 1055 if (updTime) seed->AddTimeStep(ds);
44785f3e 1056 xpos = seed->GetX();
1057 }
1058 return kTRUE;
1059}
716ccba7 1060
3e4e3c23 1061//______________________________________________________________________________
1062Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
1063{
1064 // propagate seed to given x applying material correction if requested
1065 const Double_t kEpsilon = 1e-5;
1066 Double_t xpos = seed->GetX();
1067 Int_t dir = (xpos<xToGo) ? 1:-1;
1068 Double_t xyz0[3],xyz1[3],param[7];
1069 //
b515ad7e 1070 if (AliDebugLevelClass()>1) {
1071 AliDebug(2,Form("Before Propagate to X=%f with M=%.3f MaxStep=%.4f MatCorr=%d",xToGo,mass,maxStep,matCorr));
1072 seed->AliExternalTrackParam::Print();
1073 }
3e4e3c23 1074 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
1075 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
1076 while ( (xToGo-xpos)*dir > kEpsilon){
1077 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1078 Double_t x = xpos+step;
1079 Double_t bz=GetBz(); // getting the local Bz
1080 if (!seed->PropagateTo(x,bz)) return kFALSE;
1081 double ds = 0;
1082 if (matCorr || updTime) {
1083 xyz0[0]=xyz1[0]; // global pos at the beginning of step
1084 xyz0[1]=xyz1[1];
1085 xyz0[2]=xyz1[2];
1086 seed->GetXYZ(xyz1); // // global pos at the end of step
1087 //
1088 if (matCorr) {
1089 MeanMaterialBudget(xyz0,xyz1,param);
1090 Double_t xrho=param[0]*param[4], xx0=param[1];
1091 if (dir>0) xrho = -xrho; // outward should be negative
1092 if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
1093 ds = param[4];
1094 }
1095 else { // matCorr is not requested but time integral is
1096 double d0 = xyz1[0]-xyz0[0];
1097 double d1 = xyz1[1]-xyz0[1];
1098 double d2 = xyz1[2]-xyz0[2];
1099 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1100 }
1101 }
1102 if (updTime) seed->AddTimeStep(ds);
1103 //
1104 xpos = seed->GetX();
1105 }
b515ad7e 1106 if (AliDebugLevelClass()>1) {
1107 AliDebug(2,Form("After Propagate to X=%f",xToGo));
1108 seed->AliExternalTrackParam::Print();
1109 }
3e4e3c23 1110 return kTRUE;
1111}
1112
3e4e3c23 1113//______________________________________________________________________________
1114void AliITSUTrackerGlo::FinalizeHypotheses()
1115{
1116 // select winner for each hypothesis, remove cl. sharing conflicts
15b02d69 1117 AliDebug(2,"TODO");
3e4e3c23 1118 //
1119 int nh = fHypStore.GetEntriesFast();
c03e4f8a 1120 AliITSUSeed* winner = 0;
3e4e3c23 1121 for (int ih=0;ih<nh;ih++) {
1122 AliITSUTrackHyp* hyp = (AliITSUTrackHyp*) fHypStore.UncheckedAt(ih);
c03e4f8a 1123 if (!hyp) continue; // removed
1124 if (hyp->GetESDTrack()->IsOn(AliESDtrack::kITSin)) continue; // winner already defined
1125 //
1126 if (!(winner=hyp->DefineWinner())) continue; // failed to find winner
1127 FlagSeedClusters(winner,kTRUE);
0ddbf657 1128 fCountITSin++;
68a0f687 1129 CookMCLabel(hyp);
1130 UpdateESDTrack(hyp,AliESDtrack::kITSin);
3e4e3c23 1131 }
c03e4f8a 1132 //
bdb92766 1133#ifdef _FILL_CONTROL_HISTOS_
1134 // if requested, collect cluster sharing statistics
1135 TH2* hShare = 0;
1136 if (fCHistoArr && (hShare=(TH2*)fCHistoArr->At(kHClShare))) {
1137 for (int ih=0;ih<nh;ih++) {
1138 AliITSUTrackHyp* hyp = (AliITSUTrackHyp*) fHypStore.UncheckedAt(ih);
1139 if (!hyp || !(winner=hyp->GetWinner())) continue;
1140 int lrID = 0;
1141 double pt = hyp->Pt();
1142 do {
1143 int clID = winner->GetLrCluster(lrID);
1144 if (clID<0) continue;
1145 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1146 if (!cl->IsClusterShared()) continue;
1147 hShare->Fill(pt,winner->IsFake() ? lrID+fNLrActive : lrID);
1148 } while ((winner=(AliITSUSeed*)winner->GetParent()));
1149 }
1150 }
1151#endif
1152
3e4e3c23 1153}
c51c3124 1154
1155//______________________________________________________________________________
68a0f687 1156void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp,Int_t flag)
c51c3124 1157{
1158 // update ESD track with current best hypothesis
1159 AliESDtrack* esdTr = hyp->GetESDTrack();
1160 if (!esdTr) return;
1161 AliITSUSeed* win = hyp->GetWinner();
1162 if (!win) return;
68a0f687 1163 switch (flag) {
1164 case AliESDtrack::kITSin:
1165 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1166 // TODO: set cluster info
1167 break;
1168 //
1169 case AliESDtrack::kITSout:
1170 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1171 // TODO: avoid setting friend
1172 break;
1173 //
1174 case AliESDtrack::kITSrefit:
1175 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
1176 // TODO: avoid setting cluster info
1177 break;
1178 default:
1179 AliFatal(Form("Unknown flag %d",flag));
1180 }
c51c3124 1181 //
bdb92766 1182 esdTr->SetITSLabel(hyp->GetITSLabel());
c51c3124 1183 // transfer module indices
1184 // TODO
1185}
1186
1187//______________________________________________________________________________
68a0f687 1188Bool_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest, Bool_t stopAtLastCl)
c51c3124 1189{
1190 // refit track till radius rDest
8b16dbae 1191 AliITSUTrackHyp tmpTr;
1192 //
c51c3124 1193 double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
1194 int dir,lrStart,lrStop;
1195 //
8b16dbae 1196 dir = rCurr<rDest ? 1 : -1;
1197 lrStart = fITS->FindFirstLayerID(rCurr,dir);
1198 lrStop = fITS->FindLastLayerID(rDest,dir);
15b02d69 1199 if (AliDebugLevelClass()>2) {
1200 printf("Refit %d: Lr:%d (%f) -> Lr:%d (%f)\n",dir,lrStart,rCurr, lrStop,rDest);
1201 printf("Before refit: "); trc->AliExternalTrackParam::Print();
1202 }
8b16dbae 1203 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));
1204 //
68a0f687 1205 int ncl = trc->FetchClusterInfo(fClInfo);
8b16dbae 1206 fCurrMass = trc->GetMass();
1207 tmpTr.AliKalmanTrack::operator=(*trc);
1208 tmpTr.SetChi2(0);
68a0f687 1209 int iclLr[2],nclLr,clCount=0;
c51c3124 1210 //
1211 for (int ilr=lrStart;ilr!=lrStop;ilr+=dir) {
1212 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
1213 if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
8b16dbae 1214 int ilrA2,ilrA = lr->GetActiveID();
1215 // passive layer or active w/o hits will be traversed on the way to next cluster
1216 if (!lr->IsActive() || fClInfo[ilrA2=(ilrA<<1)]<0) continue;
1217 //
8b16dbae 1218 // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
1219 nclLr=0;
1220 if (dir>0) { // clusters are stored in increasing radius order
1221 iclLr[nclLr++]=fClInfo[ilrA2++];
1222 if (fClInfo[ilrA2]>=0) iclLr[nclLr++]=fClInfo[ilrA2];
1223 }
1224 else {
1225 if ( fClInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=fClInfo[ilrA2+1];
1226 iclLr[nclLr++]=fClInfo[ilrA2];
1227 }
1228 //
b515ad7e 1229 Bool_t transportedToLayer = kFALSE;
8b16dbae 1230 for (int icl=0;icl<nclLr;icl++) {
1231 AliITSUClusterPix* clus = (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
1232 AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
e7d83d38 1233 if (!tmpTr.Rotate(sens->GetPhiTF())) {
0ddbf657 1234 AliDebug(2,Form("Failed on rotate to %f | ESDtrack#%d (MClb:%d)",sens->GetPhiTF(),fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 1235 return kFALSE;
1236 }
b515ad7e 1237 //
1238 double xClus = sens->GetXTF()+clus->GetX();
1239 if (!transportedToLayer) {
1240 if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
1241 AliDebug(2,Form("Failed to transport %d -> %d | ESDtrack#%d (MClb:%d)\n",lrStart,ilr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1242 //tmpTr.AliExternalTrackParam::Print(); trc->GetWinner()->Print("etp");
1243 return kFALSE; // go to the entrance to the layer
1244 }
1245 lrStart = ilr;
1246 transportedToLayer = kTRUE;
1247 }
1248 //
1249 if (AliDebugLevelClass()>1) {
1250 AliDebug(2,Form("Propagate to cl:%d on lr %d Need to go %.4f -> %.4f",icl,ilrA,tmpTr.GetX(),xClus));
1251 // printf("Before: "); tmpTr.AliExternalTrackParam::Print();
1252 }
1253 //
1254 if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) {
1255 AliDebug(2,Form("Failed on propagate to %f (dir=%d) | ESDtrack#%d (MClb:%d)",xClus,dir,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1256 //tmpTr.AliExternalTrackParam::Print(""); trc->GetWinner()->Print("etp");
e7d83d38 1257 return kFALSE;
1258 }
76390254 1259 //
1260#ifdef _FILL_CONTROL_HISTOS_
bdb92766 1261 int hcOffs = (1+fTrackPhase)*kHistosPhase + ilrA;
76390254 1262 double htrPt=-1;
53870004 1263 if (fCHistoArr && trc->GetLabel()>=0/* && trc->Charge()>0*/) {
76390254 1264 htrPt = tmpTr.Pt();
1265 double dy = clus->GetY()-tmpTr.GetY();
1266 double dz = clus->GetZ()-tmpTr.GetZ();
53870004 1267 ((TH2*)fCHistoArr->At(kHResY+hcOffs))->Fill(htrPt,dy);
1268 ((TH2*)fCHistoArr->At(kHResZ+hcOffs))->Fill(htrPt,dz);
76390254 1269 double errY = tmpTr.GetSigmaY2();
1270 double errZ = tmpTr.GetSigmaZ2();
53870004 1271 if (errY>0) ((TH2*)fCHistoArr->At(kHResYP+hcOffs))->Fill(htrPt,dy/Sqrt(errY));
1272 if (errZ>0) ((TH2*)fCHistoArr->At(kHResZP+hcOffs))->Fill(htrPt,dz/Sqrt(errZ));
76390254 1273 }
1274#endif
1275 //
1276 double chi2;
1277 if ( (chi2=tmpTr.Update(clus))<0 ) {
0ddbf657 1278 AliDebug(2,Form("Failed on Update | ESDtrack#%d (MClb:%d)",fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 1279 return kFALSE;
1280 }
76390254 1281#ifdef _FILL_CONTROL_HISTOS_
1282 if (htrPt>0) {
53870004 1283 ((TH2*)fCHistoArr->At(kHChi2Cl+hcOffs))->Fill(htrPt,chi2);
76390254 1284 }
1285#endif
b515ad7e 1286 if (AliDebugLevelClass()>1) {
1287 printf("AfterRefit: "); tmpTr.AliExternalTrackParam::Print();
1288 }
68a0f687 1289 if (stopAtLastCl && ++clCount==ncl) return kTRUE; // it was requested to not propagate after last update
8b16dbae 1290 }
c51c3124 1291 //
1292 }
8b16dbae 1293 // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
1294 // Still, try to go as close as possible to rDest.
1295 //
1296 if (lrStart!=lrStop) {
e7d83d38 1297 //printf("Going to last layer %d -> %d\n",lrStart,lrStop);
1298 if (!TransportToLayer(&tmpTr,lrStart,lrStop)) {
0ddbf657 1299 AliDebug(1,Form("Failed on TransportToLayer %d->%d | ESDtrack#%d (MClb:%d)",lrStart,lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1300 //tmpTr.AliExternalTrackParam::Print();
1301 //trc->GetWinner()->Print("etp");
e7d83d38 1302 return kTRUE;
1303 }
1304 if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) {
0ddbf657 1305 AliDebug(3,Form("Failed on GoToExitFromLayer %d | ESDtrack#%d (MClb:%d)",lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb));
e7d83d38 1306 return kTRUE; // go till the exit from layer
1307 }
e2d2481c 1308 //
e7d83d38 1309 //printf("On exit from last layer\n");
15b02d69 1310 //tmpTr.AliExternalTrackParam::Print();
e2d2481c 1311 // go to the destination radius
1312 if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),dir)) return kTRUE;
1313 if (!PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) return kTRUE;
8b16dbae 1314 }
1315 trc->AliKalmanTrack::operator=(tmpTr);
15b02d69 1316 if (AliDebugLevelClass()>2) {
b515ad7e 1317 printf("After refit (now at lr %d): ",lrStop); trc->AliExternalTrackParam::Print();
15b02d69 1318 }
8b16dbae 1319 return kTRUE;
c51c3124 1320}
dd2117a2 1321
1322//__________________________________________________________________
1323void AliITSUTrackerGlo::CookMCLabel(AliITSUTrackHyp* hyp)
1324{
1325 // build MC label
1326 //
1327 const int kMaxLbPerCl = 3;
1328 int lbID[kMaxLayers*6],lbStat[kMaxLayers*6];
1329 Int_t lr,nLab=0,nCl=0;
1330 AliITSUSeed *seed = hyp->GetWinner();
1331 while(seed) {
1332 int clID = seed->GetLrCluster(lr);
1333 if (clID>=0) {
1334 AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1335 nCl++;
1336 for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
1337 int trLb = cl->GetLabel(imc);
704d9398 1338 if (trLb<0) break;
dd2117a2 1339 // search this mc track in already accounted ones
1340 int iLab;
1341 for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
1342 if (iLab<nLab) lbStat[iLab]++;
1343 else {
1344 lbID[nLab] = trLb;
1345 lbStat[nLab++] = 1;
1346 }
1347 } // loop over given cluster's labels
1348 }
1349 seed = (AliITSUSeed*)seed->GetParent();
1350 } // loop over clusters
1351 //
5785a9d8 1352 AliESDtrack* esdTr = hyp->GetESDTrack();
1353 int tpcLab = esdTr ? Abs(esdTr->GetTPCLabel()) : -kDummyLabel;
41343470 1354 if (nCl && nLab) {
dd2117a2 1355 int maxLab=0,nTPCok=0;
dd2117a2 1356 for (int ilb=nLab;ilb--;) {
1357 int st = lbStat[ilb];
1358 if (lbStat[maxLab]<st) maxLab=ilb;
1359 if (tpcLab==lbID[ilb]) nTPCok += st;
1360 }
1361 hyp->SetFakeRatio(1.-float(nTPCok)/nCl);
1362 hyp->SetLabel( nTPCok==nCl ? tpcLab : -tpcLab);
bdb92766 1363 hyp->SetITSLabel( lbStat[maxLab]==nCl ? lbID[maxLab] : -lbID[maxLab]); // winner label
dd2117a2 1364 return;
1365 }
1366 //
1367 hyp->SetFakeRatio(-1.);
5785a9d8 1368 hyp->SetLabel( -tpcLab );
dd2117a2 1369 hyp->SetITSLabel( kDummyLabel );
1370}
38997c3c 1371
1372//__________________________________________________________________
1373Bool_t AliITSUTrackerGlo::AddSeedBranch(AliITSUSeed* seed)
1374{
1375 // add new prolongation branch for the seed to temporary array. The seeds are sorted according to Compare f-n
6cd80116 1376 if (fNCandidatesAdded+fNBranchesAdded>=fLayerMaxCandidates ) {
1377 AliInfo(Form("Number of candidates at layer %d for current seed reached %d, increasing buffer",fCurrActLrID,fLayerMaxCandidates));
1378 fLayerMaxCandidates = 2*(fLayerMaxCandidates+1);
1379 AliITSUSeed** tmpArr = fLayerCandidates;
1380 fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
1381 memcpy(fLayerCandidates,tmpArr,(fNCandidatesAdded+fNBranchesAdded)*sizeof(AliITSUSeed*));
1382 delete tmpArr; // delete only array, not objects
38997c3c 1383 }
1384 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded]; // note: fNCandidatesAdded is incremented after adding all branches of current seed
1385 int slot=fNBranchesAdded++;
1386 for (int slotF=slot;slotF--;) { // slotF is always slot-1
1387 AliITSUSeed* si = branches[slotF];
1388 if (si->Compare(seed)<0) break; // found the last seed with better quality
1389 // otherwise, shift the worse seed to the next slot
1390 branches[slot] = si;
1391 slot = slotF; // slot should be slotF+1
1392 }
1393 // if needed, move worse seeds
1394 branches[slot] = seed;
1395 return kTRUE;
1396 //
1397}
1398
1399//__________________________________________________________________
1400void AliITSUTrackerGlo::ValidateAllowedBranches(Int_t acceptMax)
1401{
1402 // keep allowed number of branches for current seed and disable extras
1403 int nb = Min(fNBranchesAdded,acceptMax);
1404 // if (nb<fNBranchesAdded) printf("ValidateAllowedBranches: %d of %d (%d) on lr %d\n",nb,fNBranchesAdded,fNCandidatesAdded,ilr);
1405 // disable unused branches
1406 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded];
c03e4f8a 1407#ifdef _FILL_CONTROL_HISTOS_
1408 int bestID = -1;
1409 for (int ib=0;ib<fNBranchesAdded;ib++) {
1410 AliITSUSeed* sd = branches[ib];
1411 if (!sd->ContainsFake() && (bestID<0 || sd->Compare(branches[bestID])<0) ) bestID = ib;
1412 }
1413 if (bestID>=0) {
bdb92766 1414 TH2* hb = (TH2*)fCHistoArr->At(kHBestInBranch + (1+fTrackPhase)*kHistosPhase + fCurrActLrID);
c03e4f8a 1415 if (hb) hb->Fill(branches[bestID]->Pt(), bestID);
1416 }
1417#endif
bdb92766 1418 //
38997c3c 1419 for (int ib=nb;ib<fNBranchesAdded;ib++) {
bdb92766 1420 //
1421#ifdef _FILL_CONTROL_HISTOS_
1422 if (AliDebugLevelClass()>-2 && !branches[ib]->ContainsFake() /*&& branches[ib]->GetNLayersHit()*/
1423 && (bestID<0 || branches[ib]->Compare(branches[bestID])<0 ) ) {
38997c3c 1424 printf("Suppress good branch as %d of %d |",ib,fNBranchesAdded); branches[ib]->Print();
bdb92766 1425 // printf("Survivors : \n");
1426 // for (int j=0;j<nb;j++) branches[j]->Print();
38997c3c 1427 }
bdb92766 1428#endif
38997c3c 1429 MarkSeedFree(branches[ib]);
1430 }
1431 fNCandidatesAdded += nb; // update total candidates counter
1432 fNBranchesAdded = 0; // reset branches counter
1433 //
1434}
1435
1436//__________________________________________________________________
1437void AliITSUTrackerGlo::ValidateAllowedCandidates(Int_t ilr, Int_t acceptMax)
1438{
1439 // transfer allowed number of branches to hypothesis container
1440 //
1441 // sort candidates in increasing order of chi2
6cd80116 1442 static int lastSize = 0;
1443 static int *index = 0;
1444 static float *chi2 = 0;
1445 if (fLayerMaxCandidates>lastSize) {
1446 lastSize = fLayerMaxCandidates;
1447 delete[] index;
1448 delete[] chi2;
1449 index = new int[lastSize];
1450 chi2 = new float[lastSize];
1451 }
38997c3c 1452 for (int i=fNCandidatesAdded;i--;) chi2[i] = fLayerCandidates[i]->GetChi2GloNrm();
1453 Sort(fNCandidatesAdded,chi2,index,kFALSE);
1454 //
c03e4f8a 1455#ifdef _FILL_CONTROL_HISTOS_
1456 int bestID = -1;
1457 for (int ib=0;ib<fNCandidatesAdded;ib++) {
1458 AliITSUSeed* sd = fLayerCandidates[index[ib]];
1459 if (!sd->ContainsFake() && (bestID<0 || sd->Compare(fLayerCandidates[index[bestID]])<0) ) bestID = ib;
1460 }
1461 if (bestID>=0) {
bdb92766 1462 TH2* hb = (TH2*)fCHistoArr->At(kHBestInCand + (1+fTrackPhase)*kHistosPhase + fCurrActLrID);
c03e4f8a 1463 if (hb) hb->Fill(fLayerCandidates[index[bestID]]->Pt(), bestID);
1464 }
1465#endif
1466 //
38997c3c 1467 int nb = Min(fNCandidatesAdded,acceptMax);
38997c3c 1468 //
1469 for (int ib=0;ib<nb;ib++) AddProlongationHypothesis(fLayerCandidates[index[ib]],ilr);
1470 // disable unused candidates
1471 for (int ib=nb;ib<fNCandidatesAdded;ib++) {
bdb92766 1472 //
1473#ifdef _FILL_CONTROL_HISTOS_
1474 if (AliDebugLevelClass()>-2 && !fLayerCandidates[index[ib]]->ContainsFake() /*&& fLayerCandidates[index[ib]]->GetNLayersHit()*/
1475 && (bestID<0 || fLayerCandidates[index[ib]]->Compare(fLayerCandidates[index[bestID]])<0 ) ) {
38997c3c 1476 printf("Suppress good candidate as %d of %d |",index[ib],fNCandidatesAdded); fLayerCandidates[index[ib]]->Print();
1477 }
bdb92766 1478#endif
38997c3c 1479 MarkSeedFree(fLayerCandidates[index[ib]]);
1480 }
1481 fNCandidatesAdded = 0; // reset candidates counter
1482 //
1483}
1484
6cd80116 1485//__________________________________________________________________
1486void AliITSUTrackerGlo::SaveReducedHypothesesTree(AliITSUTrackHyp* dest)
1487{
1488 // remove those hypothesis seeds which dont lead to candidates at final layer
1489 //
1490 // 1st, flag the seeds to save
1491 int lr0 = 0;
1492 for (int isd=0;isd<fCurrHyp->GetNSeeds(lr0);isd++) {
1493 AliITSUSeed* seed = fCurrHyp->RemoveSeed(lr0,isd);
1494 if (!seed) continue;
1495 seed->FlagTree(AliITSUSeed::kSave);
1496 dest->AddSeed(seed,lr0);
1497 }
1498 for (int ilr=1;ilr<fNLrActive;ilr++) {
1499 int nsd = fCurrHyp->GetNSeeds(ilr);
1500 for (int isd=0;isd<nsd;isd++) {
1501 AliITSUSeed* seed = fCurrHyp->RemoveSeed(ilr,isd);
1502 if (!seed) continue; // already discarded or saved
1503 if (seed->IsSaved()) dest->AddSeed(seed,ilr);
1504 else MarkSeedFree(seed);
1505 }
1506 }
1507 //
1508 // AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
1509}
1510
53870004 1511//__________________________________________________________________
1512void AliITSUTrackerGlo::FlagSplitClusters()
1513{
1514 // set special bit on split clusters using MC info
6cd80116 1515 for (int ilr=fNLrActive;ilr--;) {
53870004 1516 int nsplit=0;
1517 AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
1518 for (int isn=lr->GetNSensors();isn--;) {
1519 AliITSURecoSens* sens = lr->GetSensor(isn);
1520 int nCl = sens->GetNClusters();
1521 if (!nCl) continue;
1522 int cl0 = sens->GetFirstClusterId();
1523 for (int icl=nCl;icl--;) {
1524 AliITSUClusterPix *cl = (AliITSUClusterPix*)lr->GetCluster(cl0+icl);
1525 for (int icl1=icl;icl1--;) {
1526 AliITSUClusterPix *cl1 = (AliITSUClusterPix*)lr->GetCluster(cl0+icl1);
1527 if (cl->HasCommonTrack(cl1)) {
1528 if (!cl->IsSplit()) nsplit++;
1529 if (!cl1->IsSplit()) nsplit++;
1530 cl->SetSplit();
1531 cl1->SetSplit();
1532 }
1533 }
1534 }
1535 }
1536 AliInfo(Form("%4d out of %4d clusters are split",nsplit,lr->GetNClusters()));
1537 }
1538 //
1539}
1540
1541//__________________________________________________________________
1542Bool_t AliITSUTrackerGlo::ContainsSplitCluster(const AliITSUSeed* seed, Int_t maxSize)
1543{
1544 // check if the seed contains split cluster with size < maxSize
1545 int lrID,clID;
1546 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1547 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1548 if (cl->IsSplit() && cl->GetNPix()<maxSize ) return kTRUE;
1549 }
1550 return seed->GetParent() ? ContainsSplitCluster((AliITSUSeed*)seed->GetParent(),maxSize) : kFALSE;
1551 //
1552}
1553
1554//__________________________________________________________________
1555void AliITSUTrackerGlo::PrintSeedClusters(const AliITSUSeed* seed, Option_t* option)
1556{
1557 // print seeds clusters
1558 int lrID,clID;
1559 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1560 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1561 cl->Print(option);
1562 }
1563 if (seed->GetParent()) PrintSeedClusters((AliITSUSeed*)seed->GetParent(), option);
1564 //
1565}
1566
c03e4f8a 1567//__________________________________________________________________
1568void AliITSUTrackerGlo::FlagSeedClusters(const AliITSUSeed* seed, Bool_t flg)
1569{
1570 // mark used clusters
1571 int lrID,clID;
1572 while (seed) {
bdb92766 1573 if ( (clID=seed->GetLrCluster(lrID))>=0 ) ((AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID))->ModClUsage(flg);
c03e4f8a 1574 seed = (AliITSUSeed*)seed->GetParent();
1575 }
1576 //
1577}
1578
1579
1580
53870004 1581#ifdef _FILL_CONTROL_HISTOS_
1582//__________________________________________________________________
1583void AliITSUTrackerGlo::BookControlHistos()
1584{
1585 // book special control histos
bdb92766 1586 if (fCHistoArr) return;
1587 const int kNResDef=7;
1588 const double kResDef[kNResDef]={0.05,0.05,0.3, 0.05,1,0.5,1.5};
1589 fCHistoArr = new TObjArray();
1590 fCHistoArr->SetOwner(kTRUE);
1591 const double ptMax=10;
1592 const double plMax=10;
1593 const double chiMax=100;
1594 const int nptbins=50;
1595 const int nresbins=400;
1596 const int nplbins=50;
1597 const int nchbins=200;
1598 const int maxBr = 15;
1599 const int maxCand = 200;
1600 TString ttl;
1601 for (int stp=0;stp<kNTrackingPhases;stp++) {
1602 for (int ilr=0;ilr<fNLrActive;ilr++) {
1603 int hoffs = (1+stp)*kHistosPhase + ilr;
1604 double mxdf = ilr>=kNResDef ? kResDef[kNResDef-1] : kResDef[ilr];
1605 ttl = Form("S%d_residY%d",stp,ilr);
1606 TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1607 fCHistoArr->AddAtAndExpand(hdy,hoffs + kHResY);
1608 hdy->SetDirectory(0);
1609 //
1610 ttl = Form("S%d_residYPull%d",stp,ilr);
1611 TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1612 fCHistoArr->AddAtAndExpand(hdyp,hoffs + kHResYP);
1613 hdyp->SetDirectory(0);
1614 //
1615 ttl = Form("S%d_residZ%d",stp,ilr);
1616 TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1617 fCHistoArr->AddAtAndExpand(hdz,hoffs + kHResZ);
1618 hdz->SetDirectory(0);
1619 //
1620 ttl = Form("S%d_residZPull%d",stp,ilr);
1621 TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1622 hdzp->SetDirectory(0);
1623 fCHistoArr->AddAtAndExpand(hdzp,hoffs + kHResZP);
1624 //
1625 ttl = Form("S%d_chi2Cl%d",stp,ilr);
1626 TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
1627 hchi->SetDirectory(0);
1628 fCHistoArr->AddAtAndExpand(hchi,hoffs + kHChi2Cl);
1629 //
1630 ttl = Form("S%d_chi2Nrm%d",stp,ilr);
1631 TH2F* hchiN = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
1632 hchiN->SetDirectory(0);
1633 fCHistoArr->AddAtAndExpand(hchiN,hoffs + kHChi2Nrm);
1634 //
1635 if (stp==0) { // these histos make sense only for clusters2tracks stage
1636 ttl = Form("S%d_bestInBranch%d",stp,ilr);
1637 TH2* hnbr = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxBr,-0.5,maxBr-0.5);
1638 hnbr->SetDirectory(0);
1639 fCHistoArr->AddAtAndExpand(hnbr,hoffs + kHBestInBranch);
c03e4f8a 1640 //
bdb92766 1641 ttl = Form("S%d_bestInCands%d",stp,ilr);
1642 TH2* hncn = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxCand,-0.5,maxCand-0.5);
1643 hncn->SetDirectory(0);
1644 fCHistoArr->AddAtAndExpand(hncn,hoffs + kHBestInCand);
c03e4f8a 1645 //
53870004 1646 }
1647 }
1648 }
bdb92766 1649 // custom histos
1650 ttl = Form("ClSharing");
1651 TH2* hclShare = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, 2*fNLrActive,0,2*fNLrActive);
1652 hclShare->SetDirectory(0);
1653 fCHistoArr->AddAtAndExpand(hclShare,kHClShare);
53870004 1654 //
1655}
1656#endif