]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITSUTrackerGlo.cxx
Add extra histos for neutral energy fraction vs centrality (from Chris)
[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
42c3d4bd 155//_________________________________________________________________________
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
32d38de2 175//_________________________________________________________________________
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_
70cb7fe4 266 AliDebug(1,Form("Processing track %d(esd%d) | M=%.3f Pt=%.3f | MCLabel: %d",itr,trID,fCurrESDtrack->GetMass(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 //
32d38de2 683 fCurrMass = esdTr->GetMass();
684 if (fCurrMass<kPionMass*0.9) fCurrMass = kPionMass; // don't trust to mu, e identification from TPCin
c61e50c3 685 //
6cd80116 686 hyp = new AliITSUTrackHyp(fNLrActive);
687 hyp->SetESDTrack(esdTr);
688 hyp->SetUniqueID(esdID);
689 hyp->SetMass(fCurrMass);
9cdcba2c 690 hyp->SetTPCSeed( new AliExternalTrackParam(*esdTr) );
6cd80116 691 SetTrackHyp(hyp,esdID);
e1ef49ad 692 Bool_t lut = fUseMatLUT; // for propagation from TPC use TGeo material budget always
693 fUseMatLUT = kFALSE;
70cb7fe4 694 if (!TransportToLayer(hyp->GetTPCSeed(),fITS->GetNLayers(), fITS->GetLrIDActive(fNLrActive-1), fITS->GetRITSTPCRef())) hyp->SetSkip(); // propagate to outer R of ITS
e1ef49ad 695 fUseMatLUT = lut;
3e4e3c23 696 //
6cd80116 697 return hyp;
32d38de2 698 // TO DO
699}
700
701//_________________________________________________________________________
70cb7fe4 702Bool_t AliITSUTrackerGlo::TransportToLayer(AliITSUSeed* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
32d38de2 703{
70cb7fe4 704 // transport seed from layerFrom to the entrance (along the direction of the track) of layerTo or to rLim (if>0), wathever is closer
32d38de2 705 //
8b16dbae 706 //
707 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
708 //
32d38de2 709 int dir = lTo > lFrom ? 1:-1;
710 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
711 Bool_t checkFirst = kTRUE;
70cb7fe4 712 Bool_t limReached = kFALSE;
32d38de2 713 while(lFrom!=lTo) {
32d38de2 714 if (lrFr) {
70cb7fe4 715 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
e2d2481c 716 checkFirst = kFALSE;
32d38de2 717 }
718 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
719 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
720 //
721 // go the entrance of the layer, assuming no materials in between
8b16dbae 722 double xToGo = lrTo->GetR(-dir);
70cb7fe4 723 if (rLim>0) {
724 if (dir>0) {
725 if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
726 }
727 else {
728 if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
729 }
730 }
0ddbf657 731 // double xts = xToGo;
b2d3bc3a 732 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
0ddbf657 733 //printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
734 //seed->Print("etp");
70cb7fe4 735 //
b2d3bc3a 736 return kFALSE;
737 }
e1ef49ad 738#ifdef _ITSU_DEBUG_
b515ad7e 739 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
e1ef49ad 740#endif
b2d3bc3a 741 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
0ddbf657 742 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
743 //seed->Print("etp");
b2d3bc3a 744 return kFALSE;
745 }
32d38de2 746 lrFr = lrTo;
70cb7fe4 747 if (limReached) break;
32d38de2 748 }
749 return kTRUE;
750 //
751}
752
3e4e3c23 753//_________________________________________________________________________
70cb7fe4 754Bool_t AliITSUTrackerGlo::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
3e4e3c23 755{
70cb7fe4 756 // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer
3e4e3c23 757 //
8b16dbae 758 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
759 //
3e4e3c23 760 int dir = lTo > lFrom ? 1:-1;
761 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
762 Bool_t checkFirst = kTRUE;
70cb7fe4 763 Bool_t limReached = kFALSE;
3e4e3c23 764 while(lFrom!=lTo) {
3e4e3c23 765 if (lrFr) {
e2d2481c 766 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
767 checkFirst = kFALSE;
3e4e3c23 768 }
769 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
770 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
771 //
0ddbf657 772 // go the entrance of the layer, assuming no materials in between
773 double xToGo = lrTo->GetR(-dir);
70cb7fe4 774 if (rLim>0) {
775 if (dir>0) {
776 if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
777 }
778 else {
779 if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
780 }
781 }
0ddbf657 782 // double xts = xToGo;
783 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
784 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
785 // seed->Print("etp");
786 return kFALSE;
787 }
e1ef49ad 788#ifdef _ITSU_DEBUG_
b515ad7e 789 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
e1ef49ad 790#endif
b515ad7e 791 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
792 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
793 //seed->Print("etp");
794 return kFALSE;
795 }
796 lrFr = lrTo;
70cb7fe4 797 if (limReached) break;
b515ad7e 798 }
799 return kTRUE;
800 //
801}
802
803//_________________________________________________________________________
804Bool_t AliITSUTrackerGlo::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
805{
806 // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
807 //
808 if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
809 //
810 int dir = lTo > lFrom ? 1:-1;
811 AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
812 Bool_t checkFirst = kTRUE;
813 while(lFrom!=lTo) {
814 if (lrFr) {
815 if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
816 checkFirst = kFALSE;
817 }
818 AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) );
819 if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
820 //
821 // go the entrance of the layer, assuming no materials in between
822 double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
823 //
824 // double xts = xToGo;
825 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
826 // printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
827 // seed->Print("etp");
828 return kFALSE;
829 }
830 if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
831 //
e1ef49ad 832#ifdef _ITSU_DEBUG_
b515ad7e 833 AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
e1ef49ad 834#endif
0ddbf657 835 if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
836 //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
837 //seed->Print("etp");
838 return kFALSE;
839 }
3e4e3c23 840 lrFr = lrTo;
841 }
842 return kTRUE;
843 //
844}
845
e2d2481c 846//_________________________________________________________________________
847Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
848{
849 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
850 // If check is requested, do this only provided the track has not exited the layer already
851 double xToGo = lr->GetR(dir);
852 if (check) { // do we need to track till the surface of the current layer ?
853 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
e1ef49ad 854#ifdef _ITSU_DEBUG_
15b02d69 855 AliDebug(3,Form(" dir:%d Cur: %e Tgt: %e",dir,Sqrt(curR2),xToGo));
e1ef49ad 856#endif
0ddbf657 857 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
858 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 859 }
860 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
861 // go via layer to its boundary, applying material correction.
862 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 863 //
e2d2481c 864 return kTRUE;
865 //
866}
867
868//_________________________________________________________________________
869Bool_t AliITSUTrackerGlo::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
870{
871 // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
872 // If check is requested, do this only provided the track has not exited the layer already
873 double xToGo = lr->GetR(dir);
874 if (check) { // do we need to track till the surface of the current layer ?
875 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 876 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
877 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
e2d2481c 878 }
e1ef49ad 879#ifdef _ITSU_DEBUG_
b515ad7e 880 AliDebug(2,Form(" dir=%d : from R=%.4f -> R=%.4f",dir,Sqrt(seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY()), xToGo));
e1ef49ad 881#endif
e2d2481c 882 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
883 // go via layer to its boundary, applying material correction.
884 if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
e7d83d38 885 //
e2d2481c 886 return kTRUE;
887 //
888}
889
890
891//_________________________________________________________________________
892Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliITSUSeed* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
893{
894 // go to the entrance of lr in direction dir, w/o applying material corrections.
895 // If check is requested, do this only provided the track did not reach the layer already
896 double xToGo = lr->GetR(-dir);
897 if (check) { // do we need to track till the surface of the current layer ?
898 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 899 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
900 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 901 }
902 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
903 // go via layer to its boundary, applying material correction.
904 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
905 return kTRUE;
906 //
907}
908
909//_________________________________________________________________________
910Bool_t AliITSUTrackerGlo::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
911{
912 // go to the entrance of lr in direction dir, w/o applying material corrections.
913 // If check is requested, do this only provided the track did not reach the layer already
914 double xToGo = lr->GetR(-dir);
915 if (check) { // do we need to track till the surface of the current layer ?
916 double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
0ddbf657 917 if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
918 else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
e2d2481c 919 }
920 if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
921 // go via layer to its boundary, applying material correction.
922 if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
923 return kTRUE;
924 //
925}
926
32d38de2 927//_________________________________________________________________________
c61e50c3 928Bool_t AliITSUTrackerGlo::GetRoadWidth(AliITSUSeed* seed, int ilrA)
32d38de2 929{
930 // calculate road width in terms of phi and z for the track which MUST be on the external radius of the layer
931 // as well as some aux info
932 double bz = GetBz();
933 AliITSURecoLayer* lrA = fITS->GetLayerActive(ilrA);
c61e50c3 934 seed->GetXYZ(&fTrImpData[kTrXIn]); // lab position at the entrance from above
53870004 935 static AliExternalTrackParam sc; // seed copy for manipulations
943e1898 936 sc = *seed;
32d38de2 937 //
c61e50c3 938 fTrImpData[kTrPhiIn] = ATan2(fTrImpData[kTrYIn],fTrImpData[kTrXIn]);
943e1898 939 if (!sc.Rotate(fTrImpData[kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
c61e50c3 940 double dr = lrA->GetDR(); // approximate X dist at the inner radius
943e1898 941 if (!sc.GetXYZAt(sc.GetX()-dr, bz, fTrImpData + kTrXOut)) {
32d38de2 942 // special case: track does not reach inner radius, might be tangential
943e1898 943 double r = sc.GetD(0,0,bz);
32d38de2 944 double x;
bdb92766 945 if (!sc.GetXatLabR(r,x,bz,-1)) return kFALSE;
943e1898 946 dr = Abs(sc.GetX() - x);
bdb92766 947 if (!sc.GetXYZAt(x, bz, fTrImpData + kTrXOut)) return kFALSE;
32d38de2 948 }
949 //
c61e50c3 950 fTrImpData[kTrPhiOut] = ATan2(fTrImpData[kTrYOut],fTrImpData[kTrXOut]);
943e1898 951 double sgy = sc.GetSigmaY2() + dr*dr*sc.GetSigmaSnp2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(ilrA);
952 double sgz = sc.GetSigmaZ2() + dr*dr*sc.GetSigmaTgl2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(ilrA);
ee54014a 953 sgy = Sqrt(sgy)*fCurrTrackCond->GetNSigmaRoadY(ilrA);
954 sgz = Sqrt(sgz)*fCurrTrackCond->GetNSigmaRoadZ(ilrA);
5785a9d8 955 double dphi0 = 0.5*Abs(fTrImpData[kTrPhiOut]-fTrImpData[kTrPhiIn]);
956 double phi0 = 0.5*(fTrImpData[kTrPhiOut]+fTrImpData[kTrPhiIn]);
957 if ( dphi0>(0.5*Pi()) ) {
958 // special case of angles around pi
959 dphi0 = Abs(phi0);
960 phi0 += Pi();
961 }
962
963 fTrImpData[kTrPhi0] = phi0;
173b3073 964 fTrImpData[kTrZ0] = 0.5*(fTrImpData[kTrZOut]+fTrImpData[kTrZIn]);
31e2f5d4 965 dphi0 += sgy/lrA->GetR();
966 fTrImpData[kTrDPhi] = dphi0<PiOver2() ? dphi0 : PiOver2();
173b3073 967 fTrImpData[kTrDZ] = 0.5*Abs(fTrImpData[kTrZOut]-fTrImpData[kTrZIn]) + sgz;
32d38de2 968 //
969 return kTRUE;
970}
971
b8b59e05 972//________________________________________
973void AliITSUTrackerGlo::ResetSeedsPool()
974{
975 // mark all seeds in the pool as unused
0ddbf657 976 AliDebug(-1,Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
b8b59e05 977 fNFreeSeeds = 0;
978 fSeedsPool.Clear(); // seeds don't allocate memory
979}
980
981
982//________________________________________
983void AliITSUTrackerGlo::MarkSeedFree(AliITSUSeed *sd)
32d38de2 984{
b8b59e05 985 // account that this seed is "deleted"
986 int id = sd->GetPoolID();
987 if (id<0) {
988 AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
989 return;
990 }
991 // AliInfo(Form("%d %p",id, seed));
992 fSeedsPool.RemoveAt(id);
993 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
994 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
32d38de2 995}
f8832015 996
997//_________________________________________________________________________
998Int_t AliITSUTrackerGlo::CheckCluster(AliITSUSeed* track, Int_t lr, Int_t clID)
999{
1000 // Check if the cluster (in tracking frame!) is matching to track.
1001 // The track must be already propagated to sensor tracking frame.
1002 // Returns: kStopSearchOnSensor if the search on given sensor should be stopped,
1003 // kClusterMatching if the cluster is matching
1004 // kClusterMatching otherwise
1005 //
f8832015 1006 const double kTolerX = 5e-4;
c03e4f8a 1007 // AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1008 AliCluster *cl = fCurrLayer->GetCluster(clID);
f8832015 1009 //
c8d1f258 1010 Bool_t goodCl = kFALSE;
0ddbf657 1011 int currLabel = Abs(fCurrESDtrMClb);
716ccba7 1012 //
8b16dbae 1013 if (cl->GetLabel(0)>=0) {
1014 for (int i=0;i<3;i++) if (cl->GetLabel(i)>=0 && cl->GetLabel(i)==currLabel) {goodCl = kTRUE; break;}
1015 }
c8d1f258 1016 //
f8832015 1017 if (TMath::Abs(cl->GetX())>kTolerX) { // if due to the misalingment X is large, propagate track only
c8d1f258 1018 if (!track->PropagateParamOnlyTo(track->GetX()+cl->GetX(),GetBz())) {
08419930 1019 //
e1ef49ad 1020#ifdef _ITSU_DEBUG_
08419930 1021 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1022 AliDebug(2,Form("Lost good cl on L:%d failed propagation. |ESDtrack#%d (MClb:%d)",lr,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1023 track->Print("etp");
15b02d69 1024 cl->Print();
1025 }
e1ef49ad 1026#endif
c8d1f258 1027 return kStopSearchOnSensor; // propagation failed, seedT is intact
1028 }
f8832015 1029 }
3dd9c283 1030 double dy = cl->GetY()-track->GetY();
c8d1f258 1031 double dz = cl->GetZ()-track->GetZ();
0091e9f0 1032 //
f8832015 1033 double dy2 = dy*dy;
1034 double tol2 = (track->GetSigmaY2() + AliITSUReconstructor::GetRecoParam()->GetSigmaY2(lr))*
ee54014a 1035 fCurrTrackCond->GetNSigmaRoadY(lr)*fCurrTrackCond->GetNSigmaRoadY(lr); // RS TOOPTIMIZE
f8832015 1036 if (dy2>tol2) { // the clusters are sorted in Z(col) then in Y(row).
08419930 1037 //
e1ef49ad 1038#ifdef _ITSU_DEBUG_
08419930 1039 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1040 AliDebug(2,Form("Lost good cl: dy2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",dy2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1041 track->Print("etp");
15b02d69 1042 cl->Print();
53870004 1043 PrintSeedClusters(track);
15b02d69 1044 }
e1ef49ad 1045#endif
cb50e082 1046 // clusters are sorted in increasing Y and the loop goes from last (highers Y) to 1st.
1047 // if the track has too large y for this cluster (dy<0) it will be so for all other clusters of the sensor
1048 if (dy<0) return kStopSearchOnSensor; // No chance that other cluster of this sensor will match (all Y's will be even larger)
f8832015 1049 else return kClusterNotMatching; // Other clusters may match
1050 }
c8d1f258 1051 double dz2 = dz*dz;
f8832015 1052 tol2 = (track->GetSigmaZ2() + AliITSUReconstructor::GetRecoParam()->GetSigmaZ2(lr))*
ee54014a 1053 fCurrTrackCond->GetNSigmaRoadZ(lr)*fCurrTrackCond->GetNSigmaRoadZ(lr); // RS TOOPTIMIZE
c8d1f258 1054 if (dz2>tol2) {
e1ef49ad 1055#ifdef _ITSU_DEBUG_
08419930 1056 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1057 AliDebug(2,Form("Lost good cl on L:%d : dz2=%e > tol2=%e |ESDtrack#%d (MClb:%d)",lr,dz2,tol2,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1058 track->Print("etp");
15b02d69 1059 cl->Print();
53870004 1060 PrintSeedClusters(track);
15b02d69 1061 }
e1ef49ad 1062#endif
c8d1f258 1063 return kClusterNotMatching; // Other clusters may match
1064 }
f8832015 1065 //
1066 // check chi2
1067 Double_t p[2]={cl->GetY(), cl->GetZ()};
1068 Double_t cov[3]={cl->GetSigmaY2(), cl->GetSigmaYZ(), cl->GetSigmaZ2()};
1069 double chi2 = track->GetPredictedChi2(p,cov);
76390254 1070 //
ee54014a 1071 if (chi2>fCurrTrackCond->GetMaxTr2ClChi2(lr)) {
e1ef49ad 1072#ifdef _ITSU_DEBUG_
08419930 1073 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1074 AliDebug(2,Form("Lost good cl on L:%d : Chi2=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
ee54014a 1075 lr,chi2,fCurrTrackCond->GetMaxTr2ClChi2(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb));
3dd9c283 1076 track->Print("etp");
1077 cl->Print("");
53870004 1078 PrintSeedClusters(track);
3dd9c283 1079 }
e1ef49ad 1080#endif
c8d1f258 1081 return kClusterNotMatching;
1082 }
f8832015 1083 //
1084 track = NewSeedFromPool(track); // input track will be reused, use its clone for updates
943e1898 1085 if (!track->Update()) {
e1ef49ad 1086#ifdef _ITSU_DEBUG_
08419930 1087 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
53870004 1088 AliDebug(2,Form("Lost good cluster on L:%d : failed to update",lr));
38997c3c 1089 track->Print("etp");
1090 cl->Print("");
53870004 1091 PrintSeedClusters(track);
38997c3c 1092 }
e1ef49ad 1093#endif
b8b59e05 1094 MarkSeedFree(track);
c8d1f258 1095 return kClusterNotMatching;
1096 }
f8832015 1097 track->SetChi2Cl(chi2);
1098 track->SetLrClusterID(lr,clID);
08419930 1099 //
1100 if (track->GetChi2GloNrm()>fCurrTrackCond->GetMaxChi2GloNrm(lr)) {
e1ef49ad 1101#ifdef _ITSU_DEBUG_
08419930 1102 if (AliDebugLevelClass()>2 && goodCl && !track->ContainsFake()) {
1103 AliDebug(2,Form("Lost good cl on L:%d : Chi2Glo=%e > Chi2Max=%e |dy: %+.3e dz: %+.3e |ESDtrack#%d (MClb:%d)\n",
1104 lr,track->GetChi2GloNrm(),fCurrTrackCond->GetMaxChi2GloNrm(lr),dy,dz,fCurrESDtrack->GetID(),fCurrESDtrMClb));
1105 track->Print("etp");
1106 cl->Print("");
1107 PrintSeedClusters(track);
1108 }
e1ef49ad 1109#endif
08419930 1110 MarkSeedFree(track);
1111 return kClusterNotMatching;
1112 }
c03e4f8a 1113 // cl->IncreaseClusterUsage(); // do this only for winners
f8832015 1114 //
716ccba7 1115 track->SetFake(!goodCl);
1116 //
e1ef49ad 1117#ifdef _ITSU_DEBUG_
0ddbf657 1118 AliDebug(2,Form("AddCl(%d) Cl%d lr:%d: dY:%+8.4f dZ:%+8.4f (MC: %5d %5d %5d) |Chi2=%f(%c)| ",
716ccba7 1119 goodCl,clID,lr,dy,dz2,cl->GetLabel(0),cl->GetLabel(1),cl->GetLabel(2), chi2, track->IsFake() ? '-':'+'));
e1ef49ad 1120#endif
c8d1f258 1121 //
38997c3c 1122 AddSeedBranch(track);
f8832015 1123 //
1124 return kClusterMatching;
1125}
c8d1f258 1126
1127//_________________________________________________________________________
1128Bool_t AliITSUTrackerGlo::NeedToKill(AliITSUSeed *seed, Int_t flag)
1129{
1130 // check if the seed should not be discarded
1131 const UShort_t kMask = 0xffff;
1132 if (flag==kMissingCluster) {
1133 int lastChecked = seed->GetLayerID();
1134 UShort_t patt = seed->GetHitsPattern();
097e4cdf 1135 if (lastChecked) patt |= (~(kMask<<lastChecked))&fCurrTrackCond->GetAllowedLayers(); // not all layers were checked, complete unchecked ones by potential hits
42c3d4bd 1136 Bool_t seedOK = fCurrTrackCond->CheckPattern(patt);
c8d1f258 1137 return !seedOK;
1138 }
1139 return kTRUE;
1140}
44785f3e 1141
1142//______________________________________________________________________________
1143Bool_t AliITSUTrackerGlo::PropagateSeed(AliITSUSeed *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
1144{
1145 // propagate seed to given x applying material correction if requested
1146 const Double_t kEpsilon = 1e-5;
1147 Double_t xpos = seed->GetX();
1148 Int_t dir = (xpos<xToGo) ? 1:-1;
e1ef49ad 1149 Double_t xyz0[3],xyz1[3];
44785f3e 1150 //
8b16dbae 1151 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
f497470c 1152 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
44785f3e 1153 while ( (xToGo-xpos)*dir > kEpsilon){
1154 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1155 Double_t x = xpos+step;
1156 Double_t bz=GetBz(); // getting the local Bz
0ddbf657 1157 if (!seed->PropagateToX(x,bz)) return kFALSE;
8b16dbae 1158 double ds = 0;
1159 if (matCorr || updTime) {
44785f3e 1160 xyz0[0]=xyz1[0]; // global pos at the beginning of step
1161 xyz0[1]=xyz1[1];
1162 xyz0[2]=xyz1[2];
1163 seed->GetXYZ(xyz1); // // global pos at the end of step
8b16dbae 1164 if (matCorr) {
e1ef49ad 1165 Double_t xrho,xx0;
1166 ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
8b16dbae 1167 if (dir>0) xrho = -xrho; // outward should be negative
b2d3bc3a 1168 if (!seed->ApplyMaterialCorrection(xx0,xrho,mass,kFALSE)) {AliInfo("Fail2"); seed->Print("etp"); return kFALSE;}
8b16dbae 1169 }
1170 else { // matCorr is not requested but time integral is
1171 double d0 = xyz1[0]-xyz0[0];
1172 double d1 = xyz1[1]-xyz0[1];
1173 double d2 = xyz1[2]-xyz0[2];
1174 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1175 }
44785f3e 1176 }
8b16dbae 1177 if (updTime) seed->AddTimeStep(ds);
44785f3e 1178 xpos = seed->GetX();
1179 }
1180 return kTRUE;
1181}
716ccba7 1182
3e4e3c23 1183//______________________________________________________________________________
1184Bool_t AliITSUTrackerGlo::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
1185{
1186 // propagate seed to given x applying material correction if requested
1187 const Double_t kEpsilon = 1e-5;
1188 Double_t xpos = seed->GetX();
1189 Int_t dir = (xpos<xToGo) ? 1:-1;
e1ef49ad 1190 Double_t xyz0[3],xyz1[3];
3e4e3c23 1191 //
e1ef49ad 1192#ifdef _ITSU_DEBUG_
b515ad7e 1193 if (AliDebugLevelClass()>1) {
1194 AliDebug(2,Form("Before Propagate to X=%f with M=%.3f MaxStep=%.4f MatCorr=%d",xToGo,mass,maxStep,matCorr));
1195 seed->AliExternalTrackParam::Print();
1196 }
e1ef49ad 1197#endif
3e4e3c23 1198 Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
1199 if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position
1200 while ( (xToGo-xpos)*dir > kEpsilon){
1201 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1202 Double_t x = xpos+step;
1203 Double_t bz=GetBz(); // getting the local Bz
1204 if (!seed->PropagateTo(x,bz)) return kFALSE;
1205 double ds = 0;
1206 if (matCorr || updTime) {
1207 xyz0[0]=xyz1[0]; // global pos at the beginning of step
1208 xyz0[1]=xyz1[1];
1209 xyz0[2]=xyz1[2];
1210 seed->GetXYZ(xyz1); // // global pos at the end of step
1211 //
1212 if (matCorr) {
e1ef49ad 1213 Double_t xrho,xx0;
1214 ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
3e4e3c23 1215 if (dir>0) xrho = -xrho; // outward should be negative
1216 if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
3e4e3c23 1217 }
1218 else { // matCorr is not requested but time integral is
1219 double d0 = xyz1[0]-xyz0[0];
1220 double d1 = xyz1[1]-xyz0[1];
1221 double d2 = xyz1[2]-xyz0[2];
1222 ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
1223 }
1224 }
1225 if (updTime) seed->AddTimeStep(ds);
1226 //
1227 xpos = seed->GetX();
1228 }
e1ef49ad 1229#ifdef _ITSU_DEBUG_
b515ad7e 1230 if (AliDebugLevelClass()>1) {
1231 AliDebug(2,Form("After Propagate to X=%f",xToGo));
1232 seed->AliExternalTrackParam::Print();
1233 }
e1ef49ad 1234#endif
3e4e3c23 1235 return kTRUE;
1236}
1237
9cdcba2c 1238//______________________________________________________________________________
1239Bool_t AliITSUTrackerGlo::FinalizeHypothesis(AliITSUTrackHyp* hyp)
1240{
1241 // finalize single TPC track hypothesis
1242 if (!hyp || hyp->GetESDTrack()->IsOn(AliESDtrack::kITSin)) return kFALSE;
1243 AliITSUSeed* winner = 0;
1244 fCurrHyp = hyp;
70cb7fe4 1245 fCurrMass = hyp->GetMass();
d99994b2 1246 if (!(winner=hyp->DefineWinner(fCurrTrackCond->GetActiveLrInner()))) return kFALSE;
9cdcba2c 1247 CookMCLabel(hyp);
1248 //
08419930 1249#ifdef _ITSU_TUNING_MODE_ // fill tuning histos
1250 TObjArray* dest = hyp->GetLabel()>=0 ? fCHistoArrCorr : fCHistoArrFake;
1251 if (dest) {
1252 AliITSUSeed* sd = winner;
1253 double htrPt = hyp->Pt();
1254 do {
1255 int clID,lrID;
1256 if ( (clID=sd->GetLrCluster(lrID))<0 ) continue;
781c79ac 1257 ((TH2*)dest->At( GetHistoID(lrID,kHChi2Nrm ,fCurrPassID) ))->Fill(htrPt,sd->GetChi2GloNrm());
1258 ((TH2*)dest->At( GetHistoID(lrID,kHBestInBranch,fCurrPassID) ))->Fill(htrPt,sd->GetOrdBranch());
1259 ((TH2*)dest->At( GetHistoID(lrID,kHBestInCand ,fCurrPassID) ))->Fill(htrPt,sd->GetOrdCand());
08419930 1260 //
1261 if (dest==fCHistoArrFake && !sd->IsFake()) continue; // for the fake seeds fill only fake clusters part
1262 //
781c79ac 1263 ((TH2*)dest->At( GetHistoID(lrID,kHResY ,fCurrPassID) ))->Fill(htrPt,sd->GetResidY());
1264 ((TH2*)dest->At( GetHistoID(lrID,kHResZ ,fCurrPassID) ))->Fill(htrPt,sd->GetResidZ());
1265 ((TH2*)dest->At( GetHistoID(lrID,kHResYP ,fCurrPassID) ))->Fill(htrPt,sd->GetPullY());
1266 ((TH2*)dest->At( GetHistoID(lrID,kHResZP ,fCurrPassID) ))->Fill(htrPt,sd->GetPullZ());
1267 ((TH2*)dest->At( GetHistoID(lrID,kHChi2Cl ,fCurrPassID) ))->Fill(htrPt,sd->GetChi2Cl());
08419930 1268 //
1269 } while((sd=(AliITSUSeed*)sd->GetParent()));
1270 //
781c79ac 1271 ((TH2*)dest->At( GetHistoID(-1,kHChiMatch,fCurrPassID) ))->Fill(htrPt,winner->GetChi2ITSTPC());
1272 ((TH2*)dest->At( GetHistoID(-1,kHChiITSSA,fCurrPassID) ))->Fill(htrPt,winner->GetChi2ITSSA());
9cdcba2c 1273 }
9cdcba2c 1274 //
08419930 1275#endif
1276 //
1277 CheckClusterSharingConflicts(hyp);
1278 return hyp->GetWinner(); // winner might change of disappear after resolving conflicts
1279 //
9cdcba2c 1280}
1281
3e4e3c23 1282//______________________________________________________________________________
1283void AliITSUTrackerGlo::FinalizeHypotheses()
1284{
1285 // select winner for each hypothesis, remove cl. sharing conflicts
9cdcba2c 1286 AliCodeTimerAuto("",0);
3e4e3c23 1287 //
d99994b2 1288 int* esdTrackIndex = fESDIndex.GetArray();
1289 for (int ih=0;ih<fNTracksESD;ih++) FinalizeHypothesis(GetTrackHyp(esdTrackIndex[ih])); // finlize in decreasing pt order
c03e4f8a 1290 //
d99994b2 1291 for (int ih=fNTracksESD;ih--;) UpdateESDTrack(GetTrackHyp(esdTrackIndex[ih]),AliESDtrack::kITSin);
08419930 1292 //
1293}
bdb92766 1294
08419930 1295//______________________________________________________________________________
1296void AliITSUTrackerGlo::CheckClusterSharingConflicts(AliITSUTrackHyp* hyp)
1297{
1298 // remove eventual cluster sharing conflicts
1299 AliITSUSeed* winner = 0;
1300 if (!(winner=hyp->GetWinner())) return;
1301 UShort_t idH = (UShort_t)hyp->GetUniqueID()+1;
1302 int lrID,clID;
d99994b2 1303 int inLr = fCurrTrackCond->GetActiveLrInner();
08419930 1304 AliITSUSeed* winSD=winner;
1305 do {
1306 if ( (clID=winSD->GetLrCluster(lrID))<0 ) continue;
1307 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1308 Int_t refID = cl->GetRecoInfo(); // was it already referred by some track (with refID-1 if >0)
1309 if (!refID) {cl->SetRecoInfo(idH); continue;} // if not, refer from track to cluster
1310 if (refID==idH) continue; // ignore reference to itself
1311 // the cluster is already used by other track, need to resolve the conflict
1312 AliITSUTrackHyp* hypC = GetTrackHyp(refID-1); // competitor
1313 // printf("Ref to %d (%p) from %d Cl %d %d\n",refID-1,hypC, idH-1,clID,lrID);
1314 AliITSUSeed* winnerC = hypC->GetWinner();
1315 if (!winnerC) {
1316 printf("Missing winner, ERROR\n");
1317 continue;
1318 }
6091ad4e 1319 //
e1ef49ad 1320#ifdef _ITSU_DEBUG_
08419930 1321 if (AliDebugLevelClass()>1) { //------------------------ DEBUG -------------------
1322 AliInfo(Form("Shared cl#%4d on Lr:%d: Hyp%3d/Q:%6.3f vs Hyp%3d/Q:%6.3f",
1323 clID,lrID,idH-1,winner->GetQualityVar(),refID-1,winnerC->GetQualityVar()));
1324 // dump winner of hyp
1325 printf("#%4d Pt:%.3f Chi:%6.3f/%6.3f/%6.3f Ncl:%d MCits%+5d MCtpc:%+5d |",
1326 idH-1,winner->Pt(),winner->GetChi2GloNrm(),winner->GetChi2ITSTPC(),winner->GetChi2ITSSA(),
1327 winner->GetNLayersHit(),hyp->GetITSLabel(),hyp->GetESDTrack()->GetTPCLabel());
1328 int prevL=-1;
1329 AliITSUSeed* sd = winner;
1330 do {
1331 int lrs;
1332 int clIDt = sd->GetLrCluster(lrs);
1333 if (clIDt<0) continue;
1334 while( lrs>++prevL ) printf("%4s ","----");
1335 printf("%4d (%5.1f)",clIDt,sd->GetChi2Cl());
1336 } while ((sd=(AliITSUSeed*)sd->GetParent()));
1337 printf("|\n");
1338 // dump winner of hypC
1339 printf("#%4d Pt:%.3f Chi:%6.3f/%6.3f/%6.3f Ncl:%d MCits%+5d MCtpc:%+5d |",
1340 refID-1,winnerC->Pt(),winnerC->GetChi2GloNrm(),winnerC->GetChi2ITSTPC(),winnerC->GetChi2ITSSA(),
1341 winnerC->GetNLayersHit(),hypC->GetITSLabel(),hypC->GetESDTrack()->GetTPCLabel());
1342 prevL=-1;
1343 sd = winnerC;
6091ad4e 1344 do {
08419930 1345 int lrs;
1346 int clIDt = sd->GetLrCluster(lrs);
1347 if (clIDt<0) continue;
1348 while( lrs>++prevL ) printf("%4s ","----");
1349 printf("%4d (%5.1f)",clIDt,sd->GetChi2Cl());
1350 } while ((sd=(AliITSUSeed*)sd->GetParent()));
1351 printf("|\n");
1352 } //-------------------------------------------------- DEBUG -------------------
e1ef49ad 1353#endif
08419930 1354 //
1355 if (winner->GetQualityVar()<winnerC->GetQualityVar()) { // current tracks is better then competitor track
1356 FlagSeedClusters(winnerC,kFALSE,refID); // unflag cluster usage by loser
1357 cl->SetRecoInfo(idH);
1358 //
1359 //if ( hypC->GetLabel()>=0) AliInfo(Form("KILLING CORRECT: %d",refID-1));
1360 winnerC->Kill();
d99994b2 1361 if (hypC->DefineWinner(inLr)) { // find new winner instead of suppressed candidate
08419930 1362 CookMCLabel(hypC);
1363 CheckClusterSharingConflicts(hypC); // and check its sharing conflicts again
1364 if (winner->IsKilled()) break; // the current winner might have been killed during check of new winner of hypC!
1365 }
6091ad4e 1366 }
08419930 1367 else { // competitor hypC is better than the hyp
1368 FlagSeedClusters(winner,kFALSE,idH); // unflag cluster usage by loser
1369 winner->Kill();
1370 //if ( hyp->GetLabel()>=0) AliInfo(Form("KILLING CORRECT: %d",idH-1));
d99994b2 1371 if (hyp->DefineWinner(inLr)) {
08419930 1372 CookMCLabel(hyp);
1373 CheckClusterSharingConflicts(hyp);
9cdcba2c 1374 }
08419930 1375 break;
9cdcba2c 1376 }
08419930 1377 //
1378 } while ((winSD=(AliITSUSeed*)winSD->GetParent()));
1379 //
1380 return;
3e4e3c23 1381}
c51c3124 1382
1383//______________________________________________________________________________
68a0f687 1384void AliITSUTrackerGlo::UpdateESDTrack(AliITSUTrackHyp* hyp,Int_t flag)
c51c3124 1385{
1386 // update ESD track with current best hypothesis
08419930 1387 if (!hyp) return;
c51c3124 1388 AliESDtrack* esdTr = hyp->GetESDTrack();
d99994b2 1389 if (!esdTr || esdTr->IsOn(flag)) return;
c51c3124 1390 AliITSUSeed* win = hyp->GetWinner();
08419930 1391 if (!win || win->IsKilled()) return;
70cb7fe4 1392 double chiSav;
1393 //
68a0f687 1394 switch (flag) {
1395 case AliESDtrack::kITSin:
1396 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
08419930 1397 fCountITSin++;
dc443a34 1398 // set fakes cluster info
1399 {
1400 UShort_t clfake = 0;
c7e63626 1401 Int_t clSplit = 0;
dc443a34 1402 AliITSUSeed* sd = win;
c7e63626 1403 int ip=0;
dc443a34 1404 do {
c7e63626 1405 int lr, clID = sd->GetLrCluster(lr);
1406 if (sd->IsFake()) clfake |= 0x1<<lr;
1407 if (clID>=0) {
1408 esdTr->SetITSModuleIndex(ip++, sd->GetLrClusterID());
1409 AliITSUClusterPix *cl = (AliITSUClusterPix*)fITS->GetLayerActive(lr)->GetCluster(clID);
1410#ifdef _ITSU_TUNING_MODE_
1411 if (cl->IsSplit()) clSplit |= 0x1<<lr;
1412#endif
1413 }
dc443a34 1414 } while ((sd=(AliITSUSeed*)sd->GetParent()));
0e40ce66 1415 //
1416 // RS: Temporary set special flag for tracks from the afterburner
1417 if (fCurrPassID>0) clfake |= 0x1<<7;
1418 //
dc443a34 1419 esdTr->SetITSSharedMap(clfake);
c7e63626 1420 esdTr->SetITSModuleIndex(10,clSplit);
dc443a34 1421 }
097e4cdf 1422 // TEMPORARY: store iteration id
1423 esdTr->SetITSModuleIndex(11,fCurrPassID);
68a0f687 1424 break;
1425 //
1426 case AliESDtrack::kITSout:
70cb7fe4 1427 // here the stored chi2 will correspond to backward ITS-SA fit
68a0f687 1428 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
08419930 1429 fCountITSout++;
68a0f687 1430 // TODO: avoid setting friend
1431 break;
1432 //
1433 case AliESDtrack::kITSrefit:
70cb7fe4 1434
1435 // at the moment fill as a chi2 the TPC-ITS matching chi2
1436 chiSav = hyp->GetChi2();
1437 hyp->SetChi2(win->GetChi2ITSTPC());
68a0f687 1438 esdTr->UpdateTrackParams(hyp,flag); // update kinematics
70cb7fe4 1439 hyp->SetChi2(chiSav);
08419930 1440 fCountITSrefit++;
68a0f687 1441 // TODO: avoid setting cluster info
1442 break;
1443 default:
1444 AliFatal(Form("Unknown flag %d",flag));
1445 }
c51c3124 1446 //
bdb92766 1447 esdTr->SetITSLabel(hyp->GetITSLabel());
c51c3124 1448 // transfer module indices
1449 // TODO
1450}
1451
1452//______________________________________________________________________________
7ce7f445 1453Double_t AliITSUTrackerGlo::RefitTrack(AliITSUTrackHyp* trc, Double_t rDest, Int_t &nclFit, Int_t stopCond)
c51c3124 1454{
9cdcba2c 1455 // refit track till radius rDest. The cluster,mass info is taken from fCurrHyp (and its winner seed)
1456 // if stopCond<0 : propagate till last cluster then stop
1457 // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail
1458 // if stopCond>0 : rDest must be reached
8b16dbae 1459 //
c51c3124 1460 double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
1461 int dir,lrStart,lrStop;
1462 //
8b16dbae 1463 dir = rCurr<rDest ? 1 : -1;
1464 lrStart = fITS->FindFirstLayerID(rCurr,dir);
70cb7fe4 1465 lrStop = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
1466 //
e1ef49ad 1467#ifdef _ITSU_DEBUG_
15b02d69 1468 if (AliDebugLevelClass()>2) {
1469 printf("Refit %d: Lr:%d (%f) -> Lr:%d (%f)\n",dir,lrStart,rCurr, lrStop,rDest);
9cdcba2c 1470 printf("Before refit: "); trc->Print();
15b02d69 1471 }
e1ef49ad 1472#endif
8b16dbae 1473 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));
1474 //
9cdcba2c 1475 Int_t clInfo[2*AliITSUAux::kMaxLayers];
1476 Int_t nCl = fCurrHyp->FetchClusterInfo(clInfo);
1477 fCurrMass = fCurrHyp->GetMass();
7ce7f445 1478 AliITSUTrackHyp tmpTr(*(AliKalmanTrack*)trc);
9cdcba2c 1479 double chi2 = 0;
5e3587be 1480 int iclLr[2],nclLr;
1481 nclFit = 0;
c51c3124 1482 //
c7e63626 1483 int lrStop1 = lrStop+dir;
70cb7fe4 1484 for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
c51c3124 1485 AliITSURecoLayer* lr = fITS->GetLayer(ilr);
1486 if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
8b16dbae 1487 int ilrA2,ilrA = lr->GetActiveID();
1488 // passive layer or active w/o hits will be traversed on the way to next cluster
9cdcba2c 1489 if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) continue;
8b16dbae 1490 //
8b16dbae 1491 // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
1492 nclLr=0;
1493 if (dir>0) { // clusters are stored in increasing radius order
9cdcba2c 1494 iclLr[nclLr++]=clInfo[ilrA2++];
1495 if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2];
8b16dbae 1496 }
1497 else {
9cdcba2c 1498 if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1];
1499 iclLr[nclLr++]=clInfo[ilrA2];
8b16dbae 1500 }
1501 //
b515ad7e 1502 Bool_t transportedToLayer = kFALSE;
8b16dbae 1503 for (int icl=0;icl<nclLr;icl++) {
1504 AliITSUClusterPix* clus = (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
1505 AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
e7d83d38 1506 if (!tmpTr.Rotate(sens->GetPhiTF())) {
e1ef49ad 1507#ifdef _ITSU_DEBUG_
fb0eae48 1508 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1509 AliDebug(2,Form("Failed on rotate to %f | ESDtrack#%d (MClb:%d)",sens->GetPhiTF(),trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1510#endif
9cdcba2c 1511 return -1;
e7d83d38 1512 }
b515ad7e 1513 //
1514 double xClus = sens->GetXTF()+clus->GetX();
1515 if (!transportedToLayer) {
1516 if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
e1ef49ad 1517#ifdef _ITSU_DEBUG_
846c381e 1518 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1519 AliDebug(2,Form("Failed to transport %d -> %d | ESDtrack#%d (MClb:%d)\n",lrStart,ilr,trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1520#endif
9cdcba2c 1521 return -1; // go to the entrance to the layer
b515ad7e 1522 }
1523 lrStart = ilr;
1524 transportedToLayer = kTRUE;
1525 }
1526 //
e1ef49ad 1527#ifdef _ITSU_DEBUG_
b515ad7e 1528 if (AliDebugLevelClass()>1) {
1529 AliDebug(2,Form("Propagate to cl:%d on lr %d Need to go %.4f -> %.4f",icl,ilrA,tmpTr.GetX(),xClus));
b515ad7e 1530 }
e1ef49ad 1531#endif
b515ad7e 1532 //
1533 if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) {
e1ef49ad 1534#ifdef _ITSU_DEBUG_
846c381e 1535 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1536 AliDebug(2,Form("Failed on propagate to %f (dir=%d) | ESDtrack#%d (MClb:%d)",xClus,dir,trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1537#endif
9cdcba2c 1538 return -1;
e7d83d38 1539 }
76390254 1540 //
9cdcba2c 1541 Double_t p[2]={clus->GetY(), clus->GetZ()};
1542 Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()};
1543 double chi2cl = tmpTr.GetPredictedChi2(p,cov);
1544 chi2 += chi2cl;
1545 //
08419930 1546#ifdef _ITSU_TUNING_MODE_
1547 TObjArray* dest = trc->GetLabel()>=0 ? fCHistoArrCorr : fCHistoArrFake;
781c79ac 1548 if (dest && fTrackPhaseID>kClus2Tracks) {
08419930 1549 //
1550 double htrPt = tmpTr.Pt();
9cdcba2c 1551 double dy = p[0]-tmpTr.GetY();
1552 double dz = p[1]-tmpTr.GetZ();
781c79ac 1553 ((TH2*)dest->At( GetHistoID(ilrA,kHResY,0,fTrackPhaseID) ))->Fill(htrPt,dy);
1554 ((TH2*)dest->At( GetHistoID(ilrA,kHResZ,0,fTrackPhaseID) ))->Fill(htrPt,dz);
76390254 1555 double errY = tmpTr.GetSigmaY2();
1556 double errZ = tmpTr.GetSigmaZ2();
781c79ac 1557 if (errY>0) ((TH2*)dest->At( GetHistoID(ilrA,kHResYP,0,fTrackPhaseID) ))->Fill(htrPt,dy/Sqrt(errY));
1558 if (errZ>0) ((TH2*)dest->At( GetHistoID(ilrA,kHResZP,0,fTrackPhaseID) ))->Fill(htrPt,dz/Sqrt(errZ));
1559 ((TH2*)dest->At( GetHistoID(ilrA,kHChi2Cl,0,fTrackPhaseID) ))->Fill(htrPt,chi2cl);
76390254 1560 }
1561#endif
9cdcba2c 1562 //
1563 if ( !tmpTr.Update(p,cov) ) {
e1ef49ad 1564#ifdef _ITSU_DEBUG_
846c381e 1565 AliESDtrack* trESD = fCurrHyp->GetESDTrack();
9cdcba2c 1566 AliDebug(2,Form("Failed on Update | ESDtrack#%d (MClb:%d)",trESD->GetID(),trESD->GetTPCLabel()));
e1ef49ad 1567#endif
9cdcba2c 1568 return -1;
e7d83d38 1569 }
e1ef49ad 1570#ifdef _ITSU_DEBUG_
b515ad7e 1571 if (AliDebugLevelClass()>1) {
9cdcba2c 1572 printf("AfterRefit: "); tmpTr.Print();
1573 }
e1ef49ad 1574#endif
5e3587be 1575 if (++nclFit==nCl && stopCond<0) {
7ce7f445 1576 *trc = (AliKalmanTrack&)tmpTr;
9cdcba2c 1577 return chi2; // it was requested to not propagate after last update
b515ad7e 1578 }
8b16dbae 1579 }
c51c3124 1580 //
1581 }
8b16dbae 1582 // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
1583 // Still, try to go as close as possible to rDest.
1584 //
70cb7fe4 1585 if (lrStart!=lrStop) {
e7d83d38 1586 if (!TransportToLayer(&tmpTr,lrStart,lrStop)) {
e1ef49ad 1587#ifdef _ITSU_DEBUG_
2947434b 1588 AliDebug(1,Form("Failed on TransportToLayer %d->%d | ESDtrack#%d (MClb:%d), X=%f",lrStart,lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb,tmpTr.GetX()));
e1ef49ad 1589#endif
2947434b 1590 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
e7d83d38 1591 }
1592 if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) {
e1ef49ad 1593#ifdef _ITSU_DEBUG_
2947434b 1594 AliDebug(1,Form("Failed on GoToExitFromLayer %d | ESDtrack#%d (MClb:%d), X=%f",lrStop,fCurrESDtrack->GetID(),fCurrESDtrMClb,tmpTr.GetX()));
e1ef49ad 1595#endif
2947434b 1596 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
9cdcba2c 1597 }
70cb7fe4 1598 }
1599 // go to the destination radius. Note that here we don't select direction to avoid precision problems
1600 if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) {
2947434b 1601 return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
70cb7fe4 1602 }
7ce7f445 1603 *trc = (AliKalmanTrack&)tmpTr;
e1ef49ad 1604#ifdef _ITSU_DEBUG_
15b02d69 1605 if (AliDebugLevelClass()>2) {
9cdcba2c 1606 printf("After refit (now at lr %d): ",lrStop); trc->Print();
15b02d69 1607 }
e1ef49ad 1608#endif
9cdcba2c 1609 return chi2;
c51c3124 1610}
dd2117a2 1611
1612//__________________________________________________________________
1613void AliITSUTrackerGlo::CookMCLabel(AliITSUTrackHyp* hyp)
1614{
1615 // build MC label
1616 //
1617 const int kMaxLbPerCl = 3;
1618 int lbID[kMaxLayers*6],lbStat[kMaxLayers*6];
1619 Int_t lr,nLab=0,nCl=0;
1620 AliITSUSeed *seed = hyp->GetWinner();
1621 while(seed) {
1622 int clID = seed->GetLrCluster(lr);
1623 if (clID>=0) {
1624 AliCluster *cl = fITS->GetLayerActive(lr)->GetCluster(clID);
1625 nCl++;
1626 for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
1627 int trLb = cl->GetLabel(imc);
704d9398 1628 if (trLb<0) break;
dd2117a2 1629 // search this mc track in already accounted ones
1630 int iLab;
1631 for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
1632 if (iLab<nLab) lbStat[iLab]++;
1633 else {
1634 lbID[nLab] = trLb;
1635 lbStat[nLab++] = 1;
1636 }
1637 } // loop over given cluster's labels
1638 }
1639 seed = (AliITSUSeed*)seed->GetParent();
1640 } // loop over clusters
1641 //
5785a9d8 1642 AliESDtrack* esdTr = hyp->GetESDTrack();
1643 int tpcLab = esdTr ? Abs(esdTr->GetTPCLabel()) : -kDummyLabel;
41343470 1644 if (nCl && nLab) {
dd2117a2 1645 int maxLab=0,nTPCok=0;
dd2117a2 1646 for (int ilb=nLab;ilb--;) {
1647 int st = lbStat[ilb];
1648 if (lbStat[maxLab]<st) maxLab=ilb;
1649 if (tpcLab==lbID[ilb]) nTPCok += st;
1650 }
1651 hyp->SetFakeRatio(1.-float(nTPCok)/nCl);
1652 hyp->SetLabel( nTPCok==nCl ? tpcLab : -tpcLab);
bdb92766 1653 hyp->SetITSLabel( lbStat[maxLab]==nCl ? lbID[maxLab] : -lbID[maxLab]); // winner label
dd2117a2 1654 return;
1655 }
1656 //
1657 hyp->SetFakeRatio(-1.);
5785a9d8 1658 hyp->SetLabel( -tpcLab );
dd2117a2 1659 hyp->SetITSLabel( kDummyLabel );
1660}
38997c3c 1661
1662//__________________________________________________________________
1663Bool_t AliITSUTrackerGlo::AddSeedBranch(AliITSUSeed* seed)
1664{
1665 // add new prolongation branch for the seed to temporary array. The seeds are sorted according to Compare f-n
6cd80116 1666 if (fNCandidatesAdded+fNBranchesAdded>=fLayerMaxCandidates ) {
1667 AliInfo(Form("Number of candidates at layer %d for current seed reached %d, increasing buffer",fCurrActLrID,fLayerMaxCandidates));
1668 fLayerMaxCandidates = 2*(fLayerMaxCandidates+1);
1669 AliITSUSeed** tmpArr = fLayerCandidates;
1670 fLayerCandidates = new AliITSUSeed*[fLayerMaxCandidates];
1671 memcpy(fLayerCandidates,tmpArr,(fNCandidatesAdded+fNBranchesAdded)*sizeof(AliITSUSeed*));
1672 delete tmpArr; // delete only array, not objects
38997c3c 1673 }
1674 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded]; // note: fNCandidatesAdded is incremented after adding all branches of current seed
1675 int slot=fNBranchesAdded++;
1676 for (int slotF=slot;slotF--;) { // slotF is always slot-1
1677 AliITSUSeed* si = branches[slotF];
1678 if (si->Compare(seed)<0) break; // found the last seed with better quality
1679 // otherwise, shift the worse seed to the next slot
1680 branches[slot] = si;
1681 slot = slotF; // slot should be slotF+1
1682 }
1683 // if needed, move worse seeds
1684 branches[slot] = seed;
1685 return kTRUE;
1686 //
1687}
1688
1689//__________________________________________________________________
1690void AliITSUTrackerGlo::ValidateAllowedBranches(Int_t acceptMax)
1691{
1692 // keep allowed number of branches for current seed and disable extras
70cb7fe4 1693 AliCodeTimerAuto("",0);
38997c3c 1694 int nb = Min(fNBranchesAdded,acceptMax);
1695 // if (nb<fNBranchesAdded) printf("ValidateAllowedBranches: %d of %d (%d) on lr %d\n",nb,fNBranchesAdded,fNCandidatesAdded,ilr);
1696 // disable unused branches
1697 AliITSUSeed** branches = &fLayerCandidates[fNCandidatesAdded];
08419930 1698 //
1699#ifdef _ITSU_TUNING_MODE_
1700 for (int ib=0;ib<fNBranchesAdded;ib++) branches[ib]->SetOrdBranch(ib);
c03e4f8a 1701#endif
bdb92766 1702 //
08419930 1703 for (int ib=nb;ib<fNBranchesAdded;ib++) MarkSeedFree(branches[ib]);
38997c3c 1704 fNCandidatesAdded += nb; // update total candidates counter
1705 fNBranchesAdded = 0; // reset branches counter
1706 //
1707}
1708
1709//__________________________________________________________________
1710void AliITSUTrackerGlo::ValidateAllowedCandidates(Int_t ilr, Int_t acceptMax)
1711{
1712 // transfer allowed number of branches to hypothesis container
1713 //
9cdcba2c 1714 AliCodeTimerAuto("",0);
38997c3c 1715 // sort candidates in increasing order of chi2
6cd80116 1716 static int lastSize = 0;
1717 static int *index = 0;
1718 static float *chi2 = 0;
1719 if (fLayerMaxCandidates>lastSize) {
1720 lastSize = fLayerMaxCandidates;
1721 delete[] index;
1722 delete[] chi2;
1723 index = new int[lastSize];
1724 chi2 = new float[lastSize];
1725 }
38997c3c 1726 for (int i=fNCandidatesAdded;i--;) chi2[i] = fLayerCandidates[i]->GetChi2GloNrm();
1727 Sort(fNCandidatesAdded,chi2,index,kFALSE);
1728 //
9cdcba2c 1729 int nacc=0,nb=0;
097e4cdf 1730 if (ilr>fCurrTrackCond->GetActiveLrInner()) { // just take 1st acceptMax candidates
9cdcba2c 1731 nb = Min(fNCandidatesAdded,acceptMax);
1732 for (int ib=0;ib<nb;ib++) AddProlongationHypothesis(fLayerCandidates[index[ib]],ilr);
1733 }
1734 else { // for innermost layer test ITS SA backward chi2 and matching to TPC
1735 AliITSUSeed* wn0 = fCurrHyp->GetWinner(); // in principle, should be NULL
1736 for (nb=0;nb<fNCandidatesAdded;nb++) {
1737 AliITSUSeed* sd = fLayerCandidates[index[nb]];
1738 fCurrHyp->SetWinner(sd);
1739 if (!CheckBackwardMatching(sd)) MarkSeedFree(sd); // discard bad backward tracks
1740 else {
1741 AddProlongationHypothesis(sd,ilr);
1742 if (++nacc==acceptMax) {nb++; break;} // the rest will be discarded
1743 }
1744 }
1745 fCurrHyp->SetWinner(wn0); // restore original winner (NULL?)
1746 }
38997c3c 1747 // disable unused candidates
08419930 1748 for (int ib=nb;ib<fNCandidatesAdded;ib++) MarkSeedFree(fLayerCandidates[index[ib]]);
38997c3c 1749 fNCandidatesAdded = 0; // reset candidates counter
1750 //
1751}
1752
9cdcba2c 1753//__________________________________________________________________
1754Bool_t AliITSUTrackerGlo::CheckBackwardMatching(AliITSUSeed* seed)
1755{
1756 // check seed backward propagation chi2 and matching to TPC
1757 double bz0 = GetBz();
70cb7fe4 1758 double rDest = fITS->GetRITSTPCRef(); // reference radius for comparison
7ce7f445 1759 AliITSUTrackHyp trback;
1760 trback.AliExternalTrackParam::operator=(*seed);
9cdcba2c 1761 trback.ResetCovariance(10000);
5e3587be 1762 int nclFit = 0;
1763 double chi2sa = RefitTrack(&trback,rDest,nclFit,1);
9cdcba2c 1764 if (chi2sa<0) return kFALSE;
5e3587be 1765 int ndf = nclFit*2-5;
9cdcba2c 1766 if (ndf>0) chi2sa /= ndf;
5e3587be 1767 if (chi2sa>fCurrTrackCond->GetMaxITSSAChi2(nclFit)) return kFALSE;
9cdcba2c 1768 //
1769 // relate to TPC track at outer layer
1770 AliExternalTrackParam* tpcSeed = fCurrHyp->GetTPCSeed();
1771 if (!trback.Rotate(tpcSeed->GetAlpha()) || !trback.PropagateParamOnlyTo(tpcSeed->GetX(),bz0)) {
e1ef49ad 1772#ifdef _ITSU_DEBUG_
9cdcba2c 1773 if (AliDebugLevelClass()>+1 && !seed->ContainsFake()) {
1774 AliInfo(Form("Failed to propagate seed to TPC track @ X:%.3f Alpha:%.3f",tpcSeed->GetX(),tpcSeed->GetAlpha()));
1775 seed->Print("etp");
1776 trback.Print();
1777 }
e1ef49ad 1778#endif
9cdcba2c 1779 return kFALSE;
1780 }
1781 double chi2Match = trback.GetPredictedChi2(tpcSeed)/5;
1782 if (chi2Match>fCurrTrackCond->GetMaxITSTPCMatchChi2()) return kFALSE;
1783 //
1784 seed->SetChi2ITSSA(chi2sa);
1785 seed->SetChi2ITSTPC(chi2Match);
1786 return kTRUE;
1787}
1788
6cd80116 1789//__________________________________________________________________
1790void AliITSUTrackerGlo::SaveReducedHypothesesTree(AliITSUTrackHyp* dest)
1791{
1792 // remove those hypothesis seeds which dont lead to candidates at final layer
1793 //
1794 // 1st, flag the seeds to save
0e40ce66 1795 int lr0 = fCurrTrackCond->GetActiveLrInner();
6cd80116 1796 for (int isd=0;isd<fCurrHyp->GetNSeeds(lr0);isd++) {
1797 AliITSUSeed* seed = fCurrHyp->RemoveSeed(lr0,isd);
1798 if (!seed) continue;
1799 seed->FlagTree(AliITSUSeed::kSave);
1800 dest->AddSeed(seed,lr0);
1801 }
0e40ce66 1802 for (int ilr=0;ilr<fNLrActive;ilr++) {
6cd80116 1803 int nsd = fCurrHyp->GetNSeeds(ilr);
1804 for (int isd=0;isd<nsd;isd++) {
1805 AliITSUSeed* seed = fCurrHyp->RemoveSeed(ilr,isd);
1806 if (!seed) continue; // already discarded or saved
1807 if (seed->IsSaved()) dest->AddSeed(seed,ilr);
1808 else MarkSeedFree(seed);
1809 }
1810 }
1811 //
1812 // AliInfo(Form("SeedsPool: %d, BookedUpTo: %d, free: %d",fSeedsPool.GetSize(),fSeedsPool.GetEntriesFast(),fNFreeSeeds));
1813}
1814
0e40ce66 1815//__________________________________________________________________
1816void AliITSUTrackerGlo::CleanHypothesis(AliITSUTrackHyp* hyp)
1817{
1818 // clean hypothesis
1819 hyp->SetWinner(0);
1820 for (int ilr=0;ilr<fNLrActive;ilr++) {
1821 int nsd = hyp->GetNSeeds(ilr);
1822 for (int isd=0;isd<nsd;isd++) {
1823 AliITSUSeed* seed = hyp->RemoveSeed(ilr,isd);
1824 if (!seed) continue; // already discarded or saved
1825 MarkSeedFree(seed);
1826 }
1827 }
1828}
1829
53870004 1830//__________________________________________________________________
1831void AliITSUTrackerGlo::FlagSplitClusters()
1832{
1833 // set special bit on split clusters using MC info
6cd80116 1834 for (int ilr=fNLrActive;ilr--;) {
53870004 1835 int nsplit=0;
1836 AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
1837 for (int isn=lr->GetNSensors();isn--;) {
1838 AliITSURecoSens* sens = lr->GetSensor(isn);
1839 int nCl = sens->GetNClusters();
1840 if (!nCl) continue;
1841 int cl0 = sens->GetFirstClusterId();
1842 for (int icl=nCl;icl--;) {
1843 AliITSUClusterPix *cl = (AliITSUClusterPix*)lr->GetCluster(cl0+icl);
1844 for (int icl1=icl;icl1--;) {
1845 AliITSUClusterPix *cl1 = (AliITSUClusterPix*)lr->GetCluster(cl0+icl1);
1846 if (cl->HasCommonTrack(cl1)) {
1847 if (!cl->IsSplit()) nsplit++;
1848 if (!cl1->IsSplit()) nsplit++;
1849 cl->SetSplit();
1850 cl1->SetSplit();
1851 }
1852 }
1853 }
1854 }
1855 AliInfo(Form("%4d out of %4d clusters are split",nsplit,lr->GetNClusters()));
1856 }
1857 //
1858}
1859
1860//__________________________________________________________________
1861Bool_t AliITSUTrackerGlo::ContainsSplitCluster(const AliITSUSeed* seed, Int_t maxSize)
1862{
1863 // check if the seed contains split cluster with size < maxSize
1864 int lrID,clID;
1865 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1866 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1867 if (cl->IsSplit() && cl->GetNPix()<maxSize ) return kTRUE;
1868 }
1869 return seed->GetParent() ? ContainsSplitCluster((AliITSUSeed*)seed->GetParent(),maxSize) : kFALSE;
1870 //
1871}
1872
1873//__________________________________________________________________
1874void AliITSUTrackerGlo::PrintSeedClusters(const AliITSUSeed* seed, Option_t* option)
1875{
1876 // print seeds clusters
1877 int lrID,clID;
1878 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1879 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1880 cl->Print(option);
1881 }
1882 if (seed->GetParent()) PrintSeedClusters((AliITSUSeed*)seed->GetParent(), option);
1883 //
1884}
1885
c03e4f8a 1886//__________________________________________________________________
08419930 1887void AliITSUTrackerGlo::FlagSeedClusters(const AliITSUSeed* seed, Bool_t flg, UShort_t ref)
c03e4f8a 1888{
1889 // mark used clusters
1890 int lrID,clID;
1891 while (seed) {
08419930 1892 if ( (clID=seed->GetLrCluster(lrID))>=0 ) {
1893 AliITSUClusterPix* cl = (AliITSUClusterPix*)fITS->GetLayerActive(lrID)->GetCluster(clID);
1894 if (ref) { // do we need to set or delete cluster->track ref?
1895 if (flg) {
1896 if (!cl->GetRecoInfo()) cl->SetRecoInfo(ref); // set ref only if cluster already does not refer to other track, inc.counter
1897 }
1898 else {
1899 if (cl->GetRecoInfo()==ref) cl->SetRecoInfo(0); // unset reference only if it refers to ref, decrease counter
1900 }
1901 }
1902 }
c03e4f8a 1903 seed = (AliITSUSeed*)seed->GetParent();
1904 }
1905 //
1906}
1907
d99994b2 1908//__________________________________________________________________
1909void AliITSUTrackerGlo::CheckClusterUsage()
1910{
1911 // check cluster usage info
eb18cdc6 1912 printf("ClusterUsage at pass %d\n",fCurrPassID);
d99994b2 1913 for (int ilr=0;ilr<fNLrActive;ilr++) {
1914 AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
1915 int ncl = lr->GetNClusters();
1916 int nclUsed = 0;
1917 int nclUs0CntShr1 = 0;
1918 int nclUs1Cnt0 = 0;
1919 int nclUs1Shr1 = 0;
1920 int nclUseL = 0;
1921 for (int icl=0;icl<ncl;icl++) {
1922 AliITSUClusterPix *cl = (AliITSUClusterPix*)lr->GetCluster(icl);
1923 int usc = cl->GetClUsage();
1924 if (usc<1) {
1925 if (cl->IsClusterUsed() || cl->IsClusterShared()) {
1926 nclUs0CntShr1++;
1927 printf("Lr%d Cl%4d: UseCnt=0 but UseFlg=%d UseShr=%d\n",ilr,icl,cl->IsClusterUsed(),cl->IsClusterShared());
1928 }
1929 continue;
1930 }
1931 nclUsed++;
1932 if (!cl->IsClusterUsed()) {
1933 nclUs1Cnt0++;
1934 printf("Lr%d Cl%4d: UseCnt=%d but UseFlg=%d UseShr=%d\n",ilr,icl,usc, cl->IsClusterUsed(),cl->IsClusterShared());
1935 }
1936 if (cl->IsClusterShared()) {
1937 nclUs1Shr1++;
1938 printf("Lr%d Cl%4d: UseCnt=%d but UseFlg=%d UseShr=%d\n",ilr,icl,usc, cl->IsClusterUsed(),cl->IsClusterShared());
1939 }
1940 if (usc>1) {
1941 nclUseL++;
1942 printf("Lr%d Cl%4d: UseCnt=%d!! but UseFlg=%d UseShr=%d\n",ilr,icl,usc, cl->IsClusterUsed(),cl->IsClusterShared());
1943 }
1944 }
1945 printf("ClUsage Lr%d: Ncl=%5d Nusd=%5d (%5.2f) Use0CS1:%4d Use1C0:%4d Use1S1:%4d UseL:%d\n",
1946 ilr,ncl,nclUsed,ncl? float(nclUsed)/ncl : 0, nclUs0CntShr1,nclUs1Cnt0,nclUs1Shr1,nclUseL);
1947 }
1948}
1949
1950
08419930 1951#ifdef _ITSU_TUNING_MODE_
53870004 1952//__________________________________________________________________
781c79ac 1953TObjArray* AliITSUTrackerGlo::BookControlHistos(const char* pref)
53870004 1954{
1955 // book special control histos
08419930 1956 TString prefS = pref;
1957 prefS.ToLower();
781c79ac 1958 TObjArray* dest = new TObjArray;
08419930 1959 dest->SetOwner(kTRUE);
1960 //
bdb92766 1961 const int kNResDef=7;
1962 const double kResDef[kNResDef]={0.05,0.05,0.3, 0.05,1,0.5,1.5};
bdb92766 1963 const double ptMax=10;
1964 const double plMax=10;
1965 const double chiMax=100;
1966 const int nptbins=50;
1967 const int nresbins=400;
1968 const int nplbins=50;
1969 const int nchbins=200;
1970 const int maxBr = 15;
1971 const int maxCand = 200;
1972 TString ttl;
781c79ac 1973 for (int phase=0;phase<kNTrackingPhases;phase++) {
1974 for (int pass=0;pass<AliITSUReconstructor::GetRecoParam()->GetNTrackingConditions();pass++) {
bdb92766 1975 //
781c79ac 1976 for (int ilr=0;ilr<fNLrActive;ilr++) {
1977 //
1978 // ----------------- These are histos to be filled in Cluster2Tracks of each pass.
1979 // PropagateBack and RefitInward will be stored among the histos of 1st pass
1980 if (pass>0 && phase!=kClus2Tracks) continue;
1981 //
1982 double mxdf = ilr>=kNResDef ? kResDef[kNResDef-1] : kResDef[ilr];
1983 ttl = Form("Pass%d_S%d_residY%d_%s",pass,phase,ilr,pref);
1984 TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1985 dest->AddAtAndExpand(hdy,GetHistoID(ilr,kHResY,pass,phase));
1986 hdy->SetDirectory(0);
1987 //
1988 ttl = Form("Pass%d_S%d_residYPull%d_%s",pass,phase,ilr,pref);
1989 TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
1990 dest->AddAtAndExpand(hdyp,GetHistoID(ilr,kHResYP,pass,phase));
1991 hdyp->SetDirectory(0);
1992 //
1993 ttl = Form("Pass%d_S%d_residZ%d_%s",pass,phase,ilr,pref);
1994 TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nresbins,-mxdf,mxdf);
1995 dest->AddAtAndExpand(hdz,GetHistoID(ilr,kHResZ,pass,phase));
1996 hdz->SetDirectory(0);
1997 //
1998 ttl = Form("Pass%d_S%d_residZPull%d_%s",pass,phase,ilr,pref);
1999 TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax,nplbins,-plMax,plMax);
2000 hdzp->SetDirectory(0);
2001 dest->AddAtAndExpand(hdzp,GetHistoID(ilr,kHResZP,pass,phase));
2002 //
2003 ttl = Form("Pass%d_S%d_chi2Cl%d_%s",pass,phase,ilr,pref);
2004 TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
2005 hchi->SetDirectory(0);
2006 dest->AddAtAndExpand(hchi,GetHistoID(ilr,kHChi2Cl,pass,phase));
2007 //
2008 // ------------------- These histos are filled for Clusters2Tracks only
2009 if (phase!=kClus2Tracks) continue;
2010 //
2011 ttl = Form("Pass%d_S%d_chi2Nrm%d_%s",pass,phase,ilr,pref);
2012 TH2F* hchiN = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
2013 hchiN->SetDirectory(0);
2014 dest->AddAtAndExpand(hchiN,GetHistoID(ilr,kHChi2Nrm,pass,phase));
2015 //
2016 ttl = Form("Pass%d_S%d_bestInBranch%d_%s",pass,phase,ilr,pref);
bdb92766 2017 TH2* hnbr = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxBr,-0.5,maxBr-0.5);
2018 hnbr->SetDirectory(0);
781c79ac 2019 dest->AddAtAndExpand(hnbr,GetHistoID(ilr,kHBestInBranch,pass,phase));
c03e4f8a 2020 //
781c79ac 2021 ttl = Form("Pass%d_S%d_bestInCands%d_%s",pass,phase,ilr,pref);
bdb92766 2022 TH2* hncn = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, maxCand,-0.5,maxCand-0.5);
2023 hncn->SetDirectory(0);
781c79ac 2024 dest->AddAtAndExpand(hncn,GetHistoID(ilr,kHBestInCand,pass,phase));
c03e4f8a 2025 //
781c79ac 2026 } // loop over layers
2027 //
2028 // these are special histos, filled not per layer but in the end of track fit in Clusters2Tracks in EVERY pass
2029 //
2030 if (phase==kClus2Tracks) {
2031 TH2* hchiMatch = 0;
2032 ttl = Form("Pass%d_Chi2Match_%s",pass,pref);
2033 hchiMatch = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins,0.,chiMax);
2034 hchiMatch->SetDirectory(0);
2035 dest->AddAtAndExpand(hchiMatch,GetHistoID(-1,kHChiMatch,pass,phase));
2036 //
2037 TH2* hchiSA = 0;
2038 ttl = Form("Pass%d_Chi2ITSSA_%s",pass,pref);
2039 hchiSA = new TH2F(ttl.Data(),ttl.Data(),nptbins,0,ptMax, nchbins/2,0.,chiMax/2);
2040 hchiSA->SetDirectory(0);
2041 dest->AddAtAndExpand(hchiSA,GetHistoID(-1,kHChiITSSA,pass,phase));
2042 //
53870004 2043 }
781c79ac 2044 } // loop over tracking passes
2045 }// loop over tracking phases
2046 //
2047 return dest;
53870004 2048}
781c79ac 2049
2050//__________________________________________________________________
2051Int_t AliITSUTrackerGlo::GetHistoID(Int_t lr, Int_t hid, Int_t pass, Int_t phase)
2052{
2053 // get id for the requested histo
2054 if (lr<0) lr=-1;
2055 lr++;
2056 return pass*kHistosPass + phase*kHistosPhase + lr*kMaxHID + hid;
2057 //
2058}
2059
53870004 2060#endif