Line data Source code
1 : #include "AliITSSAPTracker.h"
2 : #include "AliITSSAPLayer.h"
3 : #include "AliITSRecPoint.h"
4 : #include "AliGeomManager.h"
5 : #include "AliVParticle.h"
6 : #include "AliSymMatrix.h"
7 : //
8 : #include "AliRunLoader.h"
9 : #include "AliHeader.h"
10 : #include "AliGenEventHeader.h"
11 : #include "AliStack.h"
12 : #include <TParticle.h>
13 : #include <TParticlePDG.h>
14 : #include <TFile.h>
15 : #include "AliLog.h"
16 :
17 :
18 6 : ClassImp(AliITSSAPTracker)
19 :
20 : const Float_t AliITSSAPTracker::fgkZSpanITS[AliITSSAPTracker::kMaxLrITS] =
21 : { 36. ,14.1,14.1, 38., 22.2,29.7, 51. ,43.1,48.9};
22 :
23 : const Float_t AliITSSAPTracker::fgkRLayITS[AliITSSAPTracker::kMaxLrITS] =
24 : { 2.94, 3.9,7.6, 11.04, 15.0,23.9, 29.44 ,38.0,43.0};
25 :
26 : const Float_t AliITSSAPTracker::fgkRSpanITS[AliITSSAPTracker::kMaxLrITS] = // half span in R
27 : { 0.04, 0.5,0.5, 0.5, 0.8, 0.8, 0.5 ,0.6,0.6};
28 :
29 : const Int_t AliITSSAPTracker::fgkPassivLrITS[AliITSSAPTracker::kNLrPassive] =
30 : {AliITSSAPTracker::kLrBeamPime,AliITSSAPTracker::kLrShield1,AliITSSAPTracker::kLrShield2};
31 :
32 : const Int_t AliITSSAPTracker::fgkActiveLrITS[AliITSSAPTracker::kNLrActive] =
33 : {AliITSSAPTracker::kLrSPD1,AliITSSAPTracker::kLrSPD2,
34 : AliITSSAPTracker::kLrSDD1,AliITSSAPTracker::kLrSDD2,
35 : AliITSSAPTracker::kLrSSD1,AliITSSAPTracker::kLrSSD2};
36 :
37 : const Int_t AliITSSAPTracker::fgkLr2Active[AliITSSAPTracker::kMaxLrITS] = // conversion to active lr.
38 : {-1, 0, 1, -1, 2, 3, -1, 4, 5};
39 :
40 : const Float_t AliITSSAPTracker::fgkRhoLITS[AliITSSAPTracker::kMaxLrITS] = {
41 : 0.162802, 0.321960,0.354588, 0.274995, 0.193789,0.198168, 0.435372, 0.195828,0.226940};
42 :
43 : const Float_t AliITSSAPTracker::fgkX2X0ITS[AliITSSAPTracker::kMaxLrITS] = {
44 : 0.002757, 0.011660,0.012614, 0.006488, 0.007714,0.007916, 0.012689, 0.007849,0.009128};
45 :
46 :
47 : const Double_t AliITSSAPTracker::fgkClSystYErr2[AliITSSAPTracker::kNLrActive] =
48 : {0.0010*0.0010, 0.0030*0.0030, 0.0500*0.0500, 0.0500*0.0500, 0.0020*0.0020, 0.0020*0.0020};
49 :
50 : const Double_t AliITSSAPTracker::fgkClSystZErr2[AliITSSAPTracker::kNLrActive] =
51 : {0.0050*0.0050, 0.0050*0.0050, 0.0050*0.0050, 0.0050*0.0050, 0.1000*0.1000, 0.1000*0.1000};
52 :
53 :
54 : const Int_t AliITSSAPTracker::fgkLrDefBins[AliITSSAPTracker::kNLrActive][2] = // n bins in z, phi
55 : { {20,20}, {20,20}, {20,20}, {20,20}, {20,20}, {20,20} };
56 :
57 : const Float_t AliITSSAPTracker::fgkDefMass = 0.14;
58 : const Int_t AliITSSAPTracker::fgkDummyLabel = -3141593;
59 :
60 : #ifdef _TIMING_
61 : const char* AliITSSAPTracker::fgkSWNames[AliITSSAPTracker::kNSW] = {
62 : "Total"
63 : ,"Tracklets"
64 : ,"Tracks"
65 : ,"Vertex"
66 : };
67 : #endif
68 :
69 :
70 : //______________________________________________
71 0 : AliITSSAPTracker::AliITSSAPTracker() :
72 0 : fSPD2Discard()
73 0 : ,fTracklets()
74 0 : ,fSPD1Tracklet()
75 0 : ,fBlacklist(0)
76 0 : ,fPhiShift(0.0045)
77 0 : ,fSigThetaTracklet(0.025)
78 0 : ,fSigPhiTracklet(0.08)
79 0 : ,fChi2CutTracklet(1.5)
80 0 : ,fPhiShiftSc(0.)
81 0 : ,fDThetaTrackletSc(0)
82 0 : ,fDPhiTrackletSc(0)
83 0 : ,fBz(5.0)
84 0 : ,fMaxRSPDVtx(1.5)
85 0 : ,fDPhiTol(0.)
86 0 : ,fDThSig2Inv(0.)
87 0 : ,fDPhSig2Inv(0.)
88 : //
89 0 : ,fMinPt(0.3)
90 0 : ,fCurvMax(0)
91 0 : ,fZSPD2CutMin(1e9)
92 0 : ,fZSPD2CutMax(-1e9)
93 0 : ,fMaxChi2Tr2Cl(40)
94 0 : ,fAddErr2YspdVtx(0.02*0.02)
95 0 : ,fAddErr2ZspdVtx(0.04*0.04)
96 : //
97 0 : ,fMaxDRPhi(1.0)
98 0 : ,fMaxDZ(1.0)
99 : //
100 0 : ,fMissChi2Penalty(3)
101 0 : ,fMaxMissedLayers(1)
102 0 : ,fNTracks(0)
103 0 : ,fMaxTrackletsToRunTracking(99999)
104 0 : ,fMaxVtxIter(5)
105 0 : ,fStopScaleChange(0.8)
106 0 : ,fTracks()
107 0 : ,fTrackVertex()
108 0 : ,fFitVertex(kTRUE)
109 : //
110 0 : ,fSPDVertex(0)
111 : #ifdef _CONTROLH_
112 : ,fHTrackletMC(0),fHTrackletAll(0),fHTrackletFake(0),fHTrackMC(0),fHTrackAll(0),fHTrackFake(0)
113 : ,fHVtxDiffXY(0)
114 : ,fHVtxDiffXMlt(0),fHVtxDiffYMlt(0),fHVtxDiffZMlt(0)
115 : ,fHVtxPullXMlt(0),fHVtxPullYMlt(0),fHVtxPullZMlt(0)
116 : ,fHVtxMCSPDDiffXY(0)
117 : ,fHVtxMCSPDDiffXMlt(0),fHVtxMCSPDDiffYMlt(0),fHVtxMCSPDDiffZMlt(0)
118 : ,fHVtxMCSPDPullXMlt(0),fHVtxMCSPDPullYMlt(0),fHVtxMCSPDPullZMlt(0)
119 : ,fHChi2NDFvsPT(0),fHChi2vsNC(0)
120 : ,fHVtxMltRef(0),fHVtxOKMlt(0),fHVtxDiffZ(0),fHVtxMCSPDDiffZ(0)
121 : #endif
122 0 : {
123 : // def. c-tor
124 0 : for (int i=kNLrActive;i--;) fLayers[i] = 0;
125 0 : }
126 :
127 : //______________________________________________
128 : AliITSSAPTracker::~AliITSSAPTracker()
129 0 : {
130 : // d-tor
131 0 : for (int i=0;i<kNLrActive;i++) delete fLayers[i];
132 0 : }
133 :
134 : //______________________________________________
135 : void AliITSSAPTracker::Init()
136 : {
137 : // init tracker
138 : //
139 0 : if (!AliGeomManager::GetGeometry()) {
140 0 : AliGeomManager::LoadGeometry("geometry.root");
141 0 : AliGeomManager::ApplyAlignObjsFromCDB("ITS");
142 0 : }
143 : //
144 0 : for (int i=0;i<kNLrActive;i++) {
145 0 : int iAct = fgkActiveLrITS[i];
146 0 : fLayers[i] = new AliITSSAPLayer(i,fgkZSpanITS[iAct]+1,fgkLrDefBins[i][0],fgkLrDefBins[i][1]);
147 0 : fSkipLayer[i] = kFALSE;
148 0 : fNSigma2[i] = 7*7;
149 0 : fYToler2[i] = 0.2*0.2;
150 0 : fZToler2[i] = 0.2*0.2;
151 0 : fChi2TotCut[i] = 0;
152 : }
153 0 : fMaxDRPhi = 1.;
154 0 : fMaxDZ = 1.;
155 0 : fChi2TotCut[1] = 40; // 2 cl+vtx -> NDF=1
156 0 : fChi2TotCut[2] = 40;
157 0 : fChi2TotCut[3] = 30;
158 0 : fChi2TotCut[4] = 35;
159 0 : fChi2TotCut[5] = 40;
160 : //
161 0 : fMissChi2Penalty = 3;
162 0 : fMaxMissedLayers = 1;
163 : //
164 : // auxialary precalculated variables
165 0 : if (fChi2CutTracklet<0.1) fChi2CutTracklet = 0.1;
166 0 : double scl = TMath::Sqrt(fChi2CutTracklet);
167 0 : fDThetaTrackletSc = fSigThetaTracklet*scl;
168 0 : fDPhiTrackletSc = fSigPhiTracklet*scl;
169 : //
170 0 : fDThSig2Inv = 1./(fSigThetaTracklet*fSigThetaTracklet);
171 0 : fDPhSig2Inv = 1./(fSigPhiTracklet*fSigPhiTracklet);
172 : //
173 0 : fBlacklist = new TBits(100*100);
174 : //
175 0 : SetMaxVtxIter();
176 0 : SetStopScaleChange();
177 : //
178 : #ifdef _TIMING_
179 : for (int i=kNSW;i--;) {
180 : fSW[i].Stop();
181 : fSW[i].Reset();
182 : }
183 : #endif
184 : //
185 : #ifdef _CONTROLH_
186 : BookHistos();
187 : #endif
188 0 : }
189 :
190 : //______________________________________________
191 : void AliITSSAPTracker::ProcessEvent()
192 : {
193 : // do full reconstruction
194 : #ifdef _TIMING_
195 : fSW[kSWTotal].Start(0);
196 : fSW[kSWTracklets].Start(0);
197 : #endif
198 : //
199 0 : fNTracks = 0;
200 0 : FindTracklets();
201 : //
202 0 : if (GetNTracklets()>fMaxTrackletsToRunTracking) return;
203 : //
204 : #ifdef _TIMING_
205 : fSW[kSWTracklets].Stop();
206 : fSW[kSWTracks].Start(0);
207 : #endif
208 : //
209 0 : Tracklets2Tracks();
210 0 : RefitInward();
211 : #ifdef _TIMING_
212 : fSW[kSWTracks].Stop();
213 : fSW[kSWVertex].Start(0);
214 : #endif
215 0 : if (fFitVertex) {
216 0 : if (FitTrackVertex()) {
217 : #ifdef _DEBUG_
218 : printf("FittedVertex: "); fTrackVertex.Print(();
219 : printf("SPD Vertex: "); fSPDVertex->Print();
220 : #endif
221 : }
222 0 : }
223 : //
224 : //
225 : #ifdef _TIMING_
226 : fSW[kSWVertex].Stop();
227 : fSW[kSWTotal].Stop();
228 : PrintTiming();
229 : #endif
230 : //
231 : #ifdef _CONTROLH_
232 : FillRecoStat();
233 : #endif
234 : /*
235 : PrintTracklets();
236 : PrintTracks();
237 : if (fSPDVertex) {printf("SPDvtx: "); fSPDVertex->Print();}
238 : printf("TRKVtx: "); fTrackVertex.Print();
239 : */
240 0 : }
241 :
242 :
243 : //______________________________________________
244 : void AliITSSAPTracker::Clear(Option_t*)
245 : {
246 : // reset event info
247 0 : ClearTracklets();
248 0 : ClearTracks();
249 0 : for (int i=kNLrActive;i--;) {
250 0 : fNClusters[i] = 0;
251 0 : if (fLayers[i]) fLayers[i]->Clear();
252 : }
253 0 : }
254 :
255 : //______________________________________________
256 : void AliITSSAPTracker::ClearTracklets()
257 : {
258 : // reset tracklets info
259 0 : fSPD2Discard.clear();
260 0 : fTracklets.clear();
261 0 : fSPD1Tracklet.clear();
262 0 : if (fBlacklist) fBlacklist->ResetAllBits();
263 0 : }
264 :
265 :
266 : //______________________________________________
267 : void AliITSSAPTracker::AddCluster(AliITSRecPoint* cl)
268 : {
269 : // add cluster to corresponding layer
270 0 : if (!cl->Misalign()) AliWarning("Can't misalign this cluster !");
271 0 : fLayers[cl->GetLayer()]->AddCluster(cl);
272 0 : }
273 :
274 : //______________________________________________
275 : Bool_t AliITSSAPTracker::FindTracklets()
276 : {
277 : // find SPD tracklets
278 : //
279 0 : if (!fSPDVertex) {
280 : // AliInfo("No SPD vertex set");
281 0 : return kFALSE;
282 : }
283 0 : float rv2 = fSPDVertex->GetX()*fSPDVertex->GetX()+fSPDVertex->GetY()*fSPDVertex->GetY();
284 0 : if (rv2>fMaxRSPDVtx*fMaxRSPDVtx) {
285 : #ifdef _DEBUG_
286 : AliInfo("SPD vertex is too far from beam line");
287 : fSPDVertex->Print();
288 : #endif
289 0 : return kFALSE;
290 : }
291 0 : fPhiShiftSc = fPhiShift*TMath::Abs(fBz/5.0);
292 0 : fDPhiTol = fDPhiTrackletSc + fPhiShiftSc;
293 : //
294 0 : AliITSSAPLayer &spdL1 = *fLayers[kALrSPD1];
295 0 : AliITSSAPLayer &spdL2 = *fLayers[kALrSPD2];
296 0 : spdL1.SortClusters(fSPDVertex);
297 0 : spdL2.SortClusters(fSPDVertex);
298 0 : fNClusters[0] = spdL1.GetNClusters();
299 0 : fNClusters[1] = spdL2.GetNClusters();
300 : //
301 0 : if (fNClusters[0]<1 || fNClusters[1]<1) return kFALSE;
302 : //
303 0 : fSPD2Discard.resize(fNClusters[1]);
304 0 : fSPD1Tracklet.resize(fNClusters[0]);
305 : //
306 0 : fBlacklist->SetBitNumber(TMath::Max(fNClusters[0]*fNClusters[1],10000),kFALSE); // to reserve the space
307 : //
308 : int nfound;
309 0 : do {
310 : nfound = 0;
311 0 : for (int icl2=fNClusters[1];icl2--;) if (!fSPD2Discard[icl2]) nfound += AssociateClusterOfL2(icl2);
312 0 : } while(nfound);
313 : //
314 0 : for (int itr=GetNTracklets();itr--;) CookLabel(fTracklets[itr]);
315 : //
316 : return kTRUE;
317 0 : }
318 :
319 : //______________________________________________
320 : Int_t AliITSSAPTracker::AssociateClusterOfL2(int icl2)
321 : {
322 : // find SPD1 cluster matching to SPD2 cluster icl2
323 0 : AliITSSAPLayer &spdL1 = *fLayers[kALrSPD1];
324 0 : AliITSSAPLayer &spdL2 = *fLayers[kALrSPD2];
325 0 : AliITSSAPLayer::ClsInfo* cli2 = spdL2.GetClusterInfo(icl2);
326 : // expected z at SPD1
327 0 : float zV = fSPDVertex->GetZ();
328 0 : float z2 = cli2->z - zV;
329 0 : float tg2Inv = z2/cli2->r;
330 0 : float dzt = (1.+tg2Inv*tg2Inv)*fDThetaTrackletSc;
331 0 : float dz = dzt*fgkRLayITS[kLrSPD1] + TMath::Abs(tg2Inv)*fgkRSpanITS[kLrSPD1]; // uncertainty from dTheta and from Layer1 R spread
332 0 : float zL1 = zV + tg2Inv*fgkRLayITS[kLrSPD1]; // center of expected Z1
333 0 : int nsel1 = spdL1.SelectClusters(zL1-dz,zL1+dz, cli2->phi-fDPhiTol,cli2->phi+fDPhiTol);
334 0 : if (!nsel1) {
335 0 : fSPD2Discard[icl2] = true;
336 0 : return 0; // no candidates
337 : }
338 : float chiBest = 9999;
339 0 : SPDtracklet_t trk;
340 0 : trk.id1 = -1;
341 : int icl1,nCand=0;
342 0 : while ( (icl1=spdL1.GetNextClusterInfoID())!=-1) { // loop over matching clusters of lr1
343 0 : if (IsBlacklisted(icl1,icl2)) continue;
344 0 : AliITSSAPLayer::ClsInfo* cli1 = spdL1.GetClusterInfo(icl1);
345 0 : float z1 = cli1->z - zV;
346 0 : float tg1Inv = z1/cli1->r;
347 : //
348 0 : float dTheta = (tg2Inv-tg1Inv)/(1.+tg1Inv*tg1Inv); // fast check on theta
349 0 : if (TMath::Abs(dTheta)>fDThetaTrackletSc) continue;
350 : //
351 0 : float dPhi = cli1->phi - cli2->phi; // fast check on phi
352 0 : if (dPhi>TMath::Pi()) dPhi = TMath::TwoPi()-dPhi;
353 0 : else if (dPhi<-TMath::Pi()) dPhi += TMath::TwoPi();
354 0 : double dPhiS = TMath::Abs(dPhi)-fPhiShiftSc;
355 0 : if (TMath::Abs(dPhiS)>fDPhiTrackletSc) continue;
356 : //
357 0 : float chi2 = dTheta*dTheta*fDThSig2Inv + dPhiS*dPhiS*fDPhSig2Inv; // check final chi2
358 0 : if (chi2>1.) {
359 0 : Blacklist(icl1,icl2);
360 0 : continue;
361 : }
362 0 : nCand++;
363 0 : if (chi2>chiBest) continue;
364 : // check if cl1 is already associated with better
365 0 : trk.id1 = icl1;
366 0 : trk.id2 = icl2;
367 0 : trk.dtht = dTheta;
368 0 : trk.dphi = dPhi;
369 0 : trk.chi2 = chiBest = chi2;
370 0 : }
371 : //
372 0 : if (trk.id1!=-1) { // check if there is no better icl1 candidate for icl2
373 0 : int oldId = fSPD1Tracklet[trk.id1];
374 0 : if (!oldId) { // store new tracklet
375 0 : fTracklets.push_back(trk);
376 0 : fSPD1Tracklet[trk.id1] = fTracklets.size(); // refer from clusters to tracklet (id+1)
377 0 : fSPD2Discard[icl2] = true; // mark as used
378 0 : Blacklist(trk.id1,trk.id2);
379 0 : return 1;
380 : }
381 0 : SPDtracklet_t& oldTrk = (SPDtracklet_t&)fTracklets[--oldId];
382 0 : if (oldTrk.chi2 < trk.chi2) { // previous is better
383 0 : Blacklist(trk.id1,trk.id2); // shall we blacklist new combination?
384 0 : if (nCand==1) fSPD2Discard[icl2] = true; // there was just 1 candidate and it is discarded
385 0 : return 0;
386 : }
387 : // new combination is better, overwrite the old one with new one, marking old L2 cluster free
388 0 : fSPD2Discard[oldTrk.id2] = false; // mark as free
389 0 : fSPD2Discard[icl2] = true; // mark as used
390 0 : oldTrk = trk; // new combination is better, overwrite it with new one
391 0 : Blacklist(trk.id1,trk.id2);
392 0 : return 1;
393 : }
394 : //
395 0 : fSPD2Discard[icl2] = true; // no chance to find partner for this cluster
396 0 : return 0;
397 : //
398 0 : }
399 :
400 :
401 : //______________________________________________
402 : void AliITSSAPTracker::Tracklets2Tracks()
403 : {
404 : // try to extend tracklets to outer layers
405 0 : int nTrk = GetNTracklets();
406 0 : if (!nTrk) return;
407 : //
408 0 : CalcAuxTracking(); // RS??? do we need to repeat this?
409 : //
410 0 : for (int ila=kALrSDD1;ila<kNLrActive;ila++) {
411 0 : if (fSkipLayer[ila]) continue;
412 0 : fLayers[ila]->SortClusters(0);
413 0 : fNClusters[ila] = fLayers[ila]->GetNClusters();
414 0 : }
415 : //
416 0 : fTracks.resize(nTrk);
417 :
418 : //
419 0 : for (int itr=0;itr<nTrk;itr++) {
420 0 : SPDtracklet_t& trlet = fTracklets[itr];
421 : //
422 : #ifdef _DEBUG_
423 : printf("TestTracklet %d\t|",itr);
424 : int stat = GetTrackletMCTruth(trlet);
425 : //
426 : int nmiss=0;
427 : for (int i=2;i<kNLrActive;i++) {
428 : printf("%c", (stat&(0x1<<i)) ? '*':'-');
429 : if (!(stat&(0x1<<i))) nmiss++;
430 : }
431 : printf("|\n");
432 : PrintTracklet(itr);
433 : #endif
434 : //
435 0 : float zspd2 = fLayers[kALrSPD2]->GetClusterInfo(trlet.id2)->z;
436 0 : if (zspd2<fZSPD2CutMin || zspd2>fZSPD2CutMax) continue;
437 0 : ITStrack_t &track = fTracks[fNTracks];
438 0 : if (!CreateTrack(track, trlet)) continue;
439 0 : track.trackletID = itr;
440 : Bool_t res;
441 : #ifdef _DEBUG_
442 : double xyz[3];
443 : track.paramOut.GetXYZAt(0,fBz,xyz);
444 : printf("process track pt:%f XYZ: %+.4f %+.4f %+.4f\n",track.paramOut.Pt(),xyz[0],xyz[1],xyz[2]);
445 : #endif
446 0 : for (int lrID=kLrShield1;lrID<kMaxLrITS;lrID++) {
447 0 : res = FollowToLayer(track,lrID) && IsAcceptableTrack(track);
448 0 : if (!res) break;
449 : }
450 : #ifdef _DEBUG_
451 : printf("%s:%d\n",res ? "OK" : "Fail",nmiss<=fMaxMissedLayers);
452 : #endif
453 0 : if (!res) continue;
454 0 : track.paramOut.ResetBit(kInvalidBit); // flag that outward fit succeeded
455 0 : CookLabel(track);
456 0 : fNTracks++;
457 : //
458 0 : }
459 0 : }
460 :
461 : //______________________________________________
462 : Bool_t AliITSSAPTracker::IsAcceptableTrack(const AliITSSAPTracker::ITStrack_t& /*track*/) const
463 : {
464 : // check if the track is acceptable
465 0 : return kTRUE;
466 : }
467 :
468 : //______________________________________________
469 : void AliITSSAPTracker::PrintTrack(const AliITSSAPTracker::ITStrack_t& track) const
470 : {
471 : // print track info
472 0 : printf("Chi2 = %f for %d clusters. Tracklet %d\n",track.chi2,track.ncl,track.trackletID);
473 : //
474 0 : for (int ilr=0;ilr<kNLrActive;ilr++) {
475 0 : if (track.clID[ilr]<0) continue;
476 0 : AliITSRecPoint* cl = fLayers[ilr]->GetClusterSorted(track.clID[ilr]);
477 0 : printf("L%d #%4d ",ilr,track.clID[ilr]);
478 0 : for (int i=0;i<3;i++) printf("%d ",cl->GetLabel(i)); printf("\n");
479 0 : }
480 0 : track.paramOut.Print();
481 0 : track.paramInw.Print();
482 0 : }
483 :
484 : //______________________________________________
485 : void AliITSSAPTracker::PrintTracklets() const
486 : {
487 : // print traklets info
488 0 : int ntr = fTracklets.size();
489 0 : printf("NTracklets: %d\n",ntr);
490 0 : printf("Nspd1: %4d Nspd2: %4d, Ntracklets: %d\n",fNClusters[0],fNClusters[1],ntr);
491 0 : for (int itr=0;itr<ntr;itr++) PrintTracklet(itr);
492 : //
493 0 : }
494 :
495 : //______________________________________________
496 : void AliITSSAPTracker::PrintTracklet(Int_t itr) const
497 : {
498 : // print single tracklet
499 0 : const SPDtracklet_t* trk = &fTracklets[itr];
500 0 : AliITSRecPoint* cl1 = fLayers[kALrSPD1]->GetClusterSorted(trk->id1);
501 0 : AliITSRecPoint* cl2 = fLayers[kALrSPD2]->GetClusterSorted(trk->id2);
502 0 : AliITSSAPLayer::ClsInfo_t* cli0 = fLayers[kALrSPD1]->GetClusterInfo(trk->id1);
503 0 : printf("#%3d Phi:%+.3f Eta:%+.3f Dphi:%+.3f Dtht:%+.3f Chi2:%.3f | Lbl:",
504 0 : itr,cli0->phi,
505 0 : -TMath::Log(TMath::Tan(TMath::ATan2(cli0->r,cli0->z-fSPDVertex->GetZ())/2.)),
506 0 : trk->dphi,trk->dtht,trk->chi2);
507 : int lab=-1,lb = -1;
508 0 : for (int i=0;i<3;i++) if ( (lb=cl1->GetLabel(i))>=0 ) {if (lab<0)lab=lb; printf(" %5d",lb);} printf("|");
509 0 : for (int i=0;i<3;i++) if ( (lb=cl2->GetLabel(i))>=0 ) printf(" %5d",lb);
510 0 : printf("| ->%d\n",trk->label);
511 0 : lab = TMath::Abs(trk->label);
512 : //
513 : AliStack* stack = 0;
514 0 : AliRunLoader* rl = AliRunLoader::Instance();
515 0 : if (lab>=0 && rl && (stack=rl->Stack())) {
516 0 : TParticle* mctr = stack->Particle(lab);
517 0 : if (mctr) {
518 0 : TParticlePDG* mctrPDG = mctr->GetPDG();
519 0 : if (mctrPDG) {
520 0 : double qpt = mctrPDG->Charge()>0 ? mctr->Pt() : -mctr->Pt();
521 0 : printf("MCTrack: Prim:%d Vxyz: {%+.4f %+.4f %+.4f} 1/pt: %.3f tgl: %.3f\n",
522 0 : stack->IsPhysicalPrimary(lab),
523 0 : mctr->Vx(),mctr->Vy(),mctr->Vz(),
524 0 : TMath::Abs(qpt)>0 ? 1./qpt : 9999., TMath::Tan(TMath::Pi()/2. - mctr->Theta()));
525 0 : }
526 0 : }
527 0 : }
528 0 : }
529 :
530 :
531 : //______________________________________________
532 : void AliITSSAPTracker::PrintTracks() const
533 : {
534 : // print tracks info
535 0 : printf("NTracks: %d\n",fNTracks);
536 0 : for (int itr=0;itr<fNTracks;itr++) PrintTrack(fTracks[itr]);
537 : //
538 0 : }
539 :
540 :
541 : //______________________________________________
542 : void AliITSSAPTracker::CalcAuxTracking()
543 : {
544 : // precalculate auxilarry variables for tracking
545 : //
546 : // largest track curvature to search
547 : const double ztolerEdge = 1.0;
548 0 : fCurvMax = TMath::Abs(fBz*kB2C/fMinPt);
549 : double thMin =-1e9;
550 : double thMax = 1e9;
551 0 : for (int ilA=kNLrActive-1;ilA>kALrSPD2;ilA--) {
552 0 : if (GetSkipLayer(ilA)) continue;
553 0 : int ilr=fgkActiveLrITS[ilA];
554 0 : double r = fgkRLayITS[ilr] - fgkRSpanITS[ilr];
555 0 : double dz = fgkZSpanITS[ilr]+ztolerEdge+fDThetaTrackletSc*r;
556 0 : double ri = 1./r;
557 0 : double tmin= (-dz-fSPDVertex->GetZ())*ri;
558 0 : double tmax= ( dz-fSPDVertex->GetZ())*ri;
559 0 : if (tmin>thMin) thMin = tmin;
560 0 : if (tmax<thMax) thMax = tmax;
561 0 : }
562 : double r = fgkRLayITS[kLrSPD2] + fgkRSpanITS[kLrSPD2];
563 0 : fZSPD2CutMin = fSPDVertex->GetZ()+thMin*r; // min Z of SPD2 in tracklet to consider tracking
564 0 : fZSPD2CutMax = fSPDVertex->GetZ()+thMax*r; // max Z of SPD2 in tracklet to consider tracking
565 : //
566 0 : }
567 :
568 : //______________________________________________
569 : Bool_t AliITSSAPTracker::CreateTrack(AliITSSAPTracker::ITStrack_t& track,
570 : AliITSSAPTracker::SPDtracklet_t& trlet)
571 : {
572 : // create track seed from tracklet
573 : // init track
574 0 : track.label = trlet.label;
575 : //
576 0 : AliITSSAPLayer::ClsInfo_t *cli1=fLayers[kALrSPD1]->GetClusterInfo(trlet.id1);
577 0 : AliITSSAPLayer::ClsInfo_t *cli2=fLayers[kALrSPD2]->GetClusterInfo(trlet.id2);
578 0 : AliITSRecPoint *cl1=fLayers[kALrSPD1]->GetClusterUnSorted(cli1->index);
579 0 : AliITSRecPoint *cl2=fLayers[kALrSPD2]->GetClusterUnSorted(cli2->index);
580 0 : int det1 = cl1->GetVolumeId()-fLayers[kALrSPD1]->GetVIDOffset();
581 0 : int det2 = cl2->GetVolumeId()-fLayers[kALrSPD2]->GetVIDOffset();
582 0 : AliITSSAPLayer::ITSDetInfo_t& detInfo1 = fLayers[kALrSPD1]->GetDetInfo(det1);
583 0 : AliITSSAPLayer::ITSDetInfo_t& detInfo2 = fLayers[kALrSPD2]->GetDetInfo(det2);
584 : //
585 : // crude momentun estimate
586 0 : float dx=cli1->x-cli2->x,dy=cli1->y-cli2->y,d=TMath::Sqrt(dx*dx+dy*dy);
587 0 : float qptInv = fBz ? 2*TMath::Sin(cli2->phi-cli1->phi)/d/fBz/kB2C : 0; // positive particle goes anticlockwise in B+
588 : //
589 : // we initialize the seed in the tracking frame of 1st detector
590 0 : float xv= fSPDVertex->GetX()*detInfo1.cosTF + fSPDVertex->GetY()*detInfo1.sinTF;
591 0 : float yv=-fSPDVertex->GetX()*detInfo1.sinTF + fSPDVertex->GetY()*detInfo1.cosTF;
592 0 : float zv= fSPDVertex->GetZ();
593 0 : float par[5] = {yv, zv, (float)TMath::Sin(cli1->phi-detInfo1.phiTF), (cli1->z-zv)/cli1->r, qptInv};
594 0 : double covVtx[6];
595 0 : fSPDVertex->GetCovarianceMatrix(covVtx);
596 0 : float cov[15] = {float(covVtx[0]+covVtx[2] + fAddErr2YspdVtx),
597 0 : 0, float(covVtx[5] + fAddErr2ZspdVtx),
598 : 0,0,1,
599 : 0,0,0,1,
600 : 0,0,0,0,100*100};
601 0 : AliExternalTrackParam& param = track.paramOut;
602 0 : param.Set(xv, detInfo1.phiTF, par, cov);
603 0 : track.chi2 = 0; // chi2 at 1st two point is 0
604 : // go to 1st layer, ignoring the MS (errors are anyway not defined)
605 0 : if (!param.PropagateTo(detInfo1.xTF+cl1->GetX(), fBz)) return kFALSE;
606 0 : Double_t cpar0[2]={ cl1->GetY(), cl1->GetZ()};
607 0 : Double_t ccov0[3]={ cl1->GetSigmaY2() + GetClSystYErr2(kALrSPD1), 0., cl1->GetSigmaZ2() + GetClSystZErr2(kALrSPD1)};
608 0 : if (!param.Update(cpar0,ccov0)) return kFALSE;
609 0 : if (!param.CorrectForMeanMaterial(fgkX2X0ITS[kLrSPD1],-fgkRhoLITS[kLrSPD1],fgkDefMass)) return kFALSE;
610 : // go to 2nd layer
611 0 : if (!param.Rotate(detInfo2.phiTF)) return kFALSE;
612 0 : if (!param.PropagateTo(detInfo2.xTF+cl2->GetX(), fBz)) return kFALSE;
613 0 : Double_t cpar1[2]={ cl2->GetY(), cl2->GetZ()};
614 0 : Double_t ccov1[3]={ cl2->GetSigmaY2() + GetClSystYErr2(kALrSPD2), 0., cl2->GetSigmaZ2() + GetClSystZErr2(kALrSPD2)};
615 0 : track.chi2 += param.GetPredictedChi2(cpar1,ccov1);
616 0 : if (!param.Update(cpar1,ccov1)) return kFALSE;
617 : #ifdef _CONTROLH_
618 : FillTrackingControlHistos(1,track.label,¶m,cpar1,ccov1,cl2);
619 : #endif
620 : //
621 0 : track.clID[0] = trlet.id1;
622 0 : track.clID[1] = trlet.id2;
623 0 : track.clID[2] = track.clID[3] = track.clID[4] = track.clID[5] = -1;
624 0 : track.ncl = 2;
625 0 : track.nmiss=0;
626 : //
627 0 : param.SetBit(kInvalidBit); // flag that track is not yer refitted outward
628 0 : track.paramOut.SetBit(kInvalidBit); // flag that track was not refitter inward
629 0 : return kTRUE;
630 0 : }
631 :
632 : //______________________________________________
633 : Bool_t AliITSSAPTracker::CrossPassiveLayer(AliExternalTrackParam& param, Int_t lrID)
634 : {
635 : // cross the layer, applying mat. corrections
636 0 : double xStart=param.GetX();
637 0 : double xToGo = GetXatLabRLin(param,fgkRLayITS[lrID]);
638 0 : if (xToGo<0 || !param.PropagateTo(xToGo,fBz)) return kFALSE;
639 0 : double x2x0=fgkX2X0ITS[lrID],xrho=fgkRhoLITS[lrID];
640 0 : if (xStart<xToGo) xrho = -xrho; // inward propagation
641 0 : return param.CorrectForMeanMaterial(x2x0,xrho,fgkDefMass,kFALSE);
642 : //
643 0 : }
644 :
645 : //______________________________________________
646 : Bool_t AliITSSAPTracker::FollowToLayer(AliITSSAPTracker::ITStrack_t& track, Int_t lrID)
647 : {
648 : // take track to given layer, searching hits if needed and applying mat. corrections
649 0 : int lrIDA = fgkLr2Active[lrID]; // active layer ID
650 0 : if (lrIDA<0 || fSkipLayer[lrIDA]) return CrossPassiveLayer(track.paramOut,lrID);
651 : //
652 0 : AliExternalTrackParam trCopy(track.paramOut);
653 0 : double xToGo = GetXatLabRLin(trCopy,fgkRLayITS[lrID]); // aproximate X at lrID
654 0 : if (!trCopy.PropagateTo(xToGo,fBz)) return kFALSE;
655 0 : double xyz[3];
656 0 : trCopy.GetXYZ(xyz);
657 0 : double phi=TMath::ATan2(xyz[1],xyz[0]),z=trCopy.GetZ();
658 : // we need track errors in the plane nearly tangential to crossing point
659 0 : if (!trCopy.Rotate(phi)) return kFALSE;
660 0 : double drphi = TMath::Sqrt(trCopy.GetSigmaY2()*fNSigma2[lrIDA]+fYToler2[lrIDA]);
661 0 : if (drphi>fMaxDRPhi) drphi = fMaxDRPhi;
662 0 : double dphi = drphi/fgkRLayITS[lrID];
663 0 : double dz = TMath::Sqrt(trCopy.GetSigmaZ2()*fNSigma2[lrIDA]+fZToler2[lrIDA]);
664 0 : if (dz>fMaxDZ) dz = fMaxDZ;
665 0 : AliITSSAPLayer* lrA = fLayers[lrIDA];
666 0 : int nCl = lrA->SelectClusters(z-dz,z+dz,phi-dphi,phi+dphi);
667 : Bool_t updDone = kFALSE;
668 : //
669 : #ifdef _DEBUG_
670 : printf("at Lr%d, Ncl:%d ",lrIDA,nCl);
671 : trCopy.Print();
672 : #endif
673 : //
674 0 : if (nCl) {
675 : int icl,iclBest=-1;
676 0 : double chi2Best = fMaxChi2Tr2Cl;
677 : AliITSRecPoint* bestCl = 0;
678 0 : AliExternalTrackParam bestTr;
679 : //
680 : #ifdef _DEBUG_
681 : int iclt=0;
682 : #endif
683 0 : while ( (icl=lrA->GetNextClusterInfoID())!=-1) {
684 0 : AliITSSAPLayer::ClsInfo_t *cli = lrA->GetClusterInfo(icl);
685 0 : AliITSRecPoint *cl=lrA->GetClusterUnSorted(cli->index);
686 0 : int detId = cl->GetVolumeId()-lrA->GetVIDOffset();
687 0 : AliITSSAPLayer::ITSDetInfo_t& detInfo = lrA->GetDetInfo(detId);
688 0 : trCopy = track.paramOut;
689 0 : if (!trCopy.Propagate(detInfo.phiTF, detInfo.xTF+cl->GetX(), fBz)) continue;
690 0 : double cpar[2]={ cl->GetY(), cl->GetZ()};
691 0 : double ccov[3]={ cl->GetSigmaY2() + GetClSystYErr2(lrIDA) , 0., cl->GetSigmaZ2() + GetClSystZErr2(lrIDA)};
692 0 : double chi2cl = trCopy.GetPredictedChi2(cpar,ccov);
693 : //
694 : #ifdef _DEBUG_
695 : float clXYZ[3]; cl->GetGlobalXYZ(clXYZ);
696 : double trXYZ[3]; trCopy.GetXYZ(trXYZ);
697 : Float_t xCl, alphaCl;
698 : cl->GetXAlphaRefPlane(xCl,alphaCl);
699 : //
700 : printf("cl%d Chi2:%.2f Dyz: %+e %+e Err: %e %e %e |Lb:",iclt++,chi2cl,
701 : cl->GetY()-trCopy.GetY(),cl->GetZ()-trCopy.GetZ(),
702 : TMath::Sqrt(ccov[0]),ccov[1],TMath::Sqrt(ccov[2])); //TMP
703 : for (int j=0;j<3;j++) if (cl->GetLabel(j)>=0) printf(" %d",cl->GetLabel(j)); printf("\n");
704 : printf("CL: X:%.4f Alp:%+.4f XYZ: %+.4f %+.4f %+.4f\n",xCl,alphaCl,clXYZ[0],clXYZ[1],clXYZ[2]);
705 : printf("TR: X:%.4f Alp:%+.4f XYZ: %+.4f %+.4f %+.4f\n",detInfo.xTF,detInfo.phiTF,trXYZ[0],trXYZ[1],trXYZ[2]);
706 : trCopy.Print();
707 : #endif
708 : //
709 : #ifdef _CONTROLH_
710 : FillTrackingControlHistos(lrIDA,track.label,&trCopy,cpar,ccov,cl);
711 : #endif
712 : //
713 0 : if (chi2cl>fMaxChi2Tr2Cl) continue;
714 : // SaveCandidate(lrIDA,trCopy,chi2cl,icl); // RS: do we need this?
715 0 : if (chi2cl>chi2Best) continue;
716 : chi2Best = chi2cl;
717 : iclBest = icl;
718 : bestCl = cl;
719 0 : bestTr = trCopy;
720 0 : if (nCl==1) { // in absence of competitors, do the fit on spot
721 0 : if (!bestTr.Update(cpar,ccov)) return kFALSE;
722 : updDone = kTRUE;
723 0 : }
724 0 : }
725 : #ifdef _DEBUG_
726 : printf("Lr%d -> %f\n",lrIDA,chi2Best);
727 : #endif
728 : //
729 0 : if (bestCl) {
730 0 : if (!updDone) {
731 0 : double cpar[2]={ bestCl->GetY(), bestCl->GetZ()};
732 0 : double ccov[3]={ bestCl->GetSigmaY2(), 0., bestCl->GetSigmaZ2()}; // RS: add syst errors
733 0 : if (!bestTr.Update(cpar,ccov)) return kFALSE;
734 : updDone = kTRUE;
735 0 : }
736 0 : track.paramOut = bestTr;
737 0 : track.clID[lrIDA] = iclBest;
738 0 : track.ncl++;
739 0 : track.chi2 += chi2Best;
740 0 : }
741 0 : }
742 : //
743 0 : if (!updDone) {
744 0 : if (++track.nmiss > fMaxMissedLayers) return kFALSE;
745 0 : track.paramOut = trCopy;
746 0 : track.chi2 += fMissChi2Penalty;
747 0 : }
748 : //
749 : #ifdef _CONTROLH_
750 : int ndf = 2*track.ncl-5;
751 : if (ndf>0) {
752 : fHChi2vsNC->Fill(track.ncl,track.chi2);
753 : if (lrID==kNLrActive-1) fHChi2NDFvsPT->Fill(track.paramOut.Pt(),track.chi2/ndf);
754 : }
755 : #endif
756 0 : if (track.chi2 > GetChi2TotCut(track.ncl+1)) return kFALSE;
757 : //
758 0 : return track.paramOut.CorrectForMeanMaterial(fgkX2X0ITS[lrID],-fgkRhoLITS[lrID],fgkDefMass,kFALSE);
759 : //
760 0 : }
761 :
762 : //______________________________________________
763 : void AliITSSAPTracker::CookLabel(AliITSSAPTracker::ITStrack_t& track)
764 : {
765 : // cook mc label for the track
766 0 : track.label = fgkDummyLabel;
767 0 : if (!track.ncl) return;
768 : const int kMaxLbPerCl = 3;
769 0 : int lbID[kNLrActive*6],lbStat[kNLrActive*6];
770 : Int_t nLab=0;
771 0 : for (int i=kNLrActive;i--;) {
772 0 : int clid = track.clID[i];
773 0 : if (clid<0) continue;
774 0 : AliITSRecPoint* cl = fLayers[i]->GetClusterSorted(clid);
775 0 : for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
776 0 : int trLb = cl->GetLabel(imc);
777 0 : if (trLb<0) break;
778 : // search this mc track in already accounted ones
779 : int iLab;
780 0 : for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
781 0 : if (iLab<nLab) lbStat[iLab]++;
782 : else {
783 0 : lbID[nLab] = trLb;
784 0 : lbStat[nLab++] = 1;
785 : }
786 0 : } // loop over given cluster's labels
787 0 : } // loop over all clusters
788 : //
789 0 : if (nLab) {
790 : int maxLab=0;
791 0 : for (int ilb=nLab;ilb--;) if (lbStat[maxLab]<lbStat[ilb]) maxLab=ilb;
792 0 : track.label = lbStat[maxLab]==track.ncl ? lbID[maxLab] : -lbID[maxLab];
793 0 : }
794 : //
795 0 : }
796 :
797 : //______________________________________________
798 : void AliITSSAPTracker::CookLabel(AliITSSAPTracker::SPDtracklet_t& tracklet)
799 : {
800 : // cook mc label for the tracklet
801 0 : tracklet.label = fgkDummyLabel;
802 : const int kMaxLbPerCl = 3;
803 0 : int lbID[kNLrActive*6],lbStat[kNLrActive*6];
804 : Int_t nLab=0;
805 0 : for (int i=2;i--;) {
806 0 : int clid = i ? tracklet.id2 : tracklet.id1;
807 0 : AliITSRecPoint* cl = fLayers[i]->GetClusterSorted(clid);
808 0 : for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
809 0 : int trLb = cl->GetLabel(imc);
810 0 : if (trLb<0) break;
811 : // search this mc track in already accounted ones
812 : int iLab;
813 0 : for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
814 0 : if (iLab<nLab) lbStat[iLab]++;
815 : else {
816 0 : lbID[nLab] = trLb;
817 0 : lbStat[nLab++] = 1;
818 : }
819 0 : } // loop over given cluster's labels
820 : } // loop over all clusters
821 : //
822 0 : if (nLab) {
823 : int maxLab=0;
824 0 : for (int ilb=nLab;ilb--;) if (lbStat[maxLab]<lbStat[ilb]) maxLab=ilb;
825 0 : tracklet.label = lbStat[maxLab]==2 ? lbID[maxLab] : -lbID[maxLab];
826 0 : }
827 : //
828 0 : }
829 :
830 : //______________________________________________
831 : Double_t AliITSSAPTracker::GetXatLabRLin(AliExternalTrackParam& track, double r)
832 : {
833 : // X of track circle intersection in current tracking frame, neglecting the curvature
834 : // Solution of equation (x+d)^2+(y+b*d)^2 - r^2, where x,y are current coordinates of
835 : // track and d=X-x0. b = tg(phi)
836 : //double sn=tr.GetSnp();
837 0 : double sn=track.GetSnp();
838 0 : if (TMath::Abs(sn)>kAlmost1) return -999;
839 0 : double x=track.GetX(), y=track.GetY();
840 0 : double cs2=(1.-sn)*(1.+sn), tg=sn/TMath::Sqrt(cs2);
841 0 : double t0=x+tg*y, t1=x*x+y*y-r*r, det=t0*t0-t1/cs2;
842 0 : if (det<0) return -999; // does not touch circle
843 0 : det = TMath::Sqrt(det);
844 0 : return x+(det-t0)*cs2;
845 : //
846 0 : }
847 :
848 : //______________________________________________
849 : Int_t AliITSSAPTracker::GetTrackletMCTruth(AliITSSAPTracker::SPDtracklet_t& trlet) const
850 : {
851 : int status = 0;
852 0 : AliITSSAPLayer::ClsInfo_t *cli1=fLayers[kALrSPD1]->GetClusterInfo(trlet.id1);
853 0 : AliITSSAPLayer::ClsInfo_t *cli2=fLayers[kALrSPD2]->GetClusterInfo(trlet.id2);
854 0 : AliITSRecPoint *cl1=fLayers[kALrSPD1]->GetClusterUnSorted(cli1->index);
855 0 : AliITSRecPoint *cl2=fLayers[kALrSPD2]->GetClusterUnSorted(cli2->index);
856 : //
857 : int lab = -1;
858 : //
859 0 : for (int i=0;i<3;i++) {
860 0 : int lb1 = cl1->GetLabel(i);
861 0 : if (lb1<0) continue;
862 0 : for (int j=0;j<3;j++) {
863 0 : int lb2 = cl2->GetLabel(i);
864 0 : if (lb2<0) break;
865 0 : if (lb1==lb2) {lab = lb1; break;}
866 0 : }
867 0 : if (lab>=0) break;
868 0 : }
869 0 : if (lab<0) return 0;
870 : //
871 0 : for (int ila=kALrSDD1;ila<kNLrActive;ila++) {
872 0 : for (int icl=fNClusters[ila];icl--;) {
873 0 : AliITSRecPoint *cl=fLayers[ila]->GetClusterUnSorted(icl);
874 0 : for (int i=0;i<3;i++) {
875 0 : if (cl->GetLabel(i)<0) break;
876 0 : if (cl->GetLabel(i)==lab) {status |= 0x1<<ila; break;}
877 : }
878 0 : if (status & (0x1<<ila)) break;
879 0 : }
880 : }
881 0 : return status;
882 0 : }
883 :
884 : //______________________________________________
885 : Bool_t AliITSSAPTracker::RefitInward(int itr)
886 : {
887 : // refit track inward with material correction
888 0 : ITStrack_t &track = fTracks[itr];
889 0 : AliExternalTrackParam &trout = track.paramOut;
890 0 : if (trout.TestBit(kInvalidBit)) return kFALSE;
891 0 : AliExternalTrackParam &trin = track.paramInw;
892 0 : trin = trout;
893 : int ilA = kNLrActive;
894 0 : for (;ilA--;) { // find outermost layer with cluster
895 0 : if (track.clID[ilA]<0) continue;
896 : break;
897 : }
898 0 : int ilStart = fgkActiveLrITS[ilA]; // corresponding total lr id
899 0 : AliITSSAPLayer* lrA = fLayers[ilA];
900 0 : AliITSRecPoint *cl=lrA->GetClusterSorted(track.clID[ilA]);
901 0 : AliITSSAPLayer::ITSDetInfo_t& detInfo = lrA->GetDetInfo(cl->GetVolumeId()-lrA->GetVIDOffset());
902 0 : if (!trin.RotateParamOnly(detInfo.phiTF)) return kFALSE;
903 0 : if (!trin.PropagateParamOnlyTo(detInfo.xTF+cl->GetX(), fBz)) return kFALSE;
904 : // init with outer cluster y,z and slopes, q/pt of outward track
905 0 : double par[5] = {cl->GetY(), cl->GetZ(), trin.GetSnp(), trin.GetTgl(), trin.GetSigned1Pt()};
906 0 : double cov[15] = {cl->GetSigmaY2() + GetClSystYErr2(kALrSPD1),
907 0 : 0., cl->GetSigmaZ2() + GetClSystZErr2(kALrSPD1),
908 : 0,0,1,
909 : 0,0,0,1,
910 : 0,0,0,0,100*100};
911 0 : trin.Set(double(detInfo.xTF+cl->GetX()),double(detInfo.phiTF), par, cov);
912 : // !!! no material correction is needed: errors are not defined yer
913 : //
914 0 : for (int ilr=ilStart;ilr--;) {
915 : //
916 0 : if ( (ilA=fgkLr2Active[ilr])<0 || track.clID[ilA]<0) { // either passive layer or no cluster
917 0 : if (CrossPassiveLayer(trin,ilr)) continue;
918 0 : else return kFALSE;
919 : }
920 : // there is a cluster, need to update
921 0 : lrA = fLayers[ilA];
922 0 : cl = lrA->GetClusterSorted(track.clID[ilA]);
923 0 : AliITSSAPLayer::ITSDetInfo_t& detInfo1 = lrA->GetDetInfo(cl->GetVolumeId()-lrA->GetVIDOffset());
924 0 : if (!trin.Propagate(detInfo1.phiTF, detInfo1.xTF+cl->GetX(), fBz)) return kFALSE;
925 0 : double cpar[2]={ cl->GetY(), cl->GetZ()};
926 0 : double ccov[3]={ cl->GetSigmaY2() + GetClSystYErr2(ilA) , 0., cl->GetSigmaZ2() + GetClSystZErr2(ilA)};
927 0 : if (!trin.Update(cpar,ccov)) return kFALSE;
928 : //
929 : // correct for layer materials
930 0 : if (!trin.CorrectForMeanMaterial(fgkX2X0ITS[ilr],fgkRhoLITS[ilr],fgkDefMass,kFALSE)) return kFALSE;
931 : //
932 0 : }
933 : //
934 : // now go to PCA to vertex
935 : //double dca[2],dcaCov[3];
936 0 : if (!trin.PropagateToDCA(fSPDVertex,fBz,fgkRLayITS[kLrBeamPime])) return kFALSE; //,dca,dcaCov);
937 : //
938 0 : trin.ResetBit(kInvalidBit); // flag that inward fit succeeded
939 0 : return kTRUE;
940 : //
941 0 : }
942 :
943 : //______________________________________________
944 : void AliITSSAPTracker::RefitInward()
945 : {
946 : // refit tracks inward with material correction
947 0 : for (int itr=fNTracks;itr--;) {
948 0 : if (!RefitInward(itr)) {
949 : #ifdef _DEBUG_
950 : printf("RefitInward failed for track %d\n",itr);
951 : PrintTrack(fTracks[itr]);
952 : #endif
953 : }
954 : }
955 : //
956 0 : }
957 :
958 :
959 : //______________________________________________
960 : Bool_t AliITSSAPTracker::FitTrackVertex()
961 : {
962 : // Fit the vertexTracks. The inner tracks must be already propagated to the SPD vertex.
963 : // In this case straight line extrapolation can be used
964 : //
965 : const double kTiny = 1e-9;
966 : const double kTukey2 = 6;
967 0 : fTrackVertex.SetNContributors(0); // invalidate
968 0 : if (fNTracks<3) return kFALSE;
969 0 : fTrackVertex.SetXv(fSPDVertex->GetX());
970 0 : fTrackVertex.SetYv(fSPDVertex->GetY());
971 0 : fTrackVertex.SetZv(fSPDVertex->GetZ());
972 : //
973 0 : double vtxXYZ[3];
974 0 : fSPDVertex->GetXYZ(vtxXYZ); // initial vertex
975 : double scaleSigma2=9; // initial sigma scaling
976 0 : double dz[2],covdum[3],*covt;
977 : //
978 :
979 : #ifdef _DEBUG_
980 : AliRunLoader* rl = AliRunLoader::Instance();
981 : AliHeader* hd = 0;
982 : AliGenEventHeader* hdmc=0;
983 : TArrayF vtxMC(3);
984 : if (rl && (hd=rl->GetHeader()) && (hdmc=hd->GenEventHeader())) {
985 : hdmc->PrimaryVertex(vtxMC);
986 : }
987 : #endif
988 : //
989 : int nIter = 0;
990 0 : while(nIter++<fMaxVtxIter) {
991 : int ntAcc = 0;
992 : double wghSum=0,wghChi2=0;
993 : double cxx=0,cxy=0,cxz=0,cx0=0,cyy=0,cyz=0,cy0=0,czz=0,cz0=0;
994 : //
995 0 : for (int itr=fNTracks;itr--;) {
996 : //
997 0 : AliExternalTrackParam& trc = fTracks[itr].paramInw;
998 0 : if (trc.TestBit(kInvalidBit)) continue; // the track is invalidated, skip
999 0 : trc.ResetBit(kVtUsedBit);
1000 : //
1001 0 : double *param = (double*)trc.GetParameter();
1002 0 : double *covar = (double*)trc.GetCovariance();
1003 : //
1004 0 : double x0=trc.GetX();
1005 : double &y0=param[0];
1006 0 : double &z0=param[1];
1007 0 : double sn=param[2];
1008 0 : double cs2=(1.-sn)*(1.+sn);
1009 0 : if (cs2<kAlmost0) continue;
1010 0 : double cs=TMath::Sqrt(cs2), tgp=sn/cs, tgl=param[3]/cs;
1011 : // assume straight track equation Y=y0+tgp*X, Z=z0+tgl*X in tracking frame
1012 : //
1013 0 : double alp = trc.GetAlpha();
1014 0 : sn = TMath::Sin(alp); // parameters for rotation of vertex to
1015 0 : cs = TMath::Cos(alp); // tracking frame
1016 : //
1017 0 : double &syy=covar[0], &syz=covar[1], &szz=covar[2];
1018 0 : double detI = syy*szz - syz*syz;
1019 0 : if (TMath::Abs(detI)<kAlmost0) return kFALSE;
1020 0 : detI = 1./detI;
1021 0 : double syyI = szz*detI;
1022 0 : double szzI = syy*detI;
1023 0 : double syzI =-syz*detI;
1024 : //
1025 : // determine weight of the track
1026 0 : double vlocX = vtxXYZ[0]*cs+vtxXYZ[1]*sn;
1027 0 : double vlocY =-vtxXYZ[0]*sn+vtxXYZ[1]*cs;
1028 0 : double vlocZ = vtxXYZ[2];
1029 0 : double dy = y0 + tgp*(vlocX-x0) - vlocY;
1030 0 : double dz = z0 + tgl*(vlocX-x0) - vlocZ;
1031 0 : double chi2T = 0.5*(dy*dy*syyI + dz*dz*szzI) + dy*dz*syzI;
1032 0 : double wghT = (1-chi2T/kTukey2/scaleSigma2);
1033 0 : if (wghT<kTiny) continue;
1034 0 : wghSum += wghT;
1035 0 : wghChi2 += wghT*chi2T;
1036 : //
1037 0 : syyI *= wghT;
1038 0 : syzI *= wghT;
1039 0 : szzI *= wghT;
1040 : //
1041 0 : trc.SetBit(kVtUsedBit);
1042 : //
1043 : // printf("VTXFIT Bef %d X0= %+.4f Z= %+.4f Y=%+.4f\n",itr, x0, z0, y0);
1044 : //
1045 0 : double tmpSP = sn*tgp;
1046 0 : double tmpCP = cs*tgp;
1047 0 : double tmpSC = sn+tmpCP;
1048 0 : double tmpCS =-cs+tmpSP;
1049 0 : double tmpCL = cs*tgl;
1050 0 : double tmpSL = sn*tgl;
1051 0 : double tmpYXP = y0-tgp*x0;
1052 0 : double tmpZXL = z0-tgl*x0;
1053 : //
1054 0 : double tmpCLzz = tmpCL*szzI;
1055 0 : double tmpSLzz = tmpSL*szzI;
1056 0 : double tmpSCyz = tmpSC*syzI;
1057 0 : double tmpCSyz = tmpCS*syzI;
1058 0 : double tmpCSyy = tmpCS*syyI;
1059 0 : double tmpSCyy = tmpSC*syyI;
1060 0 : double tmpSLyz = tmpSL*syzI;
1061 0 : double tmpCLyz = tmpCL*syzI;
1062 : //
1063 0 : cxx += tmpCL*(tmpCLzz+tmpSCyz+tmpSCyz)+tmpSC*tmpSCyy; // dchi^2/dx/dx
1064 0 : cxy += tmpCL*(tmpSLzz+tmpCSyz)+tmpSL*tmpSCyz+tmpSC*tmpCSyy; // dchi^2/dx/dy
1065 0 : cxz += -sn*syzI-tmpCLzz-tmpCP*syzI; // dchi^2/dx/dz
1066 0 : cx0 += -(tmpCLyz+tmpSCyy)*tmpYXP-(tmpCLzz+tmpSCyz)*tmpZXL; // RHS
1067 : //
1068 : //double cyx
1069 0 : cyy += tmpSL*(tmpSLzz+tmpCSyz+tmpCSyz)+tmpCS*tmpCSyy; // dchi^2/dy/dy
1070 0 : cyz += -(tmpCSyz+tmpSLzz); // dchi^2/dy/dz
1071 0 : cy0 += -tmpYXP*(tmpCSyy+tmpSLyz)-tmpZXL*(tmpCSyz+tmpSLzz); // RHS
1072 : //
1073 : //double czx
1074 : //double czy
1075 0 : czz += szzI; // dchi^2/dz/dz
1076 0 : cz0 += tmpZXL*szzI+tmpYXP*syzI; // RHS
1077 : //
1078 0 : ntAcc++;
1079 0 : }
1080 : //
1081 0 : if (ntAcc<2) break; // failed
1082 : //
1083 0 : double vec[3] = {cx0,cy0,cz0};
1084 0 : AliSymMatrix mat(3);
1085 0 : mat(0,0) = cxx;
1086 0 : mat(0,1) = cxy;
1087 0 : mat(0,2) = cxz;
1088 0 : mat(1,1) = cyy;
1089 0 : mat(1,2) = cyz;
1090 0 : mat(2,2) = czz;
1091 : //
1092 : #ifdef _DEBUG_
1093 : printf("MatBefore: \n"); mat.Print("d");
1094 : #endif
1095 0 : if (!mat.SolveChol(vec,kTRUE)) return kFALSE;
1096 : #ifdef _DEBUG_
1097 : printf("MatAfter : \n"); mat.Print("d");
1098 : #endif
1099 : //
1100 0 : double scaleSigma2New = wghChi2/wghSum;
1101 : //
1102 : #ifdef _DEBUG_
1103 : double dVtX = vec[0] - vtxXYZ[0];
1104 : double dVtY = vec[1] - vtxXYZ[1];
1105 : double dst2 = dVtX*dVtX+dVtY*dVtY;
1106 : double dVtZ = vec[2] - vtxXYZ[2];
1107 : printf("VTIter%d %d %d %+e %+e %e %+e %.3f %.3f %e %e %e %e %e\n",
1108 : nIter,ntAcc,fNTracks,dVtX,dVtY,dVtZ,dst2,
1109 : scaleSigma2,scaleSigma2New, wghChi2,wghSum,
1110 : vec[0]-vtxMC[0],
1111 : vec[1]-vtxMC[1],
1112 : vec[2]-vtxMC[2]
1113 : );
1114 : #endif
1115 : //
1116 0 : double vtCov[6] = {mat(0,0),mat(0,1),mat(1,1),mat(0,2),mat(1,2),mat(2,2)};
1117 0 : fTrackVertex.SetXYZ(vec);
1118 0 : fTrackVertex.SetCovarianceMatrix(vtCov);
1119 0 : fTrackVertex.SetNContributors(ntAcc);
1120 : //
1121 0 : if (scaleSigma2<1. &&
1122 0 : scaleSigma2New/scaleSigma2>fStopScaleChange) break;
1123 : scaleSigma2 = scaleSigma2New;
1124 0 : for (int i=3;i--;) vtxXYZ[i] = vec[i];
1125 : //
1126 0 : }
1127 : // calculate explicitly chi2
1128 : double chiTRC = 0;
1129 : double chiSPD = 0;
1130 : //
1131 0 : for (int itr=fNTracks;itr--;) {
1132 0 : AliExternalTrackParam& trc = fTracks[itr].paramInw;
1133 0 : if (trc.TestBit(kInvalidBit)) continue; // the track is invalidated, skip
1134 0 : AliExternalTrackParam trT(trc);
1135 0 : AliExternalTrackParam trS(trc);
1136 0 : trT.PropagateToDCA(&fTrackVertex,fBz,10,dz,covdum);
1137 0 : covt = (double*)trT.GetCovariance();
1138 0 : double detI = covt[0]*covt[2] - covt[1]*covt[1];
1139 0 : detI = 1./detI;
1140 0 : double syyI = covt[2]*detI;
1141 0 : double szzI = covt[0]*detI;
1142 0 : double syzI =-covt[1]*detI;
1143 0 : chiTRC += dz[0]*dz[0]*syyI + dz[1]*dz[1]*szzI + 2*dz[0]*dz[1]*syzI;
1144 : //
1145 0 : trS.PropagateToDCA(fSPDVertex,fBz,10,dz,covdum);
1146 0 : covt = (double*)trT.GetCovariance();
1147 0 : detI = covt[0]*covt[2] - covt[1]*covt[1];
1148 0 : detI = 1./detI;
1149 0 : syyI = covt[2]*detI;
1150 0 : szzI = covt[0]*detI;
1151 0 : syzI =-covt[1]*detI;
1152 0 : chiSPD += dz[0]*dz[0]*syyI + dz[1]*dz[1]*szzI + 2*dz[0]*dz[1]*syzI;
1153 : // printf("VTXFITChi2 Aft %d X0= %+.4f Z= %+.4f Y=:%+.4f SPD: X:%+.4f Z:%+.4f Y:%+.4f\n",itr,
1154 : // trT.GetX(), trT.GetZ(), trT.GetY(),
1155 : // trS.GetX(), trS.GetZ(), trS.GetY());
1156 0 : }
1157 : #ifdef _DEBUG_
1158 : /*
1159 : //-------------------------TMP>>>
1160 : AliRunLoader* rl = AliRunLoader::Instance();
1161 : AliHeader* hd = 0;
1162 : AliGenEventHeader* hdmc=0;
1163 : TArrayF vtxMC(3);
1164 : if (rl && (hd=rl->GetHeader()) && (hdmc=hd->GenEventHeader())) {
1165 : hdmc->PrimaryVertex(vtxMC);
1166 : }
1167 : printf("VTFIT %f %f %f %d %8.2f %8.2f %.4f %.4f %.4f %.4f %.4f %.4f\n",
1168 : vtxMC[0],vtxMC[1],vtxMC[2],
1169 : ntAcc,chiTRC,chiSPD,
1170 : fTrackVertex.GetX(),fTrackVertex.GetY(),fTrackVertex.GetZ(),
1171 : fSPDVertex->GetX(),fSPDVertex->GetY(),fSPDVertex->GetZ());
1172 : //-------------------------TMP<<<
1173 : */
1174 : #endif
1175 : //
1176 : return kTRUE;
1177 0 : }
1178 :
1179 : #ifdef _CONTROLH_
1180 : //______________________________________________
1181 : void AliITSSAPTracker::FillRecoStat()
1182 : {
1183 : // fill data for preformance study
1184 : //
1185 : AliStack* stack = 0;
1186 : AliRunLoader* rl = AliRunLoader::Instance();
1187 : if (!rl || !(stack=rl->Stack())) return;
1188 : //
1189 : TBits patternMC;
1190 : enum {kIsPrim=kNLrActive,kValidTracklet,kValidTrack,kRecDone,kBitPerTrack};
1191 : int nTrkMC = stack->GetNtrack();
1192 : patternMC.SetBitNumber(nTrkMC*kBitPerTrack,0);
1193 : //
1194 : // fill MC track patterns
1195 : for (int ilr=kNLrActive;ilr--;) {
1196 : AliITSSAPLayer *lr = fLayers[ilr];
1197 : int ncl = lr->GetNClusters();
1198 : for (int icl=ncl;icl--;) {
1199 : AliITSRecPoint* cl = lr->GetClusterUnSorted(icl);
1200 : for (int j=0;j<3;j++) {
1201 : int lb = cl->GetLabel(j);
1202 : if (lb<0 || lb>=nTrkMC) break;
1203 : patternMC.SetBitNumber(lb*kBitPerTrack+ilr,kTRUE);
1204 : }
1205 : }
1206 : }
1207 : // set reconstructability
1208 : for (int itr=nTrkMC;itr--;) {
1209 : int bitoffs = itr*kBitPerTrack;
1210 : Bool_t isPrim = stack->IsPhysicalPrimary(itr);
1211 : patternMC.SetBitNumber(bitoffs+kIsPrim,isPrim);
1212 : if (patternMC.TestBitNumber(bitoffs+kALrSPD1) && patternMC.TestBitNumber(bitoffs+kALrSPD2)) {
1213 : patternMC.SetBitNumber(bitoffs+kValidTracklet,kTRUE);
1214 : //
1215 : TParticle* mctr = stack->Particle(itr);
1216 : fHTrackletMC->Fill(mctr->Pt(),isPrim);
1217 : // check outer layers reconstructability
1218 : int nmiss = 0;
1219 : for (int il=kALrSDD1;il<=kALrSSD2;il++) if (!fSkipLayer[il] && !patternMC.TestBitNumber(bitoffs+il)) nmiss++;
1220 : if (nmiss<=fMaxMissedLayers) {
1221 : patternMC.SetBitNumber(bitoffs+kValidTrack);
1222 : fHTrackMC->Fill(mctr->Pt(),isPrim);
1223 : }
1224 : }
1225 : }
1226 : //
1227 : int nTrk = GetNTracklets();
1228 : if (!nTrk) return;
1229 : for (int itr=0;itr<nTrk;itr++) {
1230 : SPDtracklet_t& trlet = fTracklets[itr];
1231 : // PrintTracklet(itr);
1232 : //
1233 : int lbl = trlet.label;
1234 : if (lbl==fgkDummyLabel) continue;
1235 : int lblA = TMath::Abs(lbl);
1236 : int bitoffs = lblA*kBitPerTrack;
1237 : Bool_t isPrim = patternMC.TestBitNumber(bitoffs+kIsPrim);
1238 : TParticle* mctr = stack->Particle(lblA);
1239 : double pt = mctr->Pt();
1240 : fHTrackletAll->Fill(pt,isPrim);
1241 : if (lbl<0) fHTrackletFake->Fill(pt,isPrim);
1242 : }
1243 : //
1244 : nTrk = GetNTracks();
1245 : for (int itr=0;itr<nTrk;itr++) {
1246 : ITStrack_t &track = fTracks[itr];
1247 : CookLabel(track);
1248 : //
1249 : int lbl = track.label;
1250 : if (lbl==fgkDummyLabel) continue;
1251 : int lblA = TMath::Abs(lbl);
1252 : int bitoffs = lblA*kBitPerTrack;
1253 : Bool_t isPrim = patternMC.TestBitNumber(bitoffs+kIsPrim);
1254 : TParticle* mctr = stack->Particle(lblA);
1255 : double pt = mctr->Pt();
1256 : Bool_t clone = patternMC.TestBitNumber(bitoffs+kRecDone); // was the track already reconstructed?
1257 : float bn = isPrim ? (clone ? 2:1):(clone ? -1:0); // fill clones in over/underflow
1258 : fHTrackAll->Fill(pt,bn);
1259 : patternMC.SetBitNumber(bitoffs+kRecDone);
1260 : if (lbl<0) fHTrackFake->Fill(pt,bn);
1261 : }
1262 : //
1263 : AliHeader* hd = rl->GetHeader();
1264 : AliGenEventHeader* hdmc;
1265 : TArrayF vtxMC;
1266 : if (hd && (hdmc=hd->GenEventHeader())) hdmc->PrimaryVertex(vtxMC);
1267 : //
1268 : nTrk = GetNTracklets();
1269 : fHVtxMltRef->Fill(nTrk);
1270 : if (fTrackVertex.GetStatus()==1) {
1271 : if (hdmc) {
1272 : double dx = vtxMC[0]-fTrackVertex.GetX();
1273 : double dy = vtxMC[1]-fTrackVertex.GetY();
1274 : double dz = vtxMC[2]-fTrackVertex.GetZ();
1275 : fHVtxDiffXY->Fill(dx,dy);
1276 : fHVtxDiffZ->Fill(dz);
1277 : fHVtxDiffXMlt->Fill(nTrk, dx);
1278 : fHVtxDiffYMlt->Fill(nTrk, dy);
1279 : fHVtxDiffZMlt->Fill(nTrk, dz);
1280 : //
1281 : double sig[3];
1282 : fTrackVertex.GetSigmaXYZ(sig);
1283 : if (sig[0]>0) fHVtxPullXMlt->Fill(nTrk, dx/sig[0]);
1284 : if (sig[1]>0) fHVtxPullYMlt->Fill(nTrk, dy/sig[1]);
1285 : if (sig[2]>0) fHVtxPullZMlt->Fill(nTrk, dz/sig[2]);
1286 : }
1287 : fHVtxOKMlt->Fill(nTrk);
1288 : }
1289 : //
1290 : if (fSPDVertex->GetStatus()==1 && hdmc) {
1291 : double dx = vtxMC[0]-fSPDVertex->GetX();
1292 : double dy = vtxMC[1]-fSPDVertex->GetY();
1293 : double dz = vtxMC[2]-fSPDVertex->GetZ();
1294 : fHVtxMCSPDDiffXY->Fill(dx,dy);
1295 : fHVtxMCSPDDiffZ->Fill(dz);
1296 : fHVtxMCSPDDiffXMlt->Fill(nTrk, dx);
1297 : fHVtxMCSPDDiffYMlt->Fill(nTrk, dy);
1298 : fHVtxMCSPDDiffZMlt->Fill(nTrk, dz);
1299 : //
1300 : double sig[3];
1301 : fSPDVertex->GetSigmaXYZ(sig);
1302 : if (sig[0]>0) fHVtxMCSPDPullXMlt->Fill(nTrk, dx/sig[0]);
1303 : if (sig[1]>0) fHVtxMCSPDPullYMlt->Fill(nTrk, dy/sig[1]);
1304 : if (sig[2]>0) fHVtxMCSPDPullZMlt->Fill(nTrk, dz/sig[2]);
1305 : //
1306 : }
1307 : //
1308 : }
1309 :
1310 : //______________________________________________
1311 : void AliITSSAPTracker::BookHistos()
1312 : {
1313 : // book control histos
1314 : const int kNBinMlt=20, kNBPt=15, kNBDiffVtx=50, kNResBins=250,kNPullBins=50,kNChiClBins=50,kNBPullVtx=50;
1315 : const double kMinMlt=1,kMaxMlt=5000,kMinPt=0.01,kMaxPt=3, kMaxDiffVtx=0.05, kMaxResidYZ=2.5,kMaxPullYZ=10,kChiClMax=100,kMaxPullVtx=10;
1316 : //
1317 : double* axLogPt = DefLogAx(kMinPt,kMaxPt,kNBPt);
1318 : double* axLogMlt = DefLogAx(kMinMlt,kMaxMlt,kNBinMlt);
1319 :
1320 : for (int ilr=0;ilr<kNLrActive;ilr++) {
1321 : //
1322 : // ----------------- These are histos to be filled during tracking
1323 : // PropagateBack and RefitInward will be stored among the histos of 1st pass
1324 : //
1325 : int ilrS = ilr*10;
1326 : TString ttl = Form("residY%d",ilr);
1327 : TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNResBins,-kMaxResidYZ,kMaxResidYZ);
1328 : fArrHisto.AddAtAndExpand(hdy,ilrS+kHResidY);
1329 : hdy->SetDirectory(0);
1330 : //
1331 : ttl = Form("residYPull%d",ilr);
1332 : TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNPullBins,-kMaxPullYZ,kMaxPullYZ);
1333 : fArrHisto.AddAtAndExpand(hdyp,ilrS+kHPullY);
1334 : hdyp->SetDirectory(0);
1335 : //
1336 : ttl = Form("residZ%d",ilr);
1337 : TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNResBins,-kMaxResidYZ,kMaxResidYZ);
1338 : fArrHisto.AddAtAndExpand(hdz,ilrS+kHResidZ);
1339 : hdz->SetDirectory(0);
1340 : //
1341 : ttl = Form("residZPull%d",ilr);
1342 : TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNPullBins,-kMaxPullYZ,kMaxPullYZ);
1343 : hdzp->SetDirectory(0);
1344 : fArrHisto.AddAtAndExpand(hdzp,ilrS+kHPullZ);
1345 : //
1346 : ttl = Form("chi2Cl%d",ilr);
1347 : TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt, kNChiClBins,0.,kChiClMax);
1348 : hchi->SetDirectory(0);
1349 : fArrHisto.AddAtAndExpand(hchi,ilrS+kHChi2Cl);
1350 : } // loop over layers
1351 : //
1352 : fHChi2NDFvsPT = new TH2F("chi2ndfPT","chi2/ndf total vs pt",kNBPt,axLogPt, kNChiClBins,0.,kChiClMax);
1353 : fArrHisto.AddLast(fHChi2NDFvsPT);
1354 : fHChi2NDFvsPT->SetDirectory(0);
1355 : //
1356 : fHChi2vsNC = new TH2F("chi2NC","chi2 total vs NCl",kNLrActive-2,2.5,kNLrActive+0.5, kNChiClBins,0.,kChiClMax);
1357 : fArrHisto.AddLast(fHChi2vsNC);
1358 : fHChi2vsNC->SetDirectory(0);
1359 :
1360 : // SPDvertex vs MC
1361 : fHVtxMCSPDDiffXY = new TH2F("vtxMCSPDDiffXY","vtxMC-vtxSPD XY",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx,
1362 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1363 : fArrHisto.AddLast(fHVtxMCSPDDiffXY);
1364 : fHVtxMCSPDDiffXY->SetDirectory(0);
1365 : //
1366 : fHVtxMCSPDDiffZ = new TH1F("vtxMCSPDDiffZ","vtxMC-vtxSPD Z",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1367 : fArrHisto.AddLast(fHVtxMCSPDDiffZ);
1368 : fHVtxMCSPDDiffZ->SetDirectory(0);
1369 : //
1370 : fHVtxMCSPDDiffXMlt = new TH2F("VtxMCSPDDiffXMlt","vX_{MC}-vX_{SPD} vs mlt",kNBinMlt,axLogMlt,
1371 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1372 : fArrHisto.AddLast(fHVtxMCSPDDiffXMlt);
1373 : fHVtxMCSPDDiffXMlt->SetDirectory(0);
1374 : //
1375 : fHVtxMCSPDDiffYMlt = new TH2F("VtxMCSPDDiffYMlt","vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1376 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1377 : fArrHisto.AddLast(fHVtxMCSPDDiffYMlt);
1378 : fHVtxMCSPDDiffYMlt->SetDirectory(0);
1379 : //
1380 : fHVtxMCSPDDiffZMlt = new TH2F("VtxMCSPDDiffZMlt","vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1381 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1382 : fArrHisto.AddLast(fHVtxMCSPDDiffZMlt);
1383 : fHVtxMCSPDDiffZMlt->SetDirectory(0);
1384 : //
1385 : //
1386 : fHVtxMCSPDPullXMlt = new TH2F("VtxMCSPDPullXMlt","Pull vX_{MC}-vX_{SPD} vs mlt",kNBinMlt,axLogMlt,
1387 : kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1388 : fArrHisto.AddLast(fHVtxMCSPDPullXMlt);
1389 : fHVtxMCSPDPullXMlt->SetDirectory(0);
1390 : //
1391 : fHVtxMCSPDPullYMlt = new TH2F("VtxMCSPDPullYMlt","Pull vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1392 : kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1393 : fArrHisto.AddLast(fHVtxMCSPDPullYMlt);
1394 : fHVtxMCSPDPullYMlt->SetDirectory(0);
1395 : //
1396 : fHVtxMCSPDPullZMlt = new TH2F("VtxMCSPDPullZMlt","Pull vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1397 : kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1398 : fArrHisto.AddLast(fHVtxMCSPDPullZMlt);
1399 : fHVtxMCSPDPullZMlt->SetDirectory(0);
1400 : //
1401 : fHTrackletMC = new TH2F("MCRefTracklet","MCRef Tracklet",kNBPt,axLogPt, 2, -0.5, 1.5);
1402 : fHTrackletMC->SetXTitle("p_{T}");
1403 : fHTrackletMC->GetYaxis()->SetBinLabel(1,"sec");
1404 : fHTrackletMC->GetYaxis()->SetBinLabel(2,"prim");
1405 : fArrHisto.AddLast(fHTrackletMC);
1406 : fHTrackletMC->SetDirectory(0);
1407 : //
1408 : fHTrackletAll = new TH2F("TrackletAll","Tracklet All rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1409 : fHTrackletAll->SetXTitle("p_{T}");
1410 : fHTrackletAll->GetYaxis()->SetBinLabel(1,"sec");
1411 : fHTrackletAll->GetYaxis()->SetBinLabel(2,"prim");
1412 : fArrHisto.AddLast(fHTrackletAll);
1413 : fHTrackletAll->SetDirectory(0);
1414 : //
1415 : fHTrackletFake = new TH2F("TrackletFake","Tracklet Fake rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1416 : fHTrackletFake->SetXTitle("p_{T}");
1417 : fHTrackletFake->GetYaxis()->SetBinLabel(1,"sec");
1418 : fHTrackletFake->GetYaxis()->SetBinLabel(2,"prim");
1419 : fArrHisto.AddLast(fHTrackletFake);
1420 : fHTrackletFake->SetDirectory(0);
1421 : //
1422 : fHTrackMC = new TH2F("MCRefTrack","MCRef Track",kNBPt,axLogPt, 2, -0.5, 1.5);
1423 : fHTrackMC->SetXTitle("p_{T}");
1424 : fHTrackMC->GetYaxis()->SetBinLabel(1,"sec");
1425 : fHTrackMC->GetYaxis()->SetBinLabel(2,"prim");
1426 : fArrHisto.AddLast(fHTrackMC);
1427 : fHTrackMC->SetDirectory(0);
1428 : //
1429 : fHTrackAll = new TH2F("TrackAll","Track All rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1430 : fHTrackAll->SetXTitle("p_{T}");
1431 : fHTrackAll->GetYaxis()->SetBinLabel(1,"sec");
1432 : fHTrackAll->GetYaxis()->SetBinLabel(2,"prim");
1433 : fArrHisto.AddLast(fHTrackAll);
1434 : fHTrackAll->SetDirectory(0);
1435 : //
1436 : fHTrackFake = new TH2F("TrackFake","Track Fake rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1437 : fHTrackFake->SetXTitle("p_{T}");
1438 : fHTrackFake->GetYaxis()->SetBinLabel(1,"sec");
1439 : fHTrackFake->GetYaxis()->SetBinLabel(2,"prim");
1440 : fArrHisto.AddLast(fHTrackFake);
1441 : fHTrackFake->SetDirectory(0);
1442 : //
1443 : fHVtxMltRef = new TH1F("vtxRef","vtxRef",kNBinMlt,axLogMlt);
1444 : fHVtxMltRef->SetXTitle("sqrt(Ntracklets)");
1445 : fArrHisto.AddLast(fHVtxMltRef);
1446 : fHVtxMltRef->SetDirectory(0);
1447 : //
1448 : fHVtxOKMlt = new TH1F("vtxOK","vtxOK",kNBinMlt,axLogMlt);
1449 : fHVtxOKMlt->SetXTitle("sqrt(Ntracklets)");
1450 : fArrHisto.AddLast(fHVtxOKMlt);
1451 : fHVtxOKMlt->SetDirectory(0);
1452 : //
1453 : fHVtxDiffXY = new TH2F("vtxMCDiffXY","vtxMC-vtxRec XY",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx,
1454 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1455 : fArrHisto.AddLast(fHVtxDiffXY);
1456 : fHVtxDiffXY->SetDirectory(0);
1457 : //
1458 : fHVtxDiffZ = new TH1F("vtxMCDiffZ","vtxMC-vtxRec Z",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1459 : fArrHisto.AddLast(fHVtxDiffZ);
1460 : fHVtxDiffZ->SetDirectory(0);
1461 : //
1462 : fHVtxDiffXMlt = new TH2F("VtxDiffXMlt","vX_{MC}-vX_{rec} vs mlt",kNBinMlt,axLogMlt,
1463 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1464 : fArrHisto.AddLast(fHVtxDiffXMlt);
1465 : fHVtxDiffXMlt->SetDirectory(0);
1466 : //
1467 : fHVtxDiffYMlt = new TH2F("VtxDiffYMlt","vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1468 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1469 : fArrHisto.AddLast(fHVtxDiffYMlt);
1470 : fHVtxDiffYMlt->SetDirectory(0);
1471 : //
1472 : fHVtxDiffZMlt = new TH2F("VtxDiffZMlt","vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1473 : kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1474 : fArrHisto.AddLast(fHVtxDiffZMlt);
1475 : fHVtxDiffZMlt->SetDirectory(0);
1476 : //
1477 : fHVtxPullXMlt = new TH2F("VtxPullXMlt","Pull vX_{MC}-vX_{rec} vs mlt",kNBinMlt,axLogMlt,
1478 : kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1479 : fArrHisto.AddLast(fHVtxPullXMlt);
1480 : fHVtxPullXMlt->SetDirectory(0);
1481 : //
1482 : fHVtxPullYMlt = new TH2F("VtxPullYMlt","Pull vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1483 : kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1484 : fArrHisto.AddLast(fHVtxPullYMlt);
1485 : fHVtxPullYMlt->SetDirectory(0);
1486 : //
1487 : fHVtxPullZMlt = new TH2F("VtxPullZMlt","Pull vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1488 : kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1489 : fArrHisto.AddLast(fHVtxPullZMlt);
1490 : fHVtxPullZMlt->SetDirectory(0);
1491 : //
1492 : }
1493 :
1494 : //______________________________________________
1495 : void AliITSSAPTracker::SaveHistos(const char* outFName)
1496 : {
1497 : // save control histos
1498 : TString fnms = outFName;
1499 : if (fnms.IsNull()) fnms = "XXXITSTrackerControlH.root";
1500 : TFile* fl = TFile::Open(fnms.Data(),"recreate");
1501 : if (!fl) {
1502 : printf("Failed to open output file %s\n",fnms.Data());
1503 : return;
1504 : }
1505 : fArrHisto.Write();
1506 : fl->Close();
1507 : delete fl;
1508 : printf("Stored control histos in %s\n",fnms.Data());
1509 : //
1510 : }
1511 :
1512 : //______________________________________________
1513 : void AliITSSAPTracker::FillTrackingControlHistos(int lrID,int lbl,const AliExternalTrackParam* track,
1514 : const double cpar[2],const double ccov[3],
1515 : const AliITSRecPoint* bestCl)
1516 : {
1517 : // fill control histos for tracking for correct matches
1518 : Bool_t corr = kFALSE;
1519 : for (int i=0;i<3;i++) if (bestCl->GetLabel(i)==lbl) {corr=kTRUE; break;}
1520 : if (!corr) return;
1521 : double pt = track->Pt();
1522 : double dy = cpar[0]-track->GetY();
1523 : double dz = cpar[1]-track->GetZ();
1524 : double sgy = TMath::Sqrt(ccov[0]+track->GetSigmaY2());
1525 : double sgz = TMath::Sqrt(ccov[2]+track->GetSigmaZ2());
1526 : int lrIDS = lrID*10;
1527 : ((TH2F*)fArrHisto[lrIDS+kHResidY])->Fill(pt,dy);
1528 : ((TH2F*)fArrHisto[lrIDS+kHPullY])->Fill(pt,dy/sgy);
1529 : ((TH2F*)fArrHisto[lrIDS+kHResidZ])->Fill(pt,dz);
1530 : ((TH2F*)fArrHisto[lrIDS+kHPullZ])->Fill(pt,dz/sgz);
1531 : ((TH2F*)fArrHisto[lrIDS+kHChi2Cl])->Fill(pt,track->GetPredictedChi2(cpar,ccov));
1532 : //
1533 : }
1534 :
1535 : //______________________________________________
1536 : Double_t* AliITSSAPTracker::DefLogAx(double xMn,double xMx, int nbin)
1537 : {
1538 : // get array for log axis
1539 : if (xMn<=0 || xMx<=xMn || nbin<2) {
1540 : printf("Wrong axis request: %f %f %d\n",xMn,xMx,nbin);
1541 : return 0;
1542 : }
1543 : double dx = log(xMx/xMn)/nbin;
1544 : double *xax = new Double_t[nbin+1];
1545 : for (int i=0;i<=nbin;i++) xax[i]= xMn*exp(dx*i);
1546 : return xax;
1547 : }
1548 :
1549 : #endif
1550 : //
1551 : #ifdef _TIMING_
1552 : //______________________________________________
1553 : void AliITSSAPTracker::PrintTiming()
1554 : {
1555 : // print timing info
1556 : for (int i=0;i<kNSW;i++) {printf("%-10s:\t",fgkSWNames[i]); fSW[i].Print();}
1557 : }
1558 : #endif
|
