Modifications needed to use PID framework based mass during tracking and
[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"
c61e50c3 38#include "AliITSUClusterPix.h"
cb50e082 39#include "AliITSUGeomTGeo.h"
bdb92766 40#include "AliCodeTimer.h"
9cdcba2c 41#include "AliRefArray.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)
e1ef49ad 56 ,fMatLUT(0)
0091e9f0 57 ,fCurrESDtrack(0)
0ddbf657 58 ,fCurrESDtrMClb(kDummyLabel)
32d38de2 59 ,fCurrMass(kPionMass)
0ddbf657 60 ,fCountProlongationTrials(0)
61 ,fCountITSin(0)
62 ,fCountITSout(0)
63 ,fCountITSrefit(0)
d99994b2 64 ,fNTracksESD(0)
716ccba7 65 ,fHypStore(100)
6cd80116 66 ,fLayerMaxCandidates(1000)
67 ,fLayerCandidates(0)
38997c3c 68 ,fNBranchesAdded(0)
69 ,fNCandidatesAdded(0)
716ccba7 70 ,fCurrHyp(0)
6cd80116 71 ,fWorkHyp(0)
c61e50c3 72 ,fSeedsPool("AliITSUSeed",0)
b8b59e05 73 ,fFreeSeedsID(0)
9cdcba2c 74 ,fESDIndex(0)
b8b59e05 75 ,fNFreeSeeds(0)
76 ,fLastSeedID(0)
6cd80116 77 ,fNLrActive(0)
42c3d4bd 78 ,fDefTrackConds(0)
79 ,fCurrTrackCond(0)
c03e4f8a 80 ,fCurrActLrID(-1)
81 ,fCurrLayer(0)
781c79ac 82 ,fTrackPhaseID(-1)
83 ,fCurrPassID(-1)
e1ef49ad 84 ,fUseMatLUT(kFALSE)
08419930 85#ifdef _ITSU_TUNING_MODE_
86 ,fCHistoArrCorr(0)
87 ,fCHistoArrFake(0)
53870004 88#endif
32d38de2 89{
90 // Default constructor
91 if (rec) Init(rec);
92}
93
94//_________________________________________________________________________
95AliITSUTrackerGlo::~AliITSUTrackerGlo()
96{
97 // Default destructor
98 //
9cdcba2c 99 delete[] fLayerCandidates;
100 if (fWorkHyp) fWorkHyp->SetTPCSeed(0); // this hypothesis does not own the seed
101 delete fWorkHyp;
e1ef49ad 102 delete fMatLUT;
32d38de2 103 //
08419930 104#ifdef _ITSU_TUNING_MODE_
105 if (fCHistoArrCorr || fCHistoArrFake) {
76390254 106 TFile* ctrOut = TFile::Open("itsuControlHistos.root","recreate");
107 ctrOut->cd();
108 AliInfo("Storing control histos");
53870004 109 // ctrOut->WriteObject(fCHistoArr,"controlH","kSingleKey");
08419930 110 if (fCHistoArrCorr) {
111 fCHistoArrCorr->Write();
112 delete fCHistoArrCorr;
113 }
114 if (fCHistoArrFake) {
115 fCHistoArrFake->Write();
116 delete fCHistoArrFake;
117 }
76390254 118 ctrOut->Close();
119 delete ctrOut;
08419930 120 fCHistoArrCorr = 0;
121 fCHistoArrFake = 0;
76390254 122 }
123#endif
124 //
32d38de2 125}
126
127//_________________________________________________________________________
128void AliITSUTrackerGlo::Init(AliITSUReconstructor* rec)
129{
130 // init with external reconstructor
c61e50c3 131 //
3d4dc3e2 132 fITS = rec->GetITSInterface();
6cd80116 133 fNLrActive = fITS->GetNLayersActive();
6cd80116 134 fWorkHyp = new AliITSUTrackHyp(fNLrActive);
e1ef49ad 135 // create material lookup table
136 const int kNTest = 1000;
137 const double kStepsPerCM=5;
138 fMatLUT = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(),Nint(kStepsPerCM*(fITS->GetRMax()-fITS->GetRMin())));
139 double zmn = 1e6;
140 for (int ilr=fITS->GetNLayers();ilr--;) {
141 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
142 if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin());
143 if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax());
144 }
145 fMatLUT->FillData(kNTest,-zmn,zmn);
c61e50c3 146 //
6cd80116 147 if (fLayerMaxCandidates<1) fLayerMaxCandidates = 1000;
148 fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
c61e50c3 149 fSeedsPool.ExpandCreateFast(1000); // RS TOCHECK
b8b59e05 150 fFreeSeedsID.Set(1000);
9cdcba2c 151 fESDIndex.Set(1000);
c61e50c3 152 //
32d38de2 153}
154
155//_________________________________________________________________________
42c3d4bd 156void AliITSUTrackerGlo::CreateDefaultTrackCond()
157{
158 // creates default tracking conditions to be used when recoparam does not provide them
6cd80116 159
42c3d4bd 160 AliITSUTrackCond* cond = new AliITSUTrackCond();
161 //
6cd80116 162 cond->SetNLayers(fNLrActive);
163 cond->AddNewCondition(fNLrActive);
70f61d86 164 cond->AddGroupPattern( 0xffff, fNLrActive ); // require all layers hit
ee54014a 165 cond->Init();
42c3d4bd 166 //
167 fDefTrackConds.AddLast(cond);
168 //
169 AliInfo("Created conditions: ");
170 for (int i=0;i<fDefTrackConds.GetEntriesFast();i++) fDefTrackConds[i]->Print();
171 //
172}
173
174
175//_________________________________________________________________________
32d38de2 176Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv)
177{
178 //
bdb92766 179 AliCodeTimerAuto("",0);
76390254 180 SetTrackingPhase(kClus2Tracks);
181 //
08419930 182#ifdef _ITSU_TUNING_MODE_
781c79ac 183 if (!fCHistoArrCorr) fCHistoArrCorr = BookControlHistos("Corr");
184 if (!fCHistoArrFake) fCHistoArrFake = BookControlHistos("Fake");
76390254 185#endif
c03e4f8a 186 static int evID = 0;
9cdcba2c 187 static TArrayF esdTrPt(fESDIndex.GetSize());
bdb92766 188 //
42c3d4bd 189 TObjArray *trackConds = 0;
32d38de2 190 //
0ddbf657 191 fCountProlongationTrials = 0;
192 fCountITSin = 0;
193 fCountITSout = 0;
194 fCountITSrefit = 0;
195 //
b8b59e05 196 ResetSeedsPool();
d99994b2 197 fNTracksESD = esdEv->GetNumberOfTracks();
198 AliInfo(Form("Will try to find prolongations for %d tracks",fNTracksESD));
42c3d4bd 199 int nTrackCond = AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();
e1ef49ad 200 fUseMatLUT = AliITSUReconstructor::GetRecoParam()->GetUseMatLUT(fTrackPhaseID);
42c3d4bd 201 if (nTrackCond<1) {
202 if (!fDefTrackConds.GetEntriesFast()) {
203 AliInfo("No tracking conditions found in recoparams, creating default one requesting all layers hit");
204 CreateDefaultTrackCond();
205 }
206 trackConds = &fDefTrackConds;
207 nTrackCond = trackConds->GetEntriesFast();
208 }
209 else trackConds = AliITSUReconstructor::GetRecoParam()->GetTrackingConditions();
210 //
15b02d69 211 static Bool_t first = kTRUE;
212 if (first) {
213 AliITSUReconstructor::GetRecoParam()->Print();
214 first = kFALSE;
215 }
216 fHypStore.Delete();
d99994b2 217 if (fHypStore.GetSize()<fNTracksESD) fHypStore.Expand(fNTracksESD+100);
15b02d69 218 //
219 fITS->ProcessClusters();
220 //
08419930 221#ifdef _ITSU_TUNING_MODE_
53870004 222 FlagSplitClusters(); // tmp RS
bdb92766 223#endif
224 //
225 // the tracks will be reconstructed in decreasing pt order, sort them
d99994b2 226 if (fESDIndex.GetSize()<fNTracksESD) {
227 fESDIndex.Set(fNTracksESD+200);
228 esdTrPt.Set(fNTracksESD+200);
bdb92766 229 }
9cdcba2c 230 int* esdTrackIndex = fESDIndex.GetArray();
bdb92766 231 float* trPt = esdTrPt.GetArray();
891efa00 232 for (int itr=fNTracksESD;itr--;) {
233 AliESDtrack* esdTr = esdEv->GetTrack(itr);
0e40ce66 234 esdTr->SetStatus(AliESDtrack::kITSupg); // Flag all tracks as participating in the ITSup reco
891efa00 235 trPt[itr] = esdTr->Pt();
236 }
d99994b2 237 Sort(fNTracksESD,trPt,esdTrackIndex,kTRUE);
53870004 238 //
42c3d4bd 239 for (int icnd=0;icnd<nTrackCond;icnd++) {
781c79ac 240 fCurrPassID = icnd;
42c3d4bd 241 fCurrTrackCond = (AliITSUTrackCond*)trackConds->UncheckedAt(icnd);
ee54014a 242 if (!fCurrTrackCond->IsInitDone()) fCurrTrackCond->Init();
42c3d4bd 243 // select ESD tracks to propagate
d99994b2 244 for (int itr=0;itr<fNTracksESD;itr++) {
bdb92766 245 int trID = esdTrackIndex[itr];
246 fCurrESDtrack = esdEv->GetTrack(trID);
0ddbf657 247 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
248 //
0e40ce66 249
250 if (!NeedToProlong(fCurrESDtrack,trID)) continue; // are we interested in this track?
c03e4f8a 251 /*
9cdcba2c 252 // if specific tracks should be checked, create a global array int wtc[] = {evITS*10000+trID, ...};
c03e4f8a 253 Bool_t dbg = kFALSE;
254 int nwtc = sizeof(wtc)/sizeof(int);
255
256 for (int k=0;k<nwtc;k++) {
bdb92766 257 if (wtc[k]==evID*10000+trID) {
c03e4f8a 258 dbg = kTRUE;
c03e4f8a 259 printf("At esdTr: %d MC: %d\n",wtc[k],fCurrESDtrMClb);
c03e4f8a 260 break;
261 }
262 }
263 AliLog::SetClassDebugLevel("AliITSUTrackerGlo",dbg ? 10:0);
264 */
e1ef49ad 265#ifdef _ITSU_DEBUG_
1d26da6d 266 AliDebug(1,Form("Processing track %d(esd%d) | M=%.3f Pt=%.3f | MCLabel: %d",itr,trID,fCurrESDtrack->GetMassForTracking(kTRUE),fCurrESDtrack->Pt(),fCurrESDtrMClb));//RS
e1ef49ad 267#endif
bdb92766 268 FindTrack(fCurrESDtrack, trID);
42c3d4bd 269 }
270 //
e1ef49ad 271#ifdef _ITSU_DEBUG_
38997c3c 272 if (AliDebugLevelClass()>+2) {
15b02d69 273 AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
274 fHypStore.Print();
275 }
e1ef49ad 276#endif
42c3d4bd 277 FinalizeHypotheses();
891efa00 278#ifdef _ITSU_TUNING_MODE_
d99994b2 279 CheckClusterUsage(); //!!RS
891efa00 280#endif
c03e4f8a 281 //
282 //AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0); // in case wtc array was defined, uncomment this
42c3d4bd 283 }
716ccba7 284 //
d99994b2 285 AliInfo(Form("%d ITSin for %d tried TPC seeds out of %d ESD tracks\n",fCountITSin,fCountProlongationTrials,fNTracksESD));
c03e4f8a 286 evID++;
32d38de2 287 return 0;
288}
289
290//_________________________________________________________________________
c51c3124 291Int_t AliITSUTrackerGlo::PropagateBack(AliESDEvent *esdEv)
32d38de2 292{
293 //
c51c3124 294 // Do outward fits in ITS
bdb92766 295 AliCodeTimerAuto("",0);
c51c3124 296 //
e7d83d38 297 SetTrackingPhase(kPropBack);
d99994b2 298 fNTracksESD = esdEv->GetNumberOfTracks();
e1ef49ad 299 fUseMatLUT = AliITSUReconstructor::GetRecoParam()->GetUseMatLUT(fTrackPhaseID);
300#ifdef _ITSU_DEBUG_
d99994b2 301 AliDebug(1,Form("Will propagate back %d tracks",fNTracksESD));
e1ef49ad 302#endif
c51c3124 303 //
304 double bz0 = GetBz();
305 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
9cdcba2c 306 AliITSUTrackHyp dummyTr;
c51c3124 307 const double kWatchStep=10.; // for larger steps watch arc vs segment difference
308 Double_t times[AliPID::kSPECIES];
309 //
d99994b2 310 for (int itr=0;itr<fNTracksESD;itr++) {
0ddbf657 311 fCurrESDtrack = esdEv->GetTrack(itr);
312 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
c51c3124 313 // Start time integral and add distance from current position to vertex
0ddbf657 314 if (fCurrESDtrack->IsOn(AliESDtrack::kITSout)) continue; // already done
c51c3124 315 //
0ddbf657 316 fCurrESDtrack->GetXYZ(xyzTrk);
c51c3124 317 Double_t dst = 0.; // set initial track lenght, tof
318 {
319 double dxs = xyzTrk[0] - xyzVtx[0];
320 double dys = xyzTrk[1] - xyzVtx[1];
7ce7f445 321 // double dzs = xyzTrk[2] - xyzVtx[2];
c51c3124 322 // RS: for large segment steps d use approximation of cicrular arc s by
323 // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
324 // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
325 // Hence s^2/d^2 = (1+1/6 p^2)^2
326 dst = dxs*dxs + dys*dys;
327 if (dst > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
0ddbf657 328 double crv = Abs(fCurrESDtrack->GetC(bz0));
c51c3124 329 double fcarc = 1.+crv*crv*dst/6.;
330 dst *= fcarc*fcarc;
331 }
7ce7f445 332 // RS: temporary hack since we don't have SPD vertex:
333 // dst += dzs*dzs;
334 dst *= 1+fCurrESDtrack->GetTgl()*fCurrESDtrack->GetTgl();
c51c3124 335 dst = Sqrt(dst);
336 }
337 //
0ddbf657 338 fCurrESDtrack->SetStatus(AliESDtrack::kTIME);
c51c3124 339 //
0ddbf657 340 if (!fCurrESDtrack->IsOn(AliESDtrack::kITSin)) { // case of tracks w/o ITS prolongation: just set the time integration
341 dummyTr.AliExternalTrackParam::operator=(*fCurrESDtrack);
c51c3124 342 dummyTr.StartTimeIntegral();
e2d2481c 343 dummyTr.AddTimeStep(dst);
c51c3124 344 dummyTr.GetIntegratedTimes(times);
0ddbf657 345 fCurrESDtrack->SetIntegratedTimes(times);
346 fCurrESDtrack->SetIntegratedLength(dummyTr.GetIntegratedLength());
c51c3124 347 continue;
348 }
349 //
9cdcba2c 350 fCurrHyp = GetTrackHyp(itr);
70cb7fe4 351 fCurrMass = fCurrHyp->GetMass();
9cdcba2c 352 fCurrHyp->StartTimeIntegral();
353 fCurrHyp->AddTimeStep(dst);
9cdcba2c 354 fCurrHyp->ResetCovariance(10000);
5e3587be 355 int nclFit = 0;
356 double chi2 = RefitTrack(fCurrHyp,fITS->GetRMax(),nclFit);
9cdcba2c 357 if (chi2>0) { // propagate to exit from the ITS/TPC screen
5e3587be 358 int ndf = nclFit*2-5;
70cb7fe4 359 if (ndf>0) chi2 /= ndf;
9cdcba2c 360 fCurrHyp->SetChi2(chi2);
361 UpdateESDTrack(fCurrHyp,AliESDtrack::kITSout);
68a0f687 362 }
e7d83d38 363 else {
e1ef49ad 364#ifdef _ITSU_DEBUG_
0ddbf657 365 AliDebug(2,Form("Refit Failed for track %d | ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
e1ef49ad 366#endif
9cdcba2c 367 //fCurrHyp->AliExternalTrackParam::Print();
368 //fCurrHyp->GetWinner()->Print();
e7d83d38 369 }
c51c3124 370 //
371 }
32d38de2 372 //
0ddbf657 373 AliInfo(Form("%d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
d99994b2 374 fCountITSout,fCountITSin,fCountProlongationTrials,fNTracksESD));
0ddbf657 375 //
32d38de2 376 return 0;
377}
378
379//_________________________________________________________________________
e7d83d38 380Int_t AliITSUTrackerGlo::RefitInward(AliESDEvent *esdEv)
32d38de2 381{
382 //
e7d83d38 383 // refit the tracks inward, using current cov.matrix
bdb92766 384 AliCodeTimerAuto("",0);
e7d83d38 385 //
386 SetTrackingPhase(kRefitInw);
e1ef49ad 387 fUseMatLUT = AliITSUReconstructor::GetRecoParam()->GetUseMatLUT(fTrackPhaseID);
d99994b2 388 fNTracksESD = esdEv->GetNumberOfTracks();
b515ad7e 389 // AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
390
d99994b2 391 AliDebug(1,Form("Will refit inward %d tracks",fNTracksESD));
c7e63626 392 // Bool_t uselogMS = AliExternalTrackParam::GetUseLogTermMS();
393 // AliExternalTrackParam::SetUseLogTermMS(kTRUE);
e7d83d38 394 //
d99994b2 395 for (int itr=0;itr<fNTracksESD;itr++) {
0ddbf657 396 fCurrESDtrack = esdEv->GetTrack(itr);
397 fCurrESDtrMClb = fCurrESDtrack->GetLabel();
e7d83d38 398 // Start time integral and add distance from current position to vertex
0ddbf657 399 UInt_t trStat = fCurrESDtrack->GetStatus();
e7d83d38 400 if ( !(trStat & AliESDtrack::kITSout) ) continue;
401 if ( trStat & AliESDtrack::kITSrefit ) continue; // already done
402 if ( (trStat & AliESDtrack::kTPCout) && !(trStat & AliESDtrack::kTPCrefit) ) continue;
403 //
7ce7f445 404 fCurrHyp = GetTrackHyp(itr);
70cb7fe4 405 fCurrMass = fCurrHyp->GetMass();
7ce7f445 406 *fCurrHyp = *fCurrESDtrack; // fetch current ESDtrack kinematics
70cb7fe4 407 //
5e3587be 408 int nclFit = 0;
eb18cdc6 409 // fCurrHyp->ResetCovariance(1000); // in case we want to do SA refit
5e3587be 410 double chi2 = RefitTrack(fCurrHyp,fITS->GetRMin(),nclFit);
70cb7fe4 411 if (chi2>0) { // propagate up to inside radius of the beam pipe
5e3587be 412 int ndf = nclFit*2-5;
413 if (ndf>0) chi2 /= ndf;
9cdcba2c 414 fCurrHyp->SetChi2(chi2);
415 UpdateESDTrack(fCurrHyp,AliESDtrack::kITSrefit);
e7d83d38 416 }
417 else {
e1ef49ad 418#ifdef _ITSU_DEBUG_
0ddbf657 419 AliDebug(2,Form("Refit Failed for track %d |ESDtrack#%d (MClb:%d)",itr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
e1ef49ad 420#endif
e7d83d38 421 }
422 }
c7e63626 423 // AliExternalTrackParam::SetUseLogTermMS(uselogMS);
32d38de2 424 //
0ddbf657 425 AliInfo(Form("%d ITSrefit in %d ITSout in %d ITSin tracks for %d tried TPC seeds out of %d ESD tracks\n",
d99994b2 426 fCountITSrefit,fCountITSout,fCountITSin,fCountProlongationTrials,fNTracksESD));
0ddbf657 427 //
b515ad7e 428 // AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0);
32d38de2 429 return 0;
430}
431
432//_________________________________________________________________________
433Int_t AliITSUTrackerGlo::LoadClusters(TTree * treeRP)
434{
435 // read from tree (if pointer provided) or directly from the ITS reco interface
436 //
437 return fReconstructor->LoadClusters(treeRP);
438}
439
440//_________________________________________________________________________
441void AliITSUTrackerGlo::UnloadClusters()
442{
443 //
444 // To be implemented
445 //
446
447 Info("UnloadClusters","To be implemented");
448}
449//_________________________________________________________________________
450AliCluster * AliITSUTrackerGlo::GetCluster(Int_t /*index*/) const
451{
452 //
453 // To be implemented
454 //
455
456 Info("GetCluster","To be implemented");
457 return 0x0;
458}
459
460//_________________________________________________________________________
0e40ce66 461Bool_t AliITSUTrackerGlo::NeedToProlong(AliESDtrack* esdTr, Int_t esdID)
32d38de2 462{
463 // do we need to match this track to ITS?
464 //
465 static double bz = GetBz();
0e40ce66 466 //
467 // is it already prolonged
468 if (esdTr->IsOn(AliESDtrack::kITSin)) return kFALSE;
469 AliITSUTrackHyp* hyp = GetTrackHyp(esdID); // in case there is unfinished hypothesis from prev.pass, clean it
470 if (hyp) {
471 CleanHypothesis(hyp);
472 if (hyp->GetSkip()) return kFALSE; // need to skip this
473 return kTRUE;
474 }
475 //
32d38de2 476 if (!esdTr->IsOn(AliESDtrack::kTPCin) ||
477 esdTr->IsOn(AliESDtrack::kTPCout) ||
478 esdTr->IsOn(AliESDtrack::kITSin) ||
479 esdTr->GetKinkIndex(0)>0) return kFALSE;
480 //
c61e50c3 481 if (esdTr->Pt()<AliITSUReconstructor::GetRecoParam()->GetMinPtForProlongation()) return kFALSE;
32d38de2 482 //
483 float dtz[2];
484 esdTr->GetDZ(GetX(),GetY(),GetZ(),bz,dtz);
485 // if track is not V0 candidata but has large offset wrt IP, reject it. RS TOCHECK
c61e50c3 486 if ( !(esdTr->GetV0Index(0)>0 && dtz[0]>AliITSUReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation())
487 && (Abs(dtz[0])>AliITSUReconstructor::GetRecoParam()->GetMaxDForProlongation() ||
488 Abs(dtz[1])>AliITSUReconstructor::GetRecoParam()->GetMaxDZForProlongation())) return kFALSE;
32d38de2 489 //
cb50e082 490 // if (esdTr->Pt()<3) return kFALSE;//RS
32d38de2 491 return kTRUE;
492}
493
494//_________________________________________________________________________
716ccba7 495void AliITSUTrackerGlo::FindTrack(AliESDtrack* esdTr, Int_t esdID)
32d38de2 496{
497 // find prolongaion candidates finding for single seed
498 //
3e4e3c23 499 AliITSUSeed seedUC; // copy of the seed from the upper layer
500 AliITSUSeed seedT; // transient seed between the seedUC and new prolongation hypothesis
501 //
d99994b2 502 // AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
6cd80116 503 AliITSUTrackHyp* hypTr = InitHypothesis(esdTr,esdID); // initialize prolongations hypotheses tree
9cdcba2c 504 if (!hypTr || hypTr->GetSkip()) return;
505 //
d99994b2 506 double bz = GetBz();
6cd80116 507 fCurrHyp = fWorkHyp;
508 fCurrHyp->InitFrom(hypTr);
32d38de2 509 //
173b3073 510 AliITSURecoSens *hitSens[AliITSURecoSens::kNNeighbors+1];
f8832015 511 //
d99994b2 512 int ilaUp = fNLrActive; // previous active layer really checked (some may be excluded!)
32d38de2 513 //
d99994b2 514 for (int ila=fNLrActive;ila--;ilaUp=fCurrActLrID) {
515 if (fCurrTrackCond->IsLayerExcluded(ila)) continue;
c03e4f8a 516 fCurrActLrID = ila;
517 fCurrLayer = fITS->GetLayerActive(ila);
518 Bool_t noClSharing = fCurrTrackCond->GetClSharing(ila)==0;
3e4e3c23 519 //
520 // for the outermost layer the seed is created from the ESD track
6cd80116 521 int nSeedsUp = (ilaUp==fNLrActive) ? 1 : fCurrHyp->GetNSeeds(ilaUp);
38997c3c 522 int maxNBranches = fCurrTrackCond->GetMaxBranches(ila);
523 int maxNCandidates = fCurrTrackCond->GetMaxCandidates(ila);
3e4e3c23 524 //
32d38de2 525 for (int isd=0;isd<nSeedsUp;isd++) {
3e4e3c23 526 AliITSUSeed* seedU;
6cd80116 527 if (ilaUp==fNLrActive) {
3e4e3c23 528 seedU = 0;
70cb7fe4 529 seedUC.InitFromSeed(fCurrHyp->GetTPCSeed()); // init with TPC seed from ref.R
3e4e3c23 530 }
531 else {
532 seedU = fCurrHyp->GetSeed(ilaUp,isd); // seed on prev.active layer to prolong
c03e4f8a 533 if (seedU->IsKilled()) continue;
3e4e3c23 534 seedUC = *seedU; // its copy will be prolonged
b8b59e05 535 seedUC.SetParent(seedU);
3e4e3c23 536 }
943e1898 537 seedUC.ResetFMatrix(); // reset the matrix for propagation to next layer
32d38de2 538 // go till next active layer
e1ef49ad 539#ifdef _ITSU_DEBUG_
0ddbf657 540 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()));
e1ef49ad 541#endif
70cb7fe4 542 if (!TransportToLayer(&seedUC, fITS->GetLrIDActive(ilaUp), fITS->GetLrIDActive(ila)) ) { // external seed already prolonged
32d38de2 543 //
e1ef49ad 544#ifdef _ITSU_DEBUG_
0ddbf657 545 AliDebug(2,Form("Transport failed | esdID=%d (MClb:%d)",esdID,fCurrESDtrMClb));
e1ef49ad 546#endif
32d38de2 547 // Check if the seed satisfies to track definition
b8b59e05 548 if (NeedToKill(&seedUC,kTransportFailed) && seedU) KillSeed(seedU,kTRUE);
32d38de2 549 continue; // RS TODO: decide what to do with tracks stopped on higher layers w/o killing
550 }
f8832015 551 if (!GetRoadWidth(&seedUC, ila)) { // failed to find road width on the layer
b8b59e05 552 if (NeedToKill(&seedUC,kRWCheckFailed) && seedU) KillSeed(seedU,kTRUE);
32d38de2 553 continue;
554 }
5785a9d8 555 /*
556 //RS toremove
557 int mcquest = -1;
c03e4f8a 558 if (!seedUC.ContainsFake() && AliDebugLevelClass()>2) {
5785a9d8 559 mcquest = fCurrESDtrMClb;
560 seedUC.Print("etp");
561 printf("FSParams: "); for (int ip=0;ip<kNTrImpData;ip++) printf("%+e ",fTrImpData[ip]); printf("\n");
562 }
563 //
c03e4f8a 564 int nsens = fCurrLayer->FindSensors(&fTrImpData[kTrPhi0], hitSens, mcquest); // find detectors which may be hit by the track
5785a9d8 565 */
c03e4f8a 566 int nsens = fCurrLayer->FindSensors(&fTrImpData[kTrPhi0], hitSens); // find detectors which may be hit by the track
e1ef49ad 567#ifdef _ITSU_DEBUG_
0ddbf657 568 AliDebug(2,Form("Will check %d sensors on lr:%d | esdID=%d (MClb:%d)",nsens,ila,esdID,fCurrESDtrMClb));
e1ef49ad 569#endif
32d38de2 570 //
38997c3c 571 seedUC.SetLr(ila);
572 //
32d38de2 573 for (int isn=nsens;isn--;) {
c61e50c3 574 AliITSURecoSens* sens = hitSens[isn];
38997c3c 575 int ncl = sens->GetNClusters();
576 if (!ncl) continue;
577 seedT = seedUC;
f8832015 578 //
c03e4f8a 579 // We need to propagate the seed to sensor on fCurrLayer staying in the frame of the sensor from prev.layer,
943e1898 580 // since the transport matrix should be defined in this frame.
581 double xs; // X in the TF of current seed, corresponding to intersection with sensor plane
53870004 582 if (!seedT.GetTrackingXAtXAlpha(sens->GetXTF(),sens->GetPhiTF(),bz, xs)) {
e1ef49ad 583#ifdef _ITSU_DEBUG_
53870004 584 if (AliDebugLevelClass()>2) {
585 printf("Failed on GetTrackingXAtXAlpha: X=%.4f alp=%.3f\n",sens->GetXTF(),sens->GetPhiTF());
586 seedT.Print("etp");
587 }
e1ef49ad 588#endif
53870004 589 continue;
590 }
591 if (xs<seedT.GetX()) {
592 if (!PropagateSeed(&seedT,xs,fCurrMass)) continue;
593 }
594 else { // some low precision tracks may hit the sensor plane outside of the layer radius
e1ef49ad 595#ifdef _ITSU_DEBUG_
53870004 596 if (AliDebugLevelClass()>2) {
597 if (!seedT.ContainsFake()) {
598 printf("WRONG PROP on L%d, sens%d of %d: %.4f > %.4f\n",ila,isn,nsens,xs,seedT.GetX());
599 sens->Print();
600 seedT.Print("etp");
601 }
602 }
e1ef49ad 603#endif
53870004 604 if (!seedT.PropagateParamOnlyTo(xs,bz)) continue;
605 }
8b16dbae 606 // if (!seedT.PropagateToX(xs,bz)) continue;
716ccba7 607 // if (!seedT.Rotate(sens->GetPhiTF())) continue;
608 if (!seedT.RotateToAlpha(sens->GetPhiTF())) continue;
943e1898 609 //
f8832015 610 int clID0 = sens->GetFirstClusterId();
cb50e082 611 for (int icl=ncl;icl--;) { // don't change the order, clusters are sorted
c03e4f8a 612 //AliLog::SetClassDebugLevel("AliITSUTrackerGlo",10);
613 int clID = clID0 + icl;
614 if (noClSharing && fCurrLayer->GetCluster(clID)->IsClusterUsed()) continue;
f8832015 615 int res = CheckCluster(&seedT,ila,clID0+icl);
c03e4f8a 616 //AliLog::SetClassDebugLevel("AliITSUTrackerGlo", 0);
f8832015 617 //
618 if (res==kStopSearchOnSensor) break; // stop looking on this sensor
619 if (res==kClusterNotMatching) continue; // cluster does not match
620 // cluster is matching and it was added to the hypotheses tree
c61e50c3 621 }
32d38de2 622 }
38997c3c 623 // cluster search is done. Do we need to have a version of this seed skipping current layer
624 if (!NeedToKill(&seedUC,kMissingCluster)) {
b2d3bc3a 625 AliITSUSeed* seedSkp = NewSeedFromPool(&seedUC);
097e4cdf 626 if (nsens) { // to do: make penalty to account for probability to miss the cluster for good reason
627 double penalty = -fCurrTrackCond->GetMissPenalty(ila);
628 seedSkp->SetChi2Cl(penalty);
629 }
08419930 630 if (seedSkp->GetChi2GloNrm()>fCurrTrackCond->GetMaxChi2GloNrm(ila)) {
631 MarkSeedFree(seedSkp);
632 }
633 else AddSeedBranch(seedSkp);
716ccba7 634 }
38997c3c 635 // transfer the new branches of the seed to the hypothesis container
636 if (fNBranchesAdded) ValidateAllowedBranches(maxNBranches);
637 //
32d38de2 638 }
38997c3c 639 if (fNCandidatesAdded) ValidateAllowedCandidates(ila,maxNCandidates);
640 // ((TObjArray*)fCurrHyp->GetLayerSeeds(ila))->Sort();
e1ef49ad 641#ifdef _ITSU_DEBUG_
38997c3c 642 if (AliDebugLevelClass()>2) { //RS
15b02d69 643 printf(">>> All hypotheses on lr %d: \n",ila);
644 for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {printf(" #%3d ",ih); fCurrHyp->GetSeed(ila,ih)->Print();}
645 }
e1ef49ad 646#endif
15b02d69 647 //
648 /*
716ccba7 649 for (int ih=0;ih<fCurrHyp->GetNSeeds(ila);ih++) {
716ccba7 650 if (ila!=0) continue;
651 double vecL[5] = {0};
652 double matL[15] = {0};
653 AliITSUSeed* sp = fCurrHyp->GetSeed(ila,ih);
654 while(sp->GetParent()) {
655 sp->Smooth(vecL,matL);
6cd80116 656 if (sp->GetLayerID()>=fNLrActive-1) break;
716ccba7 657 sp = (AliITSUSeed*)sp->GetParent();
658 }
659 */
32d38de2 660 }
661 //
6cd80116 662 SaveReducedHypothesesTree(hypTr);
e1ef49ad 663#ifdef _ITSU_DEBUG_
6cd80116 664 if (AliDebugLevelClass()>1) {
665 printf("\nSaved hypotheses for esdTrack %d (MCLab:%d)\n",esdID,fCurrESDtrMClb);
666 hypTr->Print("l");
667 }
e1ef49ad 668#endif
6cd80116 669 fCurrHyp = 0;
d99994b2 670 // AliLog::SetClassDebugLevel("AliITSUTrackerGlo",0);
716ccba7 671 //
32d38de2 672}
673
674//_________________________________________________________________________
6cd80116 675AliITSUTrackHyp* AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esdID)
32d38de2 676{
c61e50c3 677 // init prolongaion candidates finding for single seed
6cd80116 678 AliITSUTrackHyp* hyp = GetTrackHyp(esdID);
679 if (hyp) return hyp;
716ccba7 680 //
0ddbf657 681 fCountProlongationTrials++;
682 //
1d26da6d 683 fCurrMass = esdTr->GetMassForTracking();
c61e50c3 684 //
6cd80116 685 hyp = new AliITSUTrackHyp(fNLrActive);
686 hyp->SetESDTrack(esdTr);
687 hyp->SetUniqueID(esdID);
688 hyp->SetMass(fCurrMass);
9cdcba2c 689 hyp->SetTPCSeed( new AliExternalTrackParam(*esdTr) );
6cd80116 690 SetTrackHyp(hyp,esdID);
e1ef49ad 691 Bool_t lut = fUseMatLUT; // for propagation from TPC use TGeo material budget always
692 fUseMatLUT = kFALSE;
70cb7fe4 693 if (!TransportToLayer(hyp->GetTPCSeed(),fITS->GetNLayers(), fITS->GetLrIDActive(fNLrActive-1), fITS->GetRITSTPCRef())) hyp->SetSkip(); // propagate to outer R of ITS
e1ef49ad 694 fUseMatLUT = lut;
3e4e3c23 695 //
6cd80116 696 return hyp;
32d38de2 697 // TO DO
698}
699
700//_________________________________________________________________________
70cb7fe4 701Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
32d38de2 702{
70cb7fe4 703 // transport seed from layerFrom to the entrance (along the direction of the track) of layerTo or to rLim (if>0), wathever is closer
32d38de2 704 //
8b16dbae 705 //
706 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
707 //
32d38de2 708 int dir = lTo > lFrom ? 1:-1;
709 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
710 Bool_t checkFirst = kTRUE;
70cb7fe4 711 Bool_t limReached = kFALSE;
32d38de2 712 while(lFrom!=lTo) {
32d38de2 713 if (lrFr) {
70cb7fe4 714 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
e2d2481c 715 checkFirst = kFALSE;
32d38de2 716 }
717 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
718 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
719 //
720 // go the entrance of the layer, assuming no materials in between
8b16dbae 721 double xToGo = lrTo->GetR(-dir);
70cb7fe4 722 if (rLim>0) {
723 if (dir>0) {
724 if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
725 }
726 else {
727 if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
728 }
729 }
0ddbf657 730 // double xts = xToGo;
b2d3bc3a 731 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
0ddbf657 732 //printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
733 //seed->Print("etp");
70cb7fe4 734 //
b2d3bc3a 735 return kFALSE;
736 }
e1ef49ad 737#ifdef _ITSU_DEBUG_
b515ad7e 738 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
e1ef49ad 739#endif
b2d3bc3a 740 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
0ddbf657 741 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
742 //seed->Print("etp");
b2d3bc3a 743 return kFALSE;
744 }
32d38de2 745 lrFr = lrTo;
70cb7fe4 746 if (limReached) break;
32d38de2 747 }
748 return kTRUE;
749 //
750}
751
752//_________________________________________________________________________
70cb7fe4 753Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
3e4e3c23 754{
70cb7fe4 755 // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer
3e4e3c23 756 //
8b16dbae 757 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
758 //
3e4e3c23 759 int dir = lTo > lFrom ? 1:-1;
760 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
761 Bool_t checkFirst = kTRUE;
70cb7fe4 762 Bool_t limReached = kFALSE;
3e4e3c23 763 while(lFrom!=lTo) {
3e4e3c23 764 if (lrFr) {
e2d2481c 765 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
766 checkFirst = kFALSE;
3e4e3c23 767 }
768 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
769 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
770 //
0ddbf657 771 // go the entrance of the layer, assuming no materials in between
772 double xToGo = lrTo->GetR(-dir);
70cb7fe4 773 if (rLim>0) {
774 if (dir>0) {
775 if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
776 }
777 else {
778 if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
779 }
780 }
0ddbf657 781 // double xts = xToGo;
782 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
783 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
784 // seed->Print("etp");
785 return kFALSE;
786 }
e1ef49ad 787#ifdef _ITSU_DEBUG_
b515ad7e 788 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
e1ef49ad 789#endif
b515ad7e 790 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
791 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
792 //seed->Print("etp");
793 return kFALSE;
794 }
795 lrFr = lrTo;
70cb7fe4 796 if (limReached) break;
b515ad7e 797 }
798 return kTRUE;
799 //
800}
801
802//_________________________________________________________________________
803Bool_t AliITSUTrackerGlo::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
804{
805 // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
806 //
807 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
808 //
809 int dir = lTo > lFrom ? 1:-1;
810 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
811 Bool_t checkFirst = kTRUE;
812 while(lFrom!=lTo) {
813 if (lrFr) {
814 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
815 checkFirst = kFALSE;
816 }
817 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
818 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
819 //
820 // go the entrance of the layer, assuming no materials in between
821 double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
822 //
823 // double xts = xToGo;
824 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
825 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
826 // seed->Print("etp");
827 return kFALSE;
828 }
829 if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
830 //
e1ef49ad 831#ifdef _ITSU_DEBUG_
b515ad7e 832 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
e1ef49ad 833#endif
0ddbf657 834 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
835 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
836 //seed->Print("etp");
837 return kFALSE;
838 }
3e4e3c23 839 lrFr = lrTo;
840 }
841 return kTRUE;
842 //
843}
844
845//_________________________________________________________________________
e2d2481c 846Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
847{
848 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
849 // If check is requested, do this only provided the track has not exited the layer already
850 double xToGo = lr->GetR(dir);
851 if (check) { // do we need to track till the surface of the current layer ?
852 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
e1ef49ad 853#ifdef _ITSU_DEBUG_
15b02d69 854 AliDebug(3,Form(" dir:%d Cur: %e Tgt: %e",dir,Sqrt(curR2),xToGo));
e1ef49ad 855#endif
0ddbf657 856 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
857 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 858 }
859 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
860 // go via layer to its boundary, applying material correction.
861 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 862 //
e2d2481c 863 return kTRUE;
864 //
865}
866
867//_________________________________________________________________________
868Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
869{
870 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
871 // If check is requested, do this only provided the track has not exited the layer already
872 double xToGo = lr->GetR(dir);
873 if (check) { // do we need to track till the surface of the current layer ?
874 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 875 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
876 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 877 }
e1ef49ad 878#ifdef _ITSU_DEBUG_
b515ad7e 879 AliDebug(2,Form(" dir=%d : from R=%.4f -> R=%.4f",dir,Sqrt(seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY()), xToGo));
e1ef49ad 880#endif
e2d2481c 881 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
882 // go via layer to its boundary, applying material correction.
883 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 884 //
e2d2481c 885 return kTRUE;
886 //
887}
888
889
890//_________________________________________________________________________
891Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
892{
893 // go to the entrance of lr in direction dir, w/o applying material corrections.
894 // If check is requested, do this only provided the track did not reach the layer already
895 double xToGo = lr->GetR(-dir);
896 if (check) { // do we need to track till the surface of the current layer ?
897 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 898 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
899 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 900 }
901 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
902 // go via layer to its boundary, applying material correction.
903 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
904 return kTRUE;
905 //
906}
907
908//_________________________________________________________________________
909Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
910{
911 // go to the entrance of lr in direction dir, w/o applying material corrections.
912 // If check is requested, do this only provided the track did not reach the layer already
913 double xToGo = lr->GetR(-dir);
914 if (check) { // do we need to track till the surface of the current layer ?
915 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 916 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
917 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 918 }
919 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
920 // go via layer to its boundary, applying material correction.
921 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
922 return kTRUE;
923 //
924}
925
926//_________________________________________________________________________
c61e50c3 927Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
32d38de2 928{
929 // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
930 // as well as some aux info
931 double bz = GetBz();
932 AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
c61e50c3 933 seed->GetXYZ(&fTrImpData[kTrXIn]); // lab position at the entrance from above
53870004 934 static AliExternalTrackParam sc; // seed copy for manipulations
943e1898 935 sc = *seed;
32d38de2 936 //
c61e50c3 937 fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
943e1898 938 if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
c61e50c3 939 double dr = lrA->GetDR(); // approximate X dist at the inner radius
943e1898 940 if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
32d38de2 941 // special case: track does not reach inner radius, might be tangential
943e1898 942 double r = sc.GetD(0,0,bz);
32d38de2 943 double x;
bdb92766 944 if (!sc.GetXatLabR(r,x,bz,-1)) return kFALSE;
943e1898 945 dr = Abs(sc.GetX() - x);
bdb92766 946 if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) return kFALSE;
32d38de2 947 }
948 //
c61e50c3 949 fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
943e1898 950 double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
951 double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
ee54014a 952 sgy = Sqrt(sgy)*fCurrTrackCond->GetNSigmaRoadY(ilrA);
953 sgz = Sqrt(sgz)*fCurrTrackCond->GetNSigmaRoadZ(ilrA);
5785a9d8 954 double dphi0 = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]);
955 double phi0 = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
956 if ( dphi0>(0.5*Pi()) ) {
957 // special case of angles around pi
958 dphi0 = Abs(phi0);
959 phi0 += Pi();
960 }
961
962 fTrImpData[kTrPhi0] = phi0;
173b3073 963 fTrImpData[kTrZ0] = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
31e2f5d4 964 dphi0 += sgy/lrA->GetR();
965 fTrImpData[kTrDPhi] = dphi0<PiOver2() ? dphi0 : PiOver2();
173b3073 966 fTrImpData[kTrDZ] = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn]) + sgz;
32d38de2 967 //
968 return kTRUE;
969}
970
b8b59e05 971//________________________________________
972void AliITSUTrackerGlo::ResetSeedsPool()
973{
974 // mark all seeds in the pool as unused
0ddbf657 975 AliDebug(-1,Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
b8b59e05 976 fNFreeSeeds = 0;
977 fSeedsPool.Clear(); // seeds don't allocate memory
978}
979
980
981//________________________________________
982void AliITSUTrackerGlo::MarkSeedFree(AliITSUSeed *sd)
32d38de2 983{
b8b59e05 984 // account that this seed is "deleted"
985 int id = sd->GetPoolID();
986 if (id<0) {
987 AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
988 return;
989 }
990 // AliInfo(Form("%d %p",id, seed));
991 fSeedsPool.RemoveAt(id);
992 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
993 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
32d38de2 994}
f8832015 995
996//_________________________________________________________________________
997Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
998{
999 // Check if the cluster (in tracking frame!) is matching to track.
1000 // The track must be already propagated to sensor tracking frame.
1001 // Returns: kStopSearchOnSensor if the search on given sensor should be stopped,
1002 // kClusterMatching if the cluster is matching
1003 // kClusterMatching otherwise
1004 //
f8832015 1005 const double kTolerX = 5e-4;
c03e4f8a 1006 // AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1007 AliCluster *cl = fCurrLayer->GetCluster(clID);
f8832015 1008 //
c8d1f258 1009 Bool_t goodCl = kFALSE;
0ddbf657 1010 int currLabel = Abs(fCurrESDtrMClb);
716ccba7 1011 //
8b16dbae 1012 if (cl->GetLabel(0)>=0) {
1013 for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}
1014 }
c8d1f258 1015 //
f8832015 1016 if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
c8d1f258 1017 if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
08419930 1018 //
e1ef49ad 1019#ifdef _ITSU_DEBUG_
08419930 1020 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1021 AliDebug(2,Form("Lost good cl on L:%d failed propagation. |ESDtrack#%d (MClb:%d)",lr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1022 track->Print("etp");
15b02d69 1023 cl->Print();
1024 }
e1ef49ad 1025#endif
c8d1f258 1026 return kStopSearchOnSensor; // propagation failed, seedT is intact
1027 }
f8832015 1028 }
3dd9c283 1029 double dy = cl->GetY()-track->GetY();
c8d1f258 1030 double dz = cl->GetZ()-track->GetZ();
0091e9f0 1031 //
f8832015 1032 double dy2 = dy*dy;
1033 double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
ee54014a 1034 fCurrTrackCond->GetNSigmaRoadY(lr)*fCurrTrackCond->GetNSigmaRoadY(lr); // RS TOOPTIMIZE
f8832015 1035 if (dy2>tol2) { // the clusters are sorted in Z(col) then in Y(row).
08419930 1036 //
e1ef49ad 1037#ifdef _ITSU_DEBUG_
08419930 1038 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1039 AliDebug(2,Form("Lost good cl: dy2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",dy2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1040 track->Print("etp");
15b02d69 1041 cl->Print();
53870004 1042 PrintSeedClusters(track);
15b02d69 1043 }
e1ef49ad 1044#endif
cb50e082 1045 // clusters are sorted in increasing Y and the loop goes from last (highers Y) to 1st.
1046 // if the track has too large y for this cluster (dy<0) it will be so for all other clusters of the sensor
1047 if (dy<0) return kStopSearchOnSensor; // No chance that other cluster of this sensor will match (all Y's will be even larger)
f8832015 1048 else return kClusterNotMatching; // Other clusters may match
1049 }
c8d1f258 1050 double dz2 = dz*dz;
f8832015 1051 tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
ee54014a 1052 fCurrTrackCond->GetNSigmaRoadZ(lr)*fCurrTrackCond->GetNSigmaRoadZ(lr); // RS TOOPTIMIZE
c8d1f258 1053 if (dz2>tol2) {
e1ef49ad 1054#ifdef _ITSU_DEBUG_
08419930 1055 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1056 AliDebug(2,Form("Lost good cl on L:%d : dz2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",lr,dz2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1057 track->Print("etp");
15b02d69 1058 cl->Print();
53870004 1059 PrintSeedClusters(track);
15b02d69 1060 }
e1ef49ad 1061#endif
c8d1f258 1062 return kClusterNotMatching; // Other clusters may match
1063 }
f8832015 1064 //
1065 // check chi2
1066 Double_t p[2]={cl->GetY(), cl->GetZ()};
1067 Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
1068 double chi2 = track->GetPredictedChi2(p,cov);
76390254 1069 //
ee54014a 1070 if (chi2>fCurrTrackCond->GetMaxTr2ClChi2(lr)) {
e1ef49ad 1071#ifdef _ITSU_DEBUG_
08419930 1072 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1073 AliDebug(2,Form("Lost good cl on L:%d : Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
ee54014a 1074 lr,chi2,fCurrTrackCond->GetMaxTr2ClChi2(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb));
3dd9c283 1075 track->Print("etp");
1076 cl->Print("");
53870004 1077 PrintSeedClusters(track);
3dd9c283 1078 }
e1ef49ad 1079#endif
c8d1f258 1080 return kClusterNotMatching;
1081 }
f8832015 1082 //
1083 track = NewSeedFromPool(track); // input track will be reused, use its clone for updates
943e1898 1084 if (!track->Update()) {
e1ef49ad 1085#ifdef _ITSU_DEBUG_
08419930 1086 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1087 AliDebug(2,Form("Lost good cluster on L:%d : failed to update",lr));
38997c3c 1088 track->Print("etp");
1089 cl->Print("");
53870004 1090 PrintSeedClusters(track);
38997c3c 1091 }
e1ef49ad 1092#endif
b8b59e05 1093 MarkSeedFree(track);
c8d1f258 1094 return kClusterNotMatching;
1095 }
f8832015 1096 track->SetChi2Cl(chi2);
1097 track->SetLrClusterID(lr,clID);
08419930 1098 //
1099 if (track->GetChi2GloNrm()>fCurrTrackCond->GetMaxChi2GloNrm(lr)) {
e1ef49ad 1100#ifdef _ITSU_DEBUG_
08419930 1101 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
1102 AliDebug(2,Form("Lost good cl on L:%d : Chi2Glo=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
1103 lr,track->GetChi2GloNrm(),fCurrTrackCond->GetMaxChi2GloNrm(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1104 track->Print("etp");
1105 cl->Print("");
1106 PrintSeedClusters(track);
1107 }
e1ef49ad 1108#endif
08419930 1109 MarkSeedFree(track);
1110 return kClusterNotMatching;
1111 }
c03e4f8a 1112 // cl->IncreaseClusterUsage(); // do this only for winners
f8832015 1113 //
716ccba7 1114 track->SetFake(!goodCl);
1115 //
e1ef49ad 1116#ifdef _ITSU_DEBUG_
0ddbf657 1117 AliDebug(2,Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)| ",
716ccba7 1118 goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
e1ef49ad 1119#endif
c8d1f258 1120 //
38997c3c 1121 AddSeedBranch(track);
f8832015 1122 //
1123 return kClusterMatching;
1124}
c8d1f258 1125
1126//_________________________________________________________________________
1127Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
1128{
1129 // check if the seed should not be discarded
1130 const UShort_t kMask = 0xffff;
1131 if (flag==kMissingCluster) {
1132 int lastChecked = seed->GetLayerID();
1133 UShort_t patt = seed->GetHitsPattern();
097e4cdf 1134 if (lastChecked) patt |= (~(kMask<<lastChecked))&fCurrTrackCond->GetAllowedLayers(); // not all layers were checked, complete unchecked ones by potential hits
42c3d4bd 1135 Bool_t seedOK = fCurrTrackCond->CheckPattern(patt);
c8d1f258 1136 return !seedOK;
1137 }
1138 return kTRUE;
1139}
44785f3e 1140
1141//______________________________________________________________________________
1142Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
1143{
1144 // propagate seed to given x applying material correction if requested
1145 const Double_t kEpsilon = 1e-5;
1146 Double_t xpos = seed->GetX();
1147 Int_t dir = (xpos<xToGo) ? 1:-1;
e1ef49ad 1148 Double_t xyz0[3],xyz1[3];
44785f3e 1149 //
8b16dbae 1150 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
f497470c 1151 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
44785f3e 1152 while ( (xToGo-xpos)*dir > kEpsilon){
1153 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1154 Double_t x = xpos+step;
1155 Double_t bz=GetBz(); // getting the local Bz
0ddbf657 1156 if (!seed->PropagateToX(x,bz)) return kFALSE;
8b16dbae 1157 double ds = 0;
1158 if (matCorr || updTime) {
44785f3e 1159 xyz0[0]=xyz1[0]; // global pos at the beginning of step
1160 xyz0[1]=xyz1[1];
1161 xyz0[2]=xyz1[2];
1162 seed->GetXYZ(xyz1); // // global pos at the end of step
8b16dbae 1163 if (matCorr) {
e1ef49ad 1164 Double_t xrho,xx0;
1165 ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
8b16dbae 1166 if (dir>0) xrho = -xrho; // outward should be negative
b2d3bc3a 1167 if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) {AliInfo("Fail2"); seed->Print("etp"); return kFALSE;}
8b16dbae 1168 }
1169 else { // matCorr is not requested but time integral is
1170 double d0 = xyz1[0]-xyz0[0];
1171 double d1 = xyz1[1]-xyz0[1];
1172 double d2 = xyz1[2]-xyz0[2];
1173 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1174 }
44785f3e 1175 }
8b16dbae 1176 if (updTime) seed->AddTimeStep(ds);
44785f3e 1177 xpos = seed->GetX();
1178 }
1179 return kTRUE;
1180}
716ccba7 1181
1182//______________________________________________________________________________
3e4e3c23 1183Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
1184{
1185 // propagate seed to given x applying material correction if requested
1186 const Double_t kEpsilon = 1e-5;
1187 Double_t xpos = seed->GetX();
1188 Int_t dir = (xpos<xToGo) ? 1:-1;
e1ef49ad 1189 Double_t xyz0[3],xyz1[3];
3e4e3c23 1190 //
e1ef49ad 1191#ifdef _ITSU_DEBUG_
b515ad7e 1192 if (AliDebugLevelClass()>1) {
1193 AliDebug(2,Form("Before Propagate to X=%f with M=%.3f MaxStep=%.4f MatCorr=%d",xToGo,mass,maxStep,matCorr));
1194 seed->AliExternalTrackParam::Print();
1195 }
e1ef49ad 1196#endif
3e4e3c23 1197 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
1198 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
1199 while ( (xToGo-xpos)*dir > kEpsilon){
1200 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1201 Double_t x = xpos+step;
1202 Double_t bz=GetBz(); // getting the local Bz
1203 if (!seed->PropagateTo(x,bz)) return kFALSE;
1204 double ds = 0;
1205 if (matCorr || updTime) {
1206 xyz0[0]=xyz1[0]; // global pos at the beginning of step
1207 xyz0[1]=xyz1[1];
1208 xyz0[2]=xyz1[2];
1209 seed->GetXYZ(xyz1); // // global pos at the end of step
1210 //
1211 if (matCorr) {
e1ef49ad 1212 Double_t xrho,xx0;
1213 ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
3e4e3c23 1214 if (dir>0) xrho = -xrho; // outward should be negative
1215 if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
3e4e3c23 1216 }
1217 else { // matCorr is not requested but time integral is
1218 double d0 = xyz1[0]-xyz0[0];
1219 double d1 = xyz1[1]-xyz0[1];
1220 double d2 = xyz1[2]-xyz0[2];
1221 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1222 }
1223 }
1224 if (updTime) seed->AddTimeStep(ds);
1225 //
1226 xpos = seed->GetX();
1227 }
e1ef49ad 1228#ifdef _ITSU_DEBUG_
b515ad7e 1229 if (AliDebugLevelClass()>1) {
1230 AliDebug(2,Form("After Propagate to X=%f",xToGo));
1231 seed->AliExternalTrackParam::Print();
1232 }
e1ef49ad 1233#endif
3e4e3c23 1234 return kTRUE;
1235}
1236
1237//______________________________________________________________________________
9cdcba2c 1238Bool_t AliITSUTrackerGlo::FinalizeHypothesis(AliITSUTrackHyp* hyp)
1239{
1240 // finalize single TPC track hypothesis
1241 if (!hyp || hyp->GetESDTrack()->IsOn(AliESDtrack::kITSin)) return kFALSE;
1242 AliITSUSeed* winner = 0;
1243 fCurrHyp = hyp;
70cb7fe4 1244 fCurrMass = hyp->GetMass();
d99994b2 1245 if (!(winner=hyp->DefineWinner(fCurrTrackCond->GetActiveLrInner()))) return kFALSE;
9cdcba2c 1246 CookMCLabel(hyp);
1247 //
08419930 1248#ifdef _ITSU_TUNING_MODE_ // fill tuning histos
1249 TObjArray* dest = hyp->GetLabel()>=0 ? fCHistoArrCorr : fCHistoArrFake;
1250 if (dest) {
1251 AliITSUSeed* sd = winner;
1252 double htrPt = hyp->Pt();
1253 do {
1254 int clID,lrID;
1255 if ( (clID=sd->GetLrCluster(lrID))<0 ) continue;
781c79ac 1256 ((TH2*)dest->At( GetHistoID(lrID,kHChi2Nrm ,fCurrPassID) ))->Fill(htrPt,sd->GetChi2GloNrm());
1257 ((TH2*)dest->At( GetHistoID(lrID,kHBestInBranch,fCurrPassID) ))->Fill(htrPt,sd->GetOrdBranch());
1258 ((TH2*)dest->At( GetHistoID(lrID,kHBestInCand ,fCurrPassID) ))->Fill(htrPt,sd->GetOrdCand());
08419930 1259 //
1260 if (dest==fCHistoArrFake && !sd->IsFake()) continue; // for the fake seeds fill only fake clusters part
1261 //
781c79ac 1262 ((TH2*)dest->At( GetHistoID(lrID,kHResY ,fCurrPassID) ))->Fill(htrPt,sd->GetResidY());
1263 ((TH2*)dest->At( GetHistoID(lrID,kHResZ ,fCurrPassID) ))->Fill(htrPt,sd->GetResidZ());
1264 ((TH2*)dest->At( GetHistoID(lrID,kHResYP ,fCurrPassID) ))->Fill(htrPt,sd->GetPullY());
1265 ((TH2*)dest->At( GetHistoID(lrID,kHResZP ,fCurrPassID) ))->Fill(htrPt,sd->GetPullZ());
1266 ((TH2*)dest->At( GetHistoID(lrID,kHChi2Cl ,fCurrPassID) ))->Fill(htrPt,sd->GetChi2Cl());
08419930 1267 //
1268 } while((sd=(AliITSUSeed*)sd->GetParent()));
1269 //
781c79ac 1270 ((TH2*)dest->At( GetHistoID(-1,kHChiMatch,fCurrPassID) ))->Fill(htrPt,winner->GetChi2ITSTPC());
1271 ((TH2*)dest->At( GetHistoID(-1,kHChiITSSA,fCurrPassID) ))->Fill(htrPt,winner->GetChi2ITSSA());
9cdcba2c 1272 }
9cdcba2c 1273 //
08419930 1274#endif
1275 //
1276 CheckClusterSharingConflicts(hyp);
1277 return hyp->GetWinner(); // winner might change of disappear after resolving conflicts
1278 //
9cdcba2c 1279}
1280
1281//______________________________________________________________________________
3e4e3c23 1282void AliITSUTrackerGlo::FinalizeHypotheses()
1283{
1284 // select winner for each hypothesis, remove cl. sharing conflicts
9cdcba2c 1285 AliCodeTimerAuto("",0);
3e4e3c23 1286 //
d99994b2 1287 int* esdTrackIndex = fESDIndex.GetArray();
1288 for (int ih=0;ih<fNTracksESD;ih++) FinalizeHypothesis(GetTrackHyp(esdTrackIndex[ih])); // finlize in decreasing pt order
c03e4f8a 1289 //
d99994b2 1290 for (int ih=fNTracksESD;ih--;) UpdateESDTrack(GetTrackHyp(esdTrackIndex[ih]),AliESDtrack::kITSin);
08419930 1291 //
1292}
bdb92766 1293
08419930 1294//______________________________________________________________________________
1295void AliITSUTrackerGlo::CheckClusterSharingConflicts(AliITSUTrackHyp* hyp)
1296{
1297 // remove eventual cluster sharing conflicts
1298 AliITSUSeed* winner = 0;
1299 if (!(winner=hyp->GetWinner())) return;
1300 UShort_t idH = (UShort_t)hyp->GetUniqueID()+1;
1301 int lrID,clID;
d99994b2 1302 int inLr = fCurrTrackCond->GetActiveLrInner();
08419930 1303 AliITSUSeed* winSD=winner;
1304 do {
1305 if ( (clID=winSD->GetLrCluster(lrID))<0 ) continue;
1306 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1307 Int_t refID = cl->GetRecoInfo(); // was it already referred by some track (with refID-1 if >0)
1308 if (!refID) {cl->SetRecoInfo(idH); continue;} // if not, refer from track to cluster
1309 if (refID==idH) continue; // ignore reference to itself
1310 // the cluster is already used by other track, need to resolve the conflict
1311 AliITSUTrackHyp* hypC = GetTrackHyp(refID-1); // competitor
1312 // printf("Ref to %d (%p) from %d Cl %d %d\n",refID-1,hypC, idH-1,clID,lrID);
1313 AliITSUSeed* winnerC = hypC->GetWinner();
1314 if (!winnerC) {
1315 printf("Missing winner, ERROR\n");
1316 continue;
1317 }
6091ad4e 1318 //
e1ef49ad 1319#ifdef _ITSU_DEBUG_
08419930 1320 if (AliDebugLevelClass()>1) { //------------------------ DEBUG -------------------
1321 AliInfo(Form("Shared cl#%4d on Lr:%d: Hyp%3d/Q:%6.3f vs Hyp%3d/Q:%6.3f",
1322 clID,lrID,idH-1,winner->GetQualityVar(),refID-1,winnerC->GetQualityVar()));
1323 // dump winner of hyp
1324 printf("#%4d Pt:%.3f Chi:%6.3f/%6.3f/%6.3f Ncl:%d MCits%+5d MCtpc:%+5d |",
1325 idH-1,winner->Pt(),winner->GetChi2GloNrm(),winner->GetChi2ITSTPC(),winner->GetChi2ITSSA(),
1326 winner->GetNLayersHit(),hyp->GetITSLabel(),hyp->GetESDTrack()->GetTPCLabel());
1327 int prevL=-1;
1328 AliITSUSeed* sd = winner;
1329 do {
1330 int lrs;
1331 int clIDt = sd->GetLrCluster(lrs);
1332 if (clIDt<0) continue;
1333 while( lrs>++prevL ) printf("%4s ","----");
1334 printf("%4d (%5.1f)",clIDt,sd->GetChi2Cl());
1335 } while ((sd=(AliITSUSeed*)sd->GetParent()));
1336 printf("|\n");
1337 // dump winner of hypC
1338 printf("#%4d Pt:%.3f Chi:%6.3f/%6.3f/%6.3f Ncl:%d MCits%+5d MCtpc:%+5d |",
1339 refID-1,winnerC->Pt(),winnerC->GetChi2GloNrm(),winnerC->GetChi2ITSTPC(),winnerC->GetChi2ITSSA(),
1340 winnerC->GetNLayersHit(),hypC->GetITSLabel(),hypC->GetESDTrack()->GetTPCLabel());
1341 prevL=-1;
1342 sd = winnerC;
6091ad4e 1343 do {
08419930 1344 int lrs;
1345 int clIDt = sd->GetLrCluster(lrs);
1346 if (clIDt<0) continue;
1347 while( lrs>++prevL ) printf("%4s ","----");
1348 printf("%4d (%5.1f)",clIDt,sd->GetChi2Cl());
1349 } while ((sd=(AliITSUSeed*)sd->GetParent()));
1350 printf("|\n");
1351 } //-------------------------------------------------- DEBUG -------------------
e1ef49ad 1352#endif
08419930 1353 //
1354 if (winner->GetQualityVar()<winnerC->GetQualityVar()) { // current tracks is better then competitor track
1355 FlagSeedClusters(winnerC,kFALSE,refID); // unflag cluster usage by loser
1356 cl->SetRecoInfo(idH);
1357 //
1358 //if ( hypC->GetLabel()>=0) AliInfo(Form("KILLING CORRECT: %d",refID-1));
1359 winnerC->Kill();
d99994b2 1360 if (hypC->DefineWinner(inLr)) { // find new winner instead of suppressed candidate
08419930 1361 CookMCLabel(hypC);
1362 CheckClusterSharingConflicts(hypC); // and check its sharing conflicts again
1363 if (winner->IsKilled()) break; // the current winner might have been killed during check of new winner of hypC!
1364 }
6091ad4e 1365 }
08419930 1366 else { // competitor hypC is better than the hyp
1367 FlagSeedClusters(winner,kFALSE,idH); // unflag cluster usage by loser
1368 winner->Kill();
1369 //if ( hyp->GetLabel()>=0) AliInfo(Form("KILLING CORRECT: %d",idH-1));
d99994b2 1370 if (hyp->DefineWinner(inLr)) {
08419930 1371 CookMCLabel(hyp);
1372 CheckClusterSharingConflicts(hyp);
9cdcba2c 1373 }
08419930 1374 break;
9cdcba2c 1375 }
08419930 1376 //
1377 } while ((winSD=(AliITSUSeed*)winSD->GetParent()));
1378 //
1379 return;
3e4e3c23 1380}
c51c3124 1381
1382//______________________________________________________________________________
68a0f687 1383void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp,Int_t flag)
c51c3124 1384{
1385 // update ESD track with current best hypothesis
08419930 1386 if (!hyp) return;
c51c3124 1387 AliESDtrack* esdTr = hyp->GetESDTrack();
d99994b2 1388 if (!esdTr || esdTr->IsOn(flag)) return;
c51c3124 1389 AliITSUSeed* win = hyp->GetWinner();
08419930 1390 if (!win || win->IsKilled()) return;
70cb7fe4 1391 double chiSav;
1392 //
68a0f687 1393 switch (flag) {
1394 case AliESDtrack::kITSin:
1395 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
08419930 1396 fCountITSin++;
dc443a34 1397 // set fakes cluster info
1398 {
1399 UShort_t clfake = 0;
c7e63626 1400 Int_t clSplit = 0;
dc443a34 1401 AliITSUSeed* sd = win;
c7e63626 1402 int ip=0;
dc443a34 1403 do {
c7e63626 1404 int lr, clID = sd->GetLrCluster(lr);
1405 if (sd->IsFake()) clfake |= 0x1<<lr;
1406 if (clID>=0) {
1407 esdTr->SetITSModuleIndex(ip++, sd->GetLrClusterID());
1408 AliITSUClusterPix *cl = (AliITSUClusterPix*)fITS->GetLayerActive(lr)->GetCluster(clID);
1409#ifdef _ITSU_TUNING_MODE_
1410 if (cl->IsSplit()) clSplit |= 0x1<<lr;
1411#endif
1412 }
dc443a34 1413 } while ((sd=(AliITSUSeed*)sd->GetParent()));
0e40ce66 1414 //
1415 // RS: Temporary set special flag for tracks from the afterburner
1416 if (fCurrPassID>0) clfake |= 0x1<<7;
1417 //
dc443a34 1418 esdTr->SetITSSharedMap(clfake);
c7e63626 1419 esdTr->SetITSModuleIndex(10,clSplit);
dc443a34 1420 }
097e4cdf 1421 // TEMPORARY: store iteration id
1422 esdTr->SetITSModuleIndex(11,fCurrPassID);
68a0f687 1423 break;
1424 //
1425 case AliESDtrack::kITSout:
70cb7fe4 1426 // here the stored chi2 will correspond to backward ITS-SA fit
68a0f687 1427 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
08419930 1428 fCountITSout++;
68a0f687 1429 // TODO: avoid setting friend
1430 break;
1431 //
1432 case AliESDtrack::kITSrefit:
70cb7fe4 1433
1434 // at the moment fill as a chi2 the TPC-ITS matching chi2
1435 chiSav = hyp->GetChi2();
1436 hyp->SetChi2(win->GetChi2ITSTPC());
68a0f687 1437 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
70cb7fe4 1438 hyp->SetChi2(chiSav);
08419930 1439 fCountITSrefit++;
68a0f687 1440 // TODO: avoid setting cluster info
1441 break;
1442 default:
1443 AliFatal(Form("Unknown flag %d",flag));
1444 }
c51c3124 1445 //
bdb92766 1446 esdTr->SetITSLabel(hyp->GetITSLabel());
c51c3124 1447 // transfer module indices
1448 // TODO
1449}
1450
1451//______________________________________________________________________________
7ce7f445 1452Double_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest, Int_t &nclFit, Int_t stopCond)
c51c3124 1453{
9cdcba2c 1454 // refit track till radius rDest. The cluster,mass info is taken from fCurrHyp (and its winner seed)
1455 // if stopCond<0 : propagate till last cluster then stop
1456 // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail
1457 // if stopCond>0 : rDest must be reached
8b16dbae 1458 //
c51c3124 1459 double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
1460 int dir,lrStart,lrStop;
1461 //
8b16dbae 1462 dir = rCurr<rDest ? 1 : -1;
1463 lrStart = fITS->FindFirstLayerID(rCurr,dir);
70cb7fe4 1464 lrStop = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
1465 //
e1ef49ad 1466#ifdef _ITSU_DEBUG_
15b02d69 1467 if (AliDebugLevelClass()>2) {
1468 printf("Refit %d: Lr:%d (%f) -> Lr:%d (%f)\n",dir,lrStart,rCurr, lrStop,rDest);
9cdcba2c 1469 printf("Before refit: "); trc->Print();
15b02d69 1470 }
e1ef49ad 1471#endif
8b16dbae 1472 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));
1473 //
9cdcba2c 1474 Int_t clInfo[2*AliITSUAux::kMaxLayers];
1475 Int_t nCl = fCurrHyp->FetchClusterInfo(clInfo);
1476 fCurrMass = fCurrHyp->GetMass();
7ce7f445 1477 AliITSUTrackHyp tmpTr(*(AliKalmanTrack*)trc);
9cdcba2c 1478 double chi2 = 0;
5e3587be 1479 int iclLr[2],nclLr;
1480 nclFit = 0;
c51c3124 1481 //
c7e63626 1482 int lrStop1 = lrStop+dir;
70cb7fe4 1483 for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
c51c3124 1484 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
1485 if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
8b16dbae 1486 int ilrA2,ilrA = lr->GetActiveID();
1487 // passive layer or active w/o hits will be traversed on the way to next cluster
9cdcba2c 1488 if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) continue;
8b16dbae 1489 //
8b16dbae 1490 // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
1491 nclLr=0;
1492 if (dir>0) { // clusters are stored in increasing radius order
9cdcba2c 1493 iclLr[nclLr++]=clInfo[ilrA2++];
1494 if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2];
8b16dbae 1495 }
1496 else {
9cdcba2c 1497 if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1];
1498 iclLr[nclLr++]=clInfo[ilrA2];
8b16dbae 1499 }
1500 //
b515ad7e 1501 Bool_t transportedToLayer = kFALSE;
8b16dbae 1502 for (int icl=0;icl<nclLr;icl++) {
1503 AliITSUClusterPix* clus = (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
1504 AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
e7d83d38 1505 if (!tmpTr.Rotate(sens->GetPhiTF())) {
e1ef49ad 1506#ifdef _ITSU_DEBUG_
fb0eae48 1507 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1508 AliDebug(2,Form("Failed on rotate to %f | ESDtrack#%d (MClb:%d)",sens->GetPhiTF(),trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1509#endif
9cdcba2c 1510 return -1;
e7d83d38 1511 }
b515ad7e 1512 //
1513 double xClus = sens->GetXTF()+clus->GetX();
1514 if (!transportedToLayer) {
1515 if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
e1ef49ad 1516#ifdef _ITSU_DEBUG_
846c381e 1517 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1518 AliDebug(2,Form("Failed to transport %d -> %d | ESDtrack#%d (MClb:%d)\n",lrStart,ilr,trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1519#endif
9cdcba2c 1520 return -1; // go to the entrance to the layer
b515ad7e 1521 }
1522 lrStart = ilr;
1523 transportedToLayer = kTRUE;
1524 }
1525 //
e1ef49ad 1526#ifdef _ITSU_DEBUG_
b515ad7e 1527 if (AliDebugLevelClass()>1) {
1528 AliDebug(2,Form("Propagate to cl:%d on lr %d Need to go %.4f -> %.4f",icl,ilrA,tmpTr.GetX(),xClus));
b515ad7e 1529 }
e1ef49ad 1530#endif
b515ad7e 1531 //
1532 if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) {
e1ef49ad 1533#ifdef _ITSU_DEBUG_
846c381e 1534 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1535 AliDebug(2,Form("Failed on propagate to %f (dir=%d) | ESDtrack#%d (MClb:%d)",xClus,dir,trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1536#endif
9cdcba2c 1537 return -1;
e7d83d38 1538 }
76390254 1539 //
9cdcba2c 1540 Double_t p[2]={clus->GetY(), clus->GetZ()};
1541 Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()};
1542 double chi2cl = tmpTr.GetPredictedChi2(p,cov);
1543 chi2 += chi2cl;
1544 //
08419930 1545#ifdef _ITSU_TUNING_MODE_
1546 TObjArray* dest = trc->GetLabel()>=0 ? fCHistoArrCorr : fCHistoArrFake;
781c79ac 1547 if (dest && fTrackPhaseID>kClus2Tracks) {
08419930 1548 //
1549 double htrPt = tmpTr.Pt();
9cdcba2c 1550 double dy = p[0]-tmpTr.GetY();
1551 double dz = p[1]-tmpTr.GetZ();
781c79ac 1552 ((TH2*)dest->At( GetHistoID(ilrA,kHResY,0,fTrackPhaseID) ))->Fill(htrPt,dy);
1553 ((TH2*)dest->At( GetHistoID(ilrA,kHResZ,0,fTrackPhaseID) ))->Fill(htrPt,dz);
76390254 1554 double errY = tmpTr.GetSigmaY2();
1555 double errZ = tmpTr.GetSigmaZ2();
781c79ac 1556 if (errY>0) ((TH2*)dest->At( GetHistoID(ilrA,kHResYP,0,fTrackPhaseID) ))->Fill(htrPt,dy/Sqrt(errY));
1557 if (errZ>0) ((TH2*)dest->At( GetHistoID(ilrA,kHResZP,0,fTrackPhaseID) ))->Fill(htrPt,dz/Sqrt(errZ));
1558 ((TH2*)dest->At( GetHistoID(ilrA,kHChi2Cl,0,fTrackPhaseID) ))->Fill(htrPt,chi2cl);
76390254 1559 }
1560#endif
9cdcba2c 1561 //
1562 if ( !tmpTr.Update(p,cov) ) {
e1ef49ad 1563#ifdef _ITSU_DEBUG_
846c381e 1564 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1565 AliDebug(2,Form("Failed on Update | ESDtrack#%d (MClb:%d)",trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1566#endif
9cdcba2c 1567 return -1;
e7d83d38 1568 }
e1ef49ad 1569#ifdef _ITSU_DEBUG_
b515ad7e 1570 if (AliDebugLevelClass()>1) {
9cdcba2c 1571 printf("AfterRefit: "); tmpTr.Print();
1572 }
e1ef49ad 1573#endif
5e3587be 1574 if (++nclFit==nCl && stopCond<0) {
7ce7f445 1575 *trc = (AliKalmanTrack&)tmpTr;
9cdcba2c 1576 return chi2; // it was requested to not propagate after last update
b515ad7e 1577 }
8b16dbae 1578 }
c51c3124 1579 //
1580 }
8b16dbae 1581 // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
1582 // Still, try to go as close as possible to rDest.
1583 //
70cb7fe4 1584 if (lrStart!=lrStop) {
e7d83d38 1585 if (!TransportToLayer(&tmpTr,lrStart,lrStop)) {
e1ef49ad 1586#ifdef _ITSU_DEBUG_
2947434b 1587 AliDebug(1,Form("Failed on TransportToLayer %d->%d | ESDtrack#%d (MClb:%d), X=%f",lrStart,lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb,tmpTr.GetX()));
e1ef49ad 1588#endif
2947434b 1589 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
e7d83d38 1590 }
1591 if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) {
e1ef49ad 1592#ifdef _ITSU_DEBUG_
2947434b 1593 AliDebug(1,Form("Failed on GoToExitFromLayer %d | ESDtrack#%d (MClb:%d), X=%f",lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb,tmpTr.GetX()));
e1ef49ad 1594#endif
2947434b 1595 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
9cdcba2c 1596 }
70cb7fe4 1597 }
1598 // go to the destination radius. Note that here we don't select direction to avoid precision problems
1599 if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) {
2947434b 1600 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
70cb7fe4 1601 }
7ce7f445 1602 *trc = (AliKalmanTrack&)tmpTr;
e1ef49ad 1603#ifdef _ITSU_DEBUG_
15b02d69 1604 if (AliDebugLevelClass()>2) {
9cdcba2c 1605 printf("After refit (now at lr %d): ",lrStop); trc->Print();
15b02d69 1606 }
e1ef49ad 1607#endif
9cdcba2c 1608 return chi2;
c51c3124 1609}
dd2117a2 1610
1611//__________________________________________________________________
1612void AliITSUTrackerGlo::CookMCLabel(AliITSUTrackHyp* hyp)
1613{
1614 // build MC label
1615 //
1616 const int kMaxLbPerCl = 3;
1617 int lbID[kMaxLayers*6],lbStat[kMaxLayers*6];
1618 Int_t lr,nLab=0,nCl=0;
1619 AliITSUSeed *seed = hyp->GetWinner();
1620 while(seed) {
1621 int clID = seed->GetLrCluster(lr);
1622 if (clID>=0) {
1623 AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1624 nCl++;
1625 for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
1626 int trLb = cl->GetLabel(imc);
704d9398 1627 if (trLb<0) break;
dd2117a2 1628 // search this mc track in already accounted ones
1629 int iLab;
1630 for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
1631 if (iLab<nLab) lbStat[iLab]++;
1632 else {
1633 lbID[nLab] = trLb;
1634 lbStat[nLab++] = 1;
1635 }
1636 } // loop over given cluster's labels
1637 }
1638 seed = (AliITSUSeed*)seed->GetParent();
1639 } // loop over clusters
1640 //
5785a9d8 1641 AliESDtrack* esdTr = hyp->GetESDTrack();
1642 int tpcLab = esdTr ? Abs(esdTr->GetTPCLabel()) : -kDummyLabel;
41343470 1643 if (nCl && nLab) {
dd2117a2 1644 int maxLab=0,nTPCok=0;
dd2117a2 1645 for (int ilb=nLab;ilb--;) {
1646 int st = lbStat[ilb];
1647 if (lbStat[maxLab]<st) maxLab=ilb;
1648 if (tpcLab==lbID[ilb]) nTPCok += st;
1649 }
1650 hyp->SetFakeRatio(1.-float(nTPCok)/nCl);
1651 hyp->SetLabel( nTPCok==nCl ? tpcLab : -tpcLab);
bdb92766 1652 hyp->SetITSLabel( lbStat[maxLab]==nCl ? lbID[maxLab] : -lbID[maxLab]); // winner label
dd2117a2 1653 return;
1654 }
1655 //
1656 hyp->SetFakeRatio(-1.);
5785a9d8 1657 hyp->SetLabel( -tpcLab );
dd2117a2 1658 hyp->SetITSLabel( kDummyLabel );
1659}
38997c3c 1660
1661//__________________________________________________________________
1662Bool_t AliITSUTrackerGlo::AddSeedBranch(AliITSUSeed* seed)
1663{
1664 // add new prolongation branch for the seed to temporary array. The seeds are sorted according to Compare f-n
6cd80116 1665 if (fNCandidatesAdded+fNBranchesAdded>=fLayerMaxCandidates ) {
1666 AliInfo(Form("Number of candidates at layer %d for current seed reached %d, increasing buffer",fCurrActLrID,fLayerMaxCandidates));
1667 fLayerMaxCandidates = 2*(fLayerMaxCandidates+1);
1668 AliITSUSeed** tmpArr = fLayerCandidates;
1669 fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
1670 memcpy(fLayerCandidates,tmpArr,(fNCandidatesAdded+fNBranchesAdded)*sizeof(AliITSUSeed*));
1671 delete tmpArr; // delete only array, not objects
38997c3c 1672 }
1673 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded]; // note: fNCandidatesAdded is incremented after adding all branches of current seed
1674 int slot=fNBranchesAdded++;
1675 for (int slotF=slot;slotF--;) { // slotF is always slot-1
1676 AliITSUSeed* si = branches[slotF];
1677 if (si->Compare(seed)<0) break; // found the last seed with better quality
1678 // otherwise, shift the worse seed to the next slot
1679 branches[slot] = si;
1680 slot = slotF; // slot should be slotF+1
1681 }
1682 // if needed, move worse seeds
1683 branches[slot] = seed;
1684 return kTRUE;
1685 //
1686}
1687
1688//__________________________________________________________________
1689void AliITSUTrackerGlo::ValidateAllowedBranches(Int_t acceptMax)
1690{
1691 // keep allowed number of branches for current seed and disable extras
70cb7fe4 1692 AliCodeTimerAuto("",0);
38997c3c 1693 int nb = Min(fNBranchesAdded,acceptMax);
1694 // if (nb<fNBranchesAdded) printf("ValidateAllowedBranches: %d of %d (%d) on lr %d\n",nb,fNBranchesAdded,fNCandidatesAdded,ilr);
1695 // disable unused branches
1696 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded];
08419930 1697 //
1698#ifdef _ITSU_TUNING_MODE_
1699 for (int ib=0;ib<fNBranchesAdded;ib++) branches[ib]->SetOrdBranch(ib);
c03e4f8a 1700#endif
bdb92766 1701 //
08419930 1702 for (int ib=nb;ib<fNBranchesAdded;ib++) MarkSeedFree(branches[ib]);
38997c3c 1703 fNCandidatesAdded += nb; // update total candidates counter
1704 fNBranchesAdded = 0; // reset branches counter
1705 //
1706}
1707
1708//__________________________________________________________________
1709void AliITSUTrackerGlo::ValidateAllowedCandidates(Int_t ilr, Int_t acceptMax)
1710{
1711 // transfer allowed number of branches to hypothesis container
1712 //
9cdcba2c 1713 AliCodeTimerAuto("",0);
38997c3c 1714 // sort candidates in increasing order of chi2
6cd80116 1715 static int lastSize = 0;
1716 static int *index = 0;
1717 static float *chi2 = 0;
1718 if (fLayerMaxCandidates>lastSize) {
1719 lastSize = fLayerMaxCandidates;
1720 delete[] index;
1721 delete[] chi2;
1722 index = new int[lastSize];
1723 chi2 = new float[lastSize];
1724 }
38997c3c 1725 for (int i=fNCandidatesAdded;i--;) chi2[i] = fLayerCandidates[i]->GetChi2GloNrm();
1726 Sort(fNCandidatesAdded,chi2,index,kFALSE);
1727 //
9cdcba2c 1728 int nacc=0,nb=0;
097e4cdf 1729 if (ilr>fCurrTrackCond->GetActiveLrInner()) { // just take 1st acceptMax candidates
9cdcba2c 1730 nb = Min(fNCandidatesAdded,acceptMax);
1731 for (int ib=0;ib<nb;ib++) AddProlongationHypothesis(fLayerCandidates[index[ib]],ilr);
1732 }
1733 else { // for innermost layer test ITS SA backward chi2 and matching to TPC
1734 AliITSUSeed* wn0 = fCurrHyp->GetWinner(); // in principle, should be NULL
1735 for (nb=0;nb<fNCandidatesAdded;nb++) {
1736 AliITSUSeed* sd = fLayerCandidates[index[nb]];
1737 fCurrHyp->SetWinner(sd);
1738 if (!CheckBackwardMatching(sd)) MarkSeedFree(sd); // discard bad backward tracks
1739 else {
1740 AddProlongationHypothesis(sd,ilr);
1741 if (++nacc==acceptMax) {nb++; break;} // the rest will be discarded
1742 }
1743 }
1744 fCurrHyp->SetWinner(wn0); // restore original winner (NULL?)
1745 }
38997c3c 1746 // disable unused candidates
08419930 1747 for (int ib=nb;ib<fNCandidatesAdded;ib++) MarkSeedFree(fLayerCandidates[index[ib]]);
38997c3c 1748 fNCandidatesAdded = 0; // reset candidates counter
1749 //
1750}
1751
53870004 1752//__________________________________________________________________
9cdcba2c 1753Bool_t AliITSUTrackerGlo::CheckBackwardMatching(AliITSUSeed* seed)
1754{
1755 // check seed backward propagation chi2 and matching to TPC
1756 double bz0 = GetBz();
70cb7fe4 1757 double rDest = fITS->GetRITSTPCRef(); // reference radius for comparison
7ce7f445 1758 AliITSUTrackHyp trback;
1759 trback.AliExternalTrackParam::operator=(*seed);
9cdcba2c 1760 trback.ResetCovariance(10000);
5e3587be 1761 int nclFit = 0;
1762 double chi2sa = RefitTrack(&trback,rDest,nclFit,1);
9cdcba2c 1763 if (chi2sa<0) return kFALSE;
5e3587be 1764 int ndf = nclFit*2-5;
9cdcba2c 1765 if (ndf>0) chi2sa /= ndf;
5e3587be 1766 if (chi2sa>fCurrTrackCond->GetMaxITSSAChi2(nclFit)) return kFALSE;
9cdcba2c 1767 //
1768 // relate to TPC track at outer layer
1769 AliExternalTrackParam* tpcSeed = fCurrHyp->GetTPCSeed();
1770 if (!trback.Rotate(tpcSeed->GetAlpha()) || !trback.PropagateParamOnlyTo(tpcSeed->GetX(),bz0)) {
e1ef49ad 1771#ifdef _ITSU_DEBUG_
9cdcba2c 1772 if (AliDebugLevelClass()>+1 && !seed->ContainsFake()) {
1773 AliInfo(Form("Failed to propagate seed to TPC track @ X:%.3f Alpha:%.3f",tpcSeed->GetX(),tpcSeed->GetAlpha()));
1774 seed->Print("etp");
1775 trback.Print();
1776 }
e1ef49ad 1777#endif
9cdcba2c 1778 return kFALSE;
1779 }
1780 double chi2Match = trback.GetPredictedChi2(tpcSeed)/5;
1781 if (chi2Match>fCurrTrackCond->GetMaxITSTPCMatchChi2()) return kFALSE;
1782 //
1783 seed->SetChi2ITSSA(chi2sa);
1784 seed->SetChi2ITSTPC(chi2Match);
1785 return kTRUE;
1786}
1787
1788//__________________________________________________________________
6cd80116 1789void AliITSUTrackerGlo::SaveReducedHypothesesTree(AliITSUTrackHyp* dest)
1790{
1791 // remove those hypothesis seeds which dont lead to candidates at final layer
1792 //
1793 // 1st, flag the seeds to save
0e40ce66 1794 int lr0 = fCurrTrackCond->GetActiveLrInner();
6cd80116 1795 for (int isd=0;isd<fCurrHyp->GetNSeeds(lr0);isd++) {
1796 AliITSUSeed* seed = fCurrHyp->RemoveSeed(lr0,isd);
1797 if (!seed) continue;
1798 seed->FlagTree(AliITSUSeed::kSave);
1799 dest->AddSeed(seed,lr0);
1800 }
0e40ce66 1801 for (int ilr=0;ilr<fNLrActive;ilr++) {
6cd80116 1802 int nsd = fCurrHyp->GetNSeeds(ilr);
1803 for (int isd=0;isd<nsd;isd++) {
1804 AliITSUSeed* seed = fCurrHyp->RemoveSeed(ilr,isd);
1805 if (!seed) continue; // already discarded or saved
1806 if (seed->IsSaved()) dest->AddSeed(seed,ilr);
1807 else MarkSeedFree(seed);
1808 }
1809 }
1810 //
1811 // AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
1812}
1813
1814//__________________________________________________________________
0e40ce66 1815void AliITSUTrackerGlo::CleanHypothesis(AliITSUTrackHyp* hyp)
1816{
1817 // clean hypothesis
1818 hyp->SetWinner(0);
1819 for (int ilr=0;ilr<fNLrActive;ilr++) {
1820 int nsd = hyp->GetNSeeds(ilr);
1821 for (int isd=0;isd<nsd;isd++) {
1822 AliITSUSeed* seed = hyp->RemoveSeed(ilr,isd);
1823 if (!seed) continue; // already discarded or saved
1824 MarkSeedFree(seed);
1825 }
1826 }
1827}
1828
1829//__________________________________________________________________
53870004 1830void AliITSUTrackerGlo::FlagSplitClusters()
1831{
1832 // set special bit on split clusters using MC info
6cd80116 1833 for (int ilr=fNLrActive;ilr--;) {
53870004 1834 int nsplit=0;
1835 AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
1836 for (int isn=lr->GetNSensors();isn--;) {
1837 AliITSURecoSens* sens = lr->GetSensor(isn);
1838 int nCl = sens->GetNClusters();
1839 if (!nCl) continue;
1840 int cl0 = sens->GetFirstClusterId();
1841 for (int icl=nCl;icl--;) {
1842 AliITSUClusterPix *cl = (AliITSUClusterPix*)lr->GetCluster(cl0+icl);
1843 for (int icl1=icl;icl1--;) {
1844 AliITSUClusterPix *cl1 = (AliITSUClusterPix*)lr->GetCluster(cl0+icl1);
1845 if (cl->HasCommonTrack(cl1)) {
1846 if (!cl->IsSplit()) nsplit++;
1847 if (!cl1->IsSplit()) nsplit++;
1848 cl->SetSplit();
1849 cl1->SetSplit();
1850 }
1851 }
1852 }
1853 }
1854 AliInfo(Form("%4d out of %4d clusters are split",nsplit,lr->GetNClusters()));
1855 }
1856 //
1857}
1858
1859//__________________________________________________________________
1860Bool_t AliITSUTrackerGlo::ContainsSplitCluster(const AliITSUSeed* seed, Int_t maxSize)
1861{
1862 // check if the seed contains split cluster with size < maxSize
1863 int lrID,clID;
1864 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1865 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1866 if (cl->IsSplit() && cl->GetNPix()<maxSize ) return kTRUE;
1867 }
1868 return seed->GetParent() ? ContainsSplitCluster((AliITSUSeed*)seed->GetParent(),maxSize) : kFALSE;
1869 //
1870}
1871
1872//__________________________________________________________________
1873void AliITSUTrackerGlo::PrintSeedClusters(const AliITSUSeed* seed, Option_t* option)
1874{
1875 // print seeds clusters
1876 int lrID,clID;
1877 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1878 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1879 cl->Print(option);
1880 }
1881 if (seed->GetParent()) PrintSeedClusters((AliITSUSeed*)seed->GetParent(), option);
1882 //
1883}
1884
c03e4f8a 1885//__________________________________________________________________
08419930 1886void AliITSUTrackerGlo::FlagSeedClusters(const AliITSUSeed* seed, Bool_t flg, UShort_t ref)
c03e4f8a 1887{
1888 // mark used clusters
1889 int lrID,clID;
1890 while (seed) {
08419930 1891 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1892 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1893 if (ref) { // do we need to set or delete cluster->track ref?
1894 if (flg) {
1895 if (!cl->GetRecoInfo()) cl->SetRecoInfo(ref); // set ref only if cluster already does not refer to other track, inc.counter
1896 }
1897 else {
1898 if (cl->GetRecoInfo()==ref) cl->SetRecoInfo(0); // unset reference only if it refers to ref, decrease counter
1899 }
1900 }
1901 }
c03e4f8a 1902 seed = (AliITSUSeed*)seed->GetParent();
1903 }
1904 //
1905}
1906
d99994b2 1907//__________________________________________________________________
1908void AliITSUTrackerGlo::CheckClusterUsage()
1909{
1910 // check cluster usage info
eb18cdc6 1911 printf("ClusterUsage at pass %d\n",fCurrPassID);
d99994b2 1912 for (int ilr=0;ilr<fNLrActive;ilr++) {
1913 AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
1914 int ncl = lr->GetNClusters();
1915 int nclUsed = 0;
1916 int nclUs0CntShr1 = 0;
1917 int nclUs1Cnt0 = 0;
1918 int nclUs1Shr1 = 0;
1919 int nclUseL = 0;
1920 for (int icl=0;icl<ncl;icl++) {
1921 AliITSUClusterPix *cl = (AliITSUClusterPix*)lr->GetCluster(icl);
1922 int usc = cl->GetClUsage();
1923 if (usc<1) {
1924 if (cl->IsClusterUsed() || cl->IsClusterShared()) {
1925 nclUs0CntShr1++;
1926 printf("Lr%d Cl%4d: UseCnt=0 but UseFlg=%d UseShr=%d\n",ilr,icl,cl->IsClusterUsed(),cl->IsClusterShared());
1927 }
1928 continue;
1929 }
1930 nclUsed++;
1931 if (!cl->IsClusterUsed()) {
1932 nclUs1Cnt0++;
1933 printf("Lr%d Cl%4d: UseCnt=%d but UseFlg=%d UseShr=%d\n",ilr,icl,usc, cl->IsClusterUsed(),cl->IsClusterShared());
1934 }
1935 if (cl->IsClusterShared()) {
1936 nclUs1Shr1++;
1937 printf("Lr%d Cl%4d: UseCnt=%d but UseFlg=%d UseShr=%d\n",ilr,icl,usc, cl->IsClusterUsed(),cl->IsClusterShared());
1938 }
1939 if (usc>1) {
1940 nclUseL++;
1941 printf("Lr%d Cl%4d: UseCnt=%d!! but UseFlg=%d UseShr=%d\n",ilr,icl,usc, cl->IsClusterUsed(),cl->IsClusterShared());
1942 }
1943 }
1944 printf("ClUsage Lr%d: Ncl=%5d Nusd=%5d (%5.2f) Use0CS1:%4d Use1C0:%4d Use1S1:%4d UseL:%d\n",
1945 ilr,ncl,nclUsed,ncl? float(nclUsed)/ncl : 0, nclUs0CntShr1,nclUs1Cnt0,nclUs1Shr1,nclUseL);
1946 }
1947}
1948
1949
08419930 1950#ifdef _ITSU_TUNING_MODE_
53870004 1951//__________________________________________________________________
781c79ac 1952TObjArray* AliITSUTrackerGlo::BookControlHistos(const char* pref)
53870004 1953{
1954 // book special control histos
08419930 1955 TString prefS = pref;
1956 prefS.ToLower();
781c79ac 1957 TObjArray* dest = new TObjArray;
08419930 1958 dest->SetOwner(kTRUE);
1959 //
bdb92766 1960 const int kNResDef=7;
1961 const double kResDef[kNResDef]={0.05,0.05,0.3, 0.05,1,0.5,1.5};
bdb92766 1962 const double ptMax=10;
1963 const double plMax=10;
1964 const double chiMax=100;
1965 const int nptbins=50;
1966 const int nresbins=400;
1967 const int nplbins=50;
1968 const int nchbins=200;
1969 const int maxBr = 15;
1970 const int maxCand = 200;
1971 TString ttl;
781c79ac 1972 for (int phase=0;phase<kNTrackingPhases;phase++) {
1973 for (int pass=0;pass<AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();pass++) {
bdb92766 1974 //
781c79ac 1975 for (int ilr=0;ilr<fNLrActive;ilr++) {
1976 //
1977 // ----------------- These are histos to be filled in Cluster2Tracks of each pass.
1978 // PropagateBack and RefitInward will be stored among the histos of 1st pass
1979 if (pass>0 && phase!=kClus2Tracks) continue;
1980 //
1981 double mxdf = ilr>=kNResDef ? kResDef[kNResDef-1] : kResDef[ilr];
1982 ttl = Form("Pass%d_S%d_residY%d_%s",pass,phase,ilr,pref);
1983 TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1984 dest->AddAtAndExpand(hdy,GetHistoID(ilr,kHResY,pass,phase));
1985 hdy->SetDirectory(0);
1986 //
1987 ttl = Form("Pass%d_S%d_residYPull%d_%s",pass,phase,ilr,pref);
1988 TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1989 dest->AddAtAndExpand(hdyp,GetHistoID(ilr,kHResYP,pass,phase));
1990 hdyp->SetDirectory(0);
1991 //
1992 ttl = Form("Pass%d_S%d_residZ%d_%s",pass,phase,ilr,pref);
1993 TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1994 dest->AddAtAndExpand(hdz,GetHistoID(ilr,kHResZ,pass,phase));
1995 hdz->SetDirectory(0);
1996 //
1997 ttl = Form("Pass%d_S%d_residZPull%d_%s",pass,phase,ilr,pref);
1998 TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1999 hdzp->SetDirectory(0);
2000 dest->AddAtAndExpand(hdzp,GetHistoID(ilr,kHResZP,pass,phase));
2001 //
2002 ttl = Form("Pass%d_S%d_chi2Cl%d_%s",pass,phase,ilr,pref);
2003 TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
2004 hchi->SetDirectory(0);
2005 dest->AddAtAndExpand(hchi,GetHistoID(ilr,kHChi2Cl,pass,phase));
2006 //
2007 // ------------------- These histos are filled for Clusters2Tracks only
2008 if (phase!=kClus2Tracks) continue;
2009 //
2010 ttl = Form("Pass%d_S%d_chi2Nrm%d_%s",pass,phase,ilr,pref);
2011 TH2F* hchiN = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
2012 hchiN->SetDirectory(0);
2013 dest->AddAtAndExpand(hchiN,GetHistoID(ilr,kHChi2Nrm,pass,phase));
2014 //
2015 ttl = Form("Pass%d_S%d_bestInBranch%d_%s",pass,phase,ilr,pref);
bdb92766 2016 TH2* hnbr = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxBr,-0.5,maxBr-0.5);
2017 hnbr->SetDirectory(0);
781c79ac 2018 dest->AddAtAndExpand(hnbr,GetHistoID(ilr,kHBestInBranch,pass,phase));
c03e4f8a 2019 //
781c79ac 2020 ttl = Form("Pass%d_S%d_bestInCands%d_%s",pass,phase,ilr,pref);
bdb92766 2021 TH2* hncn = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxCand,-0.5,maxCand-0.5);
2022 hncn->SetDirectory(0);
781c79ac 2023 dest->AddAtAndExpand(hncn,GetHistoID(ilr,kHBestInCand,pass,phase));
c03e4f8a 2024 //
781c79ac 2025 } // loop over layers
2026 //
2027 // these are special histos, filled not per layer but in the end of track fit in Clusters2Tracks in EVERY pass
2028 //
2029 if (phase==kClus2Tracks) {
2030 TH2* hchiMatch = 0;
2031 ttl = Form("Pass%d_Chi2Match_%s",pass,pref);
2032 hchiMatch = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
2033 hchiMatch->SetDirectory(0);
2034 dest->AddAtAndExpand(hchiMatch,GetHistoID(-1,kHChiMatch,pass,phase));
2035 //
2036 TH2* hchiSA = 0;
2037 ttl = Form("Pass%d_Chi2ITSSA_%s",pass,pref);
2038 hchiSA = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins/2,0.,chiMax/2);
2039 hchiSA->SetDirectory(0);
2040 dest->AddAtAndExpand(hchiSA,GetHistoID(-1,kHChiITSSA,pass,phase));
2041 //
53870004 2042 }
781c79ac 2043 } // loop over tracking passes
2044 }// loop over tracking phases
2045 //
2046 return dest;
53870004 2047}
781c79ac 2048
2049//__________________________________________________________________
2050Int_t AliITSUTrackerGlo::GetHistoID(Int_t lr, Int_t hid, Int_t pass, Int_t phase)
2051{
2052 // get id for the requested histo
2053 if (lr<0) lr=-1;
2054 lr++;
2055 return pass*kHistosPass + phase*kHistosPhase + lr*kMaxHID + hid;
2056 //
2057}
2058
53870004 2059#endif