]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/TPCLib/AliHLTTPCCATracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
effC++ and warnings
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCCATracker.cxx
CommitLineData
21b6a334 1/** \class AliHLTTPCCATracker
2<pre>
3//_____________________________________________________________
4// AliHLTTPCCATracker
5//
6//
7// Author: Ivan Kisel
8// Copyright &copy ALICE HLT Group
9</pre>
10*/
11
12#define WRITETRACKS 1
13
14#include <sys/time.h>
15
16
17#include "AliHLTTPCRootTypes.h"
18#include "AliHLTTPCSpacePointData.h"
19#include "AliHLTTPCLogging.h"
20#include "AliHLTTPCVertex.h"
21#include "AliHLTTPCConfMapTrack.h"
22#include "AliHLTTPCConfMapPoint.h"
23#include "AliHLTTPCTrackArray.h"
24#include "AliHLTTPCTransform.h"
25#include "AliHLTTPCCATracker.h"
26#include "AliHLTTPCTrackSegmentData.h"
27
28#include "TApplication.h"
29#include "TStyle.h"
30#include "TCanvas.h"
31#include "TFrame.h"
32#include "TMarker.h"
33#include "TPaveText.h"
34#include "TEllipse.h"
35#include "TText.h"
44815588 36#include "TROOT.h"
21b6a334 37
38#include <iostream>
39
40#if __GNUC__ >= 3
41using namespace std;
42#endif
43
44static TList *listHisto;
45static TH1F *h_NClusters;
46static TH2F *h_ClustersXY[10];
47
48static Int_t h_Event;
49
50//#define DRAW
51
52#ifdef DRAW
53
54static TApplication *myapp;
55static TCanvas *YX, *ZX;
56
57static bool ask = false;
58
59#endif //DRAW
60
61ClassImp(AliHLTTPCCATracker)
62
63// ----------------------------------------------------------------------------------
64AliHLTTPCCATracker::AliHLTTPCCATracker()
65{
66 //Default constructor
67 fVertex = NULL;
68 fTrack = NULL;
69 fHit = NULL;
70 fVolume = NULL;
71 fRow = NULL;
72 fBench = (Bool_t)true;
73 fVertexConstraint = (Bool_t)true;
74 fParamSet[0]=0;
75 fParamSet[1]=0;
76 //JMT 2006/11/13
77 fOutputPtr= NULL;
78 fOutputNTracks=0;
79 fOutputSize=0;
80}
81
82// ----------------------------------------------------------------------------------
83AliHLTTPCCATracker::~AliHLTTPCCATracker()
84{
85 // Destructor.
86 if(fVolume) {
87 delete [] fVolume;
88 }
89 if(fRow) {
90 delete [] fRow;
91 }
92 if(fHit) {
93 delete [] fHit;
94 }
95 if(fTrack) {
96 delete fTrack;
97 }
98 //JMT 2006/11/13
99 if(fOutputPtr){
100 fOutputPtr = NULL;
101 delete fOutputPtr;
102 }
103
104}
105
106// ----------------------------------------------------------------------------------
107void AliHLTTPCCATracker::CACreateHistos()
108{
109
110 // IK histogramming
111 static bool first_call = true;
112 if (first_call){
113 first_call = false;
114
115 h_Event = 0;
116
117 TDirectory *curdir = gDirectory;
118 TDirectory *histodir = gROOT->mkdir("histodir");
119 histodir->cd();
120
121 h_NClusters = new TH1F("h_NClusters", "Number of clusters in event", 100, 0.0, 1000.0);
122
123 h_ClustersXY[0] = new TH2F("h_ClustersXY_0", "Clusters in XY plane Ev 0", 100, 50.0, 300.0, 100, -100.0, 100.0);
124 h_ClustersXY[1] = new TH2F("h_ClustersXY_1", "Clusters in XY plane Ev 1", 100, 50.0, 300.0, 100, -100.0, 100.0);
125 h_ClustersXY[2] = new TH2F("h_ClustersXY_2", "Clusters in XY plane Ev 2", 100, 50.0, 300.0, 100, -100.0, 100.0);
126 h_ClustersXY[3] = new TH2F("h_ClustersXY_3", "Clusters in XY plane Ev 3", 100, 50.0, 300.0, 100, -100.0, 100.0);
127 h_ClustersXY[4] = new TH2F("h_ClustersXY_4", "Clusters in XY plane Ev 4", 100, 50.0, 300.0, 100, -100.0, 100.0);
128 h_ClustersXY[5] = new TH2F("h_ClustersXY_5", "Clusters in XY plane Ev 5", 100, 50.0, 300.0, 100, -100.0, 100.0);
129 h_ClustersXY[6] = new TH2F("h_ClustersXY_6", "Clusters in XY plane Ev 6", 100, 50.0, 300.0, 100, -100.0, 100.0);
130 h_ClustersXY[7] = new TH2F("h_ClustersXY_7", "Clusters in XY plane Ev 7", 100, 50.0, 300.0, 100, -100.0, 100.0);
131 h_ClustersXY[8] = new TH2F("h_ClustersXY_8", "Clusters in XY plane Ev 8", 100, 50.0, 300.0, 100, -100.0, 100.0);
132 h_ClustersXY[9] = new TH2F("h_ClustersXY_9", "Clusters in XY plane Ev 9", 100, 50.0, 300.0, 100, -100.0, 100.0);
133
134 listHisto = gDirectory->GetList();
135 curdir->cd();
136
137#ifdef DRAW
138
139 myapp = new TApplication("myapp", 0, 0);
140
141 gStyle->SetCanvasBorderMode(0);
142 gStyle->SetCanvasBorderSize(1);
143 gStyle->SetCanvasColor(0);
144
145
146 YX = new TCanvas ("YX", "YX (Pad-Row) window", -10, 0, 680, 400);
147 YX->Range(-50.0, 80.0, 50.0, 250.0);
148 YX->Clear();
149 YX->Draw();
150
151 ZX = new TCanvas ("ZX", "ZX window", -700, 0, 680, 400);
152 ZX->Range(-250.0, 80.0, 250.0, 250.0);
153 ZX->Clear();
154 ZX->Draw();
155
156#ifdef XXX
157 char symbol;
158 if (ask){
159 do{
160 std::cin.get(symbol);
161 if (symbol == 'r')
162 ask = false;
163 } while (symbol != '\n');
164 }
165#endif //XXX
166
167#endif //DRAW
168
169 }
170}
171
172// ----------------------------------------------------------------------------------
173void AliHLTTPCCATracker::CAWriteHistos()
174{
175 // IK Open the output file and write histogramms
176
177 TFile *outfile = new TFile("ConfMapper_histo.root", "RECREATE");
178 TIter hiter(listHisto);
179 while (TObject *obj = hiter()) obj->Write();
180 outfile->Close();
181}
182
183// ----------------------------------------------------------------------------------
184Bool_t AliHLTTPCCATracker::ReadHits(UInt_t count, AliHLTTPCSpacePointData* hits )
185{
186
187 //read hits
188 Int_t nhit=(Int_t)count;
189
190 vector<CAHit> vec_hits;
191 vec_hits.clear();
192
193 for (Int_t i=0;i<nhit;i++)
194 {
195 fHit[i+fClustersUnused].Reset();
196 fHit[i+fClustersUnused].ReadHits(&(hits[i]));
197 }
198
199 fClustersUnused += nhit;
200
201 return true;
202}
203
204// ----------------------------------------------------------------------------------
205void AliHLTTPCCATracker::CAInitialize()
206{
207 vec_hits.clear();
208 vec_patch_tracks.clear();
209 vec_slice_tracks.clear();
210
211 return;
212}
213
214// ----------------------------------------------------------------------------------
215void AliHLTTPCCATracker::CAReadPatchHits(Int_t patch, UInt_t count, AliHLTTPCSpacePointData* hits )
216{
217 //read hits
218 Int_t nhit=(Int_t)count;
219
220 patch_first_hit_ind = vec_hits.size();
221 patch_last_hit_ind = patch_first_hit_ind + nhit - 1;
222
223 for (Int_t i=0;i<nhit;i++)
224 {
225 AliHLTTPCSpacePointData* pSP = &(hits[i]);
226
227 CAHit hit;
228 hit.x = pSP->fX;
229 hit.y = pSP->fY;
230 hit.z = pSP->fZ;
231
232 hit.errx = 1.0;//pSP->fSigmaY2;
233 if (patch <= 3)
234 hit.erry = 10.0;//1.0;//pSP->fSigmaY2;
235 else
236 hit.erry = 10.0;//pSP->fSigmaY2;
237 hit.errz = 1.0;//pSP->fSigmaZ2;
238
239 hit.index = i;
240 hit.counter = 0;
241
242 vec_hits.push_back(hit);
243 }
244
245#ifdef DRAW
246
247 TMarker *mhit = new TMarker(0.0, 0.0, 6);
248 mhit->SetMarkerColor(kBlue);
249
250 for (Int_t i = patch_first_hit_ind; i <= patch_last_hit_ind; i++)
251 {
252 YX->cd();
253 mhit->DrawMarker(vec_hits[i].y, vec_hits[i].x);
254
255 ZX->cd();
256 mhit->DrawMarker(vec_hits[i].z, vec_hits[i].x);
257 }
258
259 delete mhit;
260
261 YX->cd(); YX->Update(); YX->Draw();
262 ZX->cd(); ZX->Update(); ZX->Draw();
263
264#endif //DRAW
265
266 return;
267}
268
269// ----------------------------------------------------------------------------------
270void AliHLTTPCCATracker::CAFindPatchTracks(Int_t patch)
271{
272
273#ifdef DRAW
274
275 TLine *line = new TLine();
276 line->SetLineColor(kBlack);
277
278#endif //DRAW
279
280 Double_t dist01, dist12, dist02, mindist, chi2, minchi2;
281 CAHit *ph0, *ph1, *ph2, *phbest, *next;
282
283 Double_t x, y, z, ty, tz, cov_y, cov_ty, cov_yty, cov_z, cov_tz, cov_ztz, ye, ze, dx, dy, dz, sy2, sz2;
284
285 for (Int_t i1 = patch_first_hit_ind; i1 < patch_last_hit_ind; i1++) // <nhit-1
286 {
287 ph1 = &(vec_hits[i1]);
288 next = &(vec_hits[i1+1]);
289
290 mindist = 20.0;
291 minchi2 = 1.0;
292 phbest = NULL;
293
294 if (ph1->counter != 0){ // the previous hit exists
295 x = ph1->x;
296 y = ph1->y;
297 z = ph1->z;
298
299 dx = (ph1->x-ph0->x); // assume different x !!!
300 if (abs(dx) < 0.0001) dx = 0.0001;
301 ty = (ph1->y-ph0->y)/dx;
302 tz = (ph1->z-ph0->z)/dx;
303
304 cov_y = ph1->erry*ph1->erry;
305 cov_ty = (ph0->erry*ph0->erry + ph1->erry*ph1->erry)/(dx*dx);
306 cov_yty = (ph1->erry*ph1->erry)/dx;
307
308 cov_z = ph1->errz*ph1->errz;
309 cov_tz = (ph0->errz*ph0->errz + ph1->errz*ph1->errz)/(dx*dx);
310 cov_ztz = (ph1->errz*ph1->errz)/dx;
311
312 }
313
314 for (Int_t i2 = i1+1; i2 <= patch_last_hit_ind; i2++)
315 {
316 ph2 = &(vec_hits[i2]);
317 if (ph2->x < ph1->x) continue; // no backward direction!
318
319 dist12 = sqrt((ph1->x-ph2->x)*(ph1->x-ph2->x) + (ph1->y-ph2->y)*(ph1->y-ph2->y) + (ph1->z-ph2->z)*(ph1->z-ph2->z));
320 if (dist12 < mindist){
321 if (ph1->counter == 0){
322 mindist = dist12;
323 phbest = ph2;
324 }
325 else{
326 // calculate chi2 distance between the hit and the line
327 dx = (ph2->x-x); // assume different x !!!
328 if (abs(dx) < 0.0001) dx = 0.0001;
329
330 ye = y + ty*dx;
331 ze = z + tz*dx;
332
333 sy2 = cov_y + 2.0*cov_yty*dx + cov_ty*dx*dx;
334 sz2 = cov_z + 2.0*cov_ztz*dx + cov_tz*dx*dx;
335
336 dy = ph2->y - ye;
337 dz = ph2->z - ze;
338
339 chi2 = (dy*dy)/(sy2+ph2->erry*ph2->erry) + (dz*dz)/(sz2+ph2->errz*ph2->errz);
340
341 // check the straight line model
342 if (chi2 < minchi2){
343 mindist = dist12;
344 phbest = ph2;
345 }
346
347 }
348 }
349 }
350
351 if (phbest != NULL){// closest hit found
352
353 // exchange jbest with the next hit
354 Double_t x_tmp = next->x;
355 Double_t y_tmp = next->y;
356 Double_t z_tmp = next->z;
357
358 Double_t errx_tmp = next->errx;
359 Double_t erry_tmp = next->erry;
360 Double_t errz_tmp = next->errz;
361
362 Int_t index_tmp = next->index;
363 Int_t counter_tmp = next->counter;
364
365 next->x = phbest->x;
366 next->y = phbest->y;
367 next->z = phbest->z;
368
369 next->errx = phbest->errx;
370 next->erry = phbest->erry;
371 next->errz = phbest->errz;
372
373 next->index = phbest->index;
374 next->counter = phbest->counter;
375
376 phbest->x = x_tmp;
377 phbest->y = y_tmp;
378 phbest->z = z_tmp;
379
380 phbest->errx = errx_tmp;
381 phbest->erry = erry_tmp;
382 phbest->errz = errz_tmp;
383
384 phbest->index = index_tmp;
385 phbest->counter = counter_tmp;
386
387 if (ph1->counter == 0)
388 ph1->counter = 1;
389 next->counter = ph1->counter+1;
390
391 ph0 = ph1;
392 dist01 = mindist;
393
394#ifdef DRAW
395
396 YX->cd();
397 line->DrawLine(ph1->y, ph1->x, next->y, next->x);
398
399 ZX->cd();
400 line->DrawLine(ph1->z, ph1->x, next->z, next->x);
401
402#endif //DRAW
403
404 }
405
406 }
407
408#ifdef DRAW
409
410 delete line;
411
412 YX->cd(); YX->Update(); YX->Draw();
413 ZX->cd(); ZX->Update(); ZX->Draw();
414
415#endif //DRAW
416
417 CATrack tr;
418 tr.nhits = 0;
419 tr.chi2 = 0.0;
420 tr.ndf = -4;
421 tr.vec_ihits.clear();
422
423 patch_first_track_ind = vec_patch_tracks.size();
424 patch_last_track_ind = patch_first_track_ind - 1;
425
426 CAHit *my_phf;
427
428 // collect tracks
429 for (Int_t ihr = patch_last_hit_ind; ihr >= patch_first_hit_ind; ihr--){
430
431 CAHit* phit = &(vec_hits[ihr]);
432
433 if ((phit->counter < 5)&&(tr.nhits == 0))
434 continue; // do not collect short tracks
435
436 if ((phit->counter != 0)&&(tr.nhits == 0))// store the first hit
437 my_phf = phit;
438
439 if (phit->counter != 0){//add hit to the track
440 tr.nhits++;
441 tr.vec_ihits.push_back(ihr);
442
443 if (tr.nhits == 2){ // calculate initial track parameters
444 CAHit* phit0 = &(vec_hits[tr.vec_ihits[0]]);
445 CAHit* phit1 = &(vec_hits[tr.vec_ihits[1]]);
446
447 tr.x = phit1->x;
448 tr.y = phit1->y;
449 tr.z = phit1->z;
450
451 Double_t dx = (phit1->x-phit0->x); // assume different x !!!
452 if (abs(dx) < 0.0001) dx = 0.0001;
453 tr.ty = (phit1->y-phit0->y)/dx;
454 tr.tz = (phit1->z-phit0->z)/dx;
455
456 tr.cov_y = phit1->erry*phit1->erry;
457 tr.cov_ty = (phit0->erry*phit0->erry + phit1->erry*phit1->erry)/(dx*dx);
458 tr.cov_yty = (phit1->erry*phit1->erry)/dx;
459
460 tr.cov_z = phit1->errz*phit1->errz;
461 tr.cov_tz = (phit0->errz*phit0->errz + phit1->errz*phit1->errz)/(dx*dx);
462 tr.cov_ztz = (phit1->errz*phit1->errz)/dx;
463 }
464
465 if (tr.nhits > 2){
466 // propagate the track
467 Double_t dx = (phit->x-tr.x); // assume different x !!!
468 if (abs(dx) < 0.0001) dx = 0.0001;
469
470 Double_t ye = tr.y + tr.ty*dx;
471 Double_t ze = tr.z + tr.tz*dx;
472
473 Double_t cov_y = tr.cov_y + 2.0*tr.cov_yty*dx + tr.cov_ty*dx*dx;
474 Double_t cov_yty = tr.cov_yty + tr.cov_ty*dx;
475 Double_t cov_ty = tr.cov_ty;
476
477 Double_t cov_z = tr.cov_z + 2.0*tr.cov_ztz*dx + tr.cov_tz*dx*dx;
478 Double_t cov_ztz = tr.cov_ztz + tr.cov_tz*dx;
479 Double_t cov_tz = tr.cov_tz;
480
481
482 Double_t dy = phit->y - ye;
483 Double_t dz = phit->z - ze;
484
485 // add the measurement
486
487 Double_t w = 1.0/(phit->erry*phit->erry + cov_y);
488
489 tr.y = ye + cov_y*dy*w;
490 tr.ty = tr.ty + cov_yty*dy*w;
491
492 tr.cov_y = cov_y - cov_y*cov_y*w;
493 tr.cov_yty = cov_yty - cov_y*cov_yty*w;
494 tr.cov_ty = cov_ty -cov_yty*cov_yty*w;
495
496 tr.chi2 += dy*dy*w;
497 tr.ndf += 1;
498
499 w = 1.0/(phit->errz*phit->errz + cov_z);
500
501 tr.z = ze + cov_z*dz*w;
502 tr.tz = tr.tz + cov_ztz*dz*w;
503
504 tr.cov_z = cov_z - cov_z*cov_z*w;
505 tr.cov_ztz = cov_ztz - cov_z*cov_ztz*w;
506 tr.cov_tz = cov_tz -cov_ztz*cov_ztz*w;
507
508 tr.x += dx;
509
510 tr.chi2 += dz*dz*w;
511 tr.ndf += 1;
512 }
513 }
514 if (phit->counter == 1){// store the track
515
516 tr.good = 1;
517 tr.used = 0;
518 tr.next = -1;
519 tr.patch = patch;
520
521 Double_t trdist = sqrt((my_phf->x-phit->x)*(my_phf->x-phit->x)+
522 (my_phf->y-phit->y)*(my_phf->y-phit->y)+
523 (my_phf->z-phit->z)*(my_phf->z-phit->z));
524
525 if ((trdist > 1.0)&&(tr.chi2/tr.ndf < 10.0)){
526 vec_patch_tracks.push_back(tr);
527 patch_last_track_ind++;
528 }
529
530 // reset the track
531 tr.nhits = 0;
532 tr.chi2 = 0.0;
533 tr.ndf = -4;
534 tr.vec_ihits.clear();
535 }
536 }
537
538 // sort tracks according (currently) to number of hits
539 //sort(vec_patch_tracks.begin(), vec_patch_tracks.end(), compareNhits); //don't need to sort here
540
541
542#ifdef DRAW
543
544 TLine *trline = new TLine();
545 trline->SetLineColor(kGreen);
546
547 for (Int_t itr = patch_first_track_ind; itr <= patch_last_track_ind; itr++){
548
549 CATrack *ptr = &(vec_patch_tracks[itr]);
550
551 CAHit* phitf = &(vec_hits[ptr->vec_ihits.back()]);
552 CAHit* phitl = &(vec_hits[ptr->vec_ihits.front()]);
553
554 YX->cd();
555 Double_t yf = ptr->y + ptr->ty*(phitf->x-ptr->x);
556 Double_t yl = ptr->y + ptr->ty*(phitl->x-ptr->x);
557 trline->DrawLine(yf, phitf->x, yl, phitl->x);
558
559 ZX->cd();
560 Double_t zf = ptr->z + ptr->tz*(phitf->x-ptr->x);
561 Double_t zl = ptr->z + ptr->tz*(phitl->x-ptr->x);
562 trline->DrawLine(zf, phitf->x, zl, phitl->x);
563
564 }
565
566 delete trline;
567
568 YX->cd(); YX->Update(); YX->Draw();
569 ZX->cd(); ZX->Update(); ZX->Draw();
570
571#endif //DRAW
572
573 return;
574}
575
576// ----------------------------------------------------------------------------------
577void AliHLTTPCCATracker::CAFindSliceTracks()
578{
579
580 Int_t ntracks = vec_patch_tracks.size();
581
582 // merge patch tracks into slice tracks
583 // start with the longest tracks in the first patches
584
585 // sort tracks according (currently) to patch and within patch - number of hits
586 sort(vec_patch_tracks.begin(), vec_patch_tracks.end(), compareCATracks);
587
588 for (Int_t itr1 = 0; itr1 < ntracks-1; itr1++){
589
590 CATrack *ptr1 = &(vec_patch_tracks[itr1]);
591
592 Double_t maxdisty = 10.0;
593 Double_t maxdistz = 5.0;
594
595 Double_t mindist = 10.0;
596
597 Int_t next = -1;
598#ifdef XXX
599 for (Int_t itr2 = itr1+1; itr2 < ntracks; itr2++){
600
601 CATrack *ptr2 = &(vec_patch_tracks[itr2]);
602
603 if (ptr2->patch == ptr1->patch) continue; // the same patch - no prolongation
604 if (ptr2->patch > ptr1->patch+1) break; // sorted tracks - no prolongation skippeng over patches
605
606 if (ptr2->used == 1) continue; // already matched
607
608 CAHit* phitf2 = &(vec_hits[ptr2->vec_ihits.back()]);
609 CAHit* phitl2 = &(vec_hits[ptr2->vec_ihits.front()]);
610
611 // extrapolate tr1 to the both ends of tr2
612
613 Double_t dyf = ptr1->y + ptr1->ty*(phitf2->x-ptr1->x) - phitf2->y;
614 Double_t dyl = ptr1->y + ptr1->ty*(phitl2->x-ptr1->x) - phitl2->y;
615
616 Double_t dzf = ptr1->z + ptr1->tz*(phitf2->x-ptr1->x) - phitf2->z;
617 Double_t dzl = ptr1->z + ptr1->tz*(phitl2->x-ptr1->x) - phitl2->z;
618
619 // roughly similar tracks ?
620 if ((abs(dyf) > maxdisty)||(abs(dyl) > maxdisty)||(abs(dzf) > maxdistz)||(abs(dzl) > maxdistz)) continue;
621
622 // roughly parallel tracks ?
623 Double_t dist = sqrt((dyf-dyl)*(dyf-dyl)+(dzf-dzl)*(dzf-dzl));
624 if (dist > mindist) continue;
625
626 mindist = dist;
627 next = itr2;
628
629 }
630#endif
631
632 // found track prolongations?
633 if (next != -1){
634 CATrack *ptr2 = &(vec_patch_tracks[next]);
635
636 ptr1->next = next;
637 ptr2->used = 1;
638
639
640 }
641
642 }
643
644 CATrack tr;
645 tr.nhits = 0;
646 tr.chi2 = 0.0;
647 tr.ndf = -4;
648 tr.vec_ihits.clear();
649
650
651 //collect tracks
652 for (Int_t itr = 0; itr < ntracks; itr++){
653
654 CATrack *ptr = &(vec_patch_tracks[itr]);
655 if (ptr->used) continue; // start with a track not used in prolongations
656
657 Int_t first = 1;
658 do{
659 if (first == 1){
660 first = 0;
661 }
662 else{
663 ptr = &(vec_patch_tracks[ptr->next]);
664 }
665
666 for (Int_t ih = 0; ih < ptr->nhits; ih++){
667 tr.vec_ihits.push_back(ptr->vec_ihits[ih]);
668 tr.nhits++;
669 }
670 }while(ptr->next != -1);
671
672 //sort hits according to increasing x
673 sort(tr.vec_ihits.begin(), tr.vec_ihits.end(), compareCAHitsX);
674
675 vec_slice_tracks.push_back(tr);
676
677 tr.nhits = 0;
678 tr.chi2 = 0.0;
679 tr.ndf = -4;
680 tr.vec_ihits.clear();
681
682 }
683
684#ifdef DRAW
685
686 TLine *trline = new TLine();
687 trline->SetLineColor(kRed);
688
689#endif //DRAW
690
691#if WRITETRACKS
692 UInt_t size = 0;
693 UInt_t nTracks = 0;
694#endif
695
696
697 for (vector<CATrack>::iterator trIt = vec_slice_tracks.begin(); trIt != vec_slice_tracks.end(); trIt++){
698
699 CAHit* phit0 = &(vec_hits[trIt->vec_ihits[trIt->nhits-1]]);
700 CAHit* phit1 = &(vec_hits[trIt->vec_ihits[trIt->nhits-2]]);
701
702 trIt->x = phit1->x;
703 trIt->y = phit1->y;
704 trIt->z = phit1->z;
705
706 Double_t dx = (phit1->x-phit0->x); // assume different x !!!
707 if (abs(dx) < 0.0001) dx = 0.0001;
708 trIt->ty = (phit1->y-phit0->y)/dx;
709 trIt->tz = (phit1->z-phit0->z)/dx;
710
711 trIt->cov_y = phit1->erry*phit1->erry;
712 trIt->cov_ty = (phit0->erry*phit0->erry + phit1->erry*phit1->erry)/(dx*dx);
713 trIt->cov_yty = (phit1->erry*phit1->erry)/dx;
714
715 trIt->cov_z = phit1->errz*phit1->errz;
716 trIt->cov_tz = (phit0->errz*phit0->errz + phit1->errz*phit1->errz)/(dx*dx);
717 trIt->cov_ztz = (phit1->errz*phit1->errz)/dx;
718
719 for (Int_t ih = trIt->nhits-3; ih >= 0; ih--){
720
721 CAHit* phit = &(vec_hits[trIt->vec_ihits[ih]]);
722
723 // propagate the track
724 Double_t dx = (phit->x-trIt->x); // assume different x !!!
725 if (abs(dx) < 0.0001) dx = 0.0001;
726
727 Double_t ye = trIt->y + trIt->ty*dx;
728 Double_t ze = trIt->z + trIt->tz*dx;
729
730 Double_t cov_y = trIt->cov_y + 2.0*trIt->cov_yty*dx + trIt->cov_ty*dx*dx;
731 Double_t cov_yty = trIt->cov_yty + trIt->cov_ty*dx;
732 Double_t cov_ty = trIt->cov_ty;
733
734 Double_t cov_z = trIt->cov_z + 2.0*trIt->cov_ztz*dx + trIt->cov_tz*dx*dx;
735 Double_t cov_ztz = trIt->cov_ztz + trIt->cov_tz*dx;
736 Double_t cov_tz = trIt->cov_tz;
737
738 Double_t dy = phit->y - ye;
739 Double_t dz = phit->z - ze;
740
741 // add the measurement
742
743 Double_t w = 1.0/(phit->erry*phit->erry + cov_y);
744
745 trIt->y = ye + cov_y*dy*w;
746 trIt->ty = trIt->ty + cov_yty*dy*w;
747
748 trIt->cov_y = cov_y - cov_y*cov_y*w;
749 trIt->cov_yty = cov_yty - cov_y*cov_yty*w;
750 trIt->cov_ty = cov_ty -cov_yty*cov_yty*w;
751
752 trIt->chi2 += dy*dy*w;
753 trIt->ndf += 1;
754
755 w = 1.0/(phit->errz*phit->errz + cov_z);
756
757 trIt->z = ze + cov_z*dz*w;
758 trIt->tz = trIt->tz + cov_ztz*dz*w;
759
760 trIt->cov_z = cov_z - cov_z*cov_z*w;
761 trIt->cov_ztz = cov_ztz - cov_z*cov_ztz*w;
762 trIt->cov_tz = cov_tz -cov_ztz*cov_ztz*w;
763
764 trIt->chi2 += dz*dz*w;
765 trIt->ndf += 1;
766
767 trIt->x += dx;
768 }
769
770 trIt->good = 1;
771 trIt->used = 0;
772 trIt->next = -1;
773 trIt->patch = -1;
774
775 //if (trIt->chi2/trIt->ndf < 1.0)
776 //trIt->good = 0;
777
778
779 // JMT 2006/11/13 Write Tracks to container
780#if WRITETRACKS
781 CAHit* firstHit = &(vec_hits[trIt->vec_ihits[0]]);
782 CAHit* lastHit = &(vec_hits[trIt->vec_ihits[trIt->nhits-1]]);
783
784 Float_t xFirst = (Float_t) firstHit->x;
785 Float_t yFirst = (Float_t) trIt->y + trIt->ty*(xFirst-trIt->x);
786 float_t zFirst = (Float_t) trIt->z + trIt->tz*(xFirst-trIt->x);
787
788 Float_t xLast = (Float_t) lastHit->x;
789 Float_t yLast = (Float_t) trIt->y + trIt->ty*(xLast-trIt->x);
790 Float_t zLast = (Float_t) trIt->z + trIt->tz*(xLast-trIt->x);
791
792 fOutputPtr->fX = xFirst;
793 fOutputPtr->fY = yFirst;
794 fOutputPtr->fZ = zFirst;
795 fOutputPtr->fPt = -9876.0;
796 fOutputPtr->fPterr = 1;
797 fOutputPtr->fLastX = xLast;
798 fOutputPtr->fLastY = yLast;
799 fOutputPtr->fLastZ = zLast;
800 fOutputPtr->fPsi = atan(trIt->ty);
801 fOutputPtr->fTgl = trIt->tz;
802 fOutputPtr->fPsierr = 1;
803 fOutputPtr->fTglerr = 1;
804 fOutputPtr->fCharge = 1;
805 fOutputPtr->fNPoints = 0;
806
807 Byte_t *tmpP = (Byte_t *)fOutputPtr;
808
809 tmpP += sizeof(AliHLTTPCTrackSegmentData); //+fOutputPtr->fNPoints*sizeof(UInt_t);
810 size += sizeof(AliHLTTPCTrackSegmentData); //+fOutputPtr->fNPoints*sizeof(UInt_t);
811 fOutputPtr = (AliHLTTPCTrackSegmentData*)tmpP;
812 nTracks++;
813#endif
814
815#ifdef DRAW
816 CAHit* phitf = &(vec_hits[trIt->vec_ihits[0]]);
817 CAHit* phitl = &(vec_hits[trIt->vec_ihits[trIt->nhits-1]]);
818
819 Double_t xf = phitf->x;
820 Double_t yf = trIt->y + trIt->ty*(xf-trIt->x);
821 Double_t zf = trIt->z + trIt->tz*(xf-trIt->x);
822
823 Double_t xl = phitl->x;
824 Double_t yl = trIt->y + trIt->ty*(xl-trIt->x);
825 Double_t zl = trIt->z + trIt->tz*(xl-trIt->x);
826
827 YX->cd();
828 trline->DrawLine(yf, xf, yl, xl);
829
830 ZX->cd();
831 trline->DrawLine(zf, xf, zl, xl);
832
833#endif //DRAW
834
835 }
836
837#if WRITETRACKS
838 fOutputNTracks = nTracks;
839 fOutputSize = size;
840 cout << "NTRACKS=" << nTracks << endl;
841 cout << "SIZEoF=" << sizeof(AliHLTTPCTrackSegmentData) << endl;
842 cout << "SIZE=" << fOutputSize << endl;
843#endif
844
845#ifdef DRAW
846
847 delete trline;
848
849 YX->cd(); YX->Update(); YX->Draw();
850 ZX->cd(); ZX->Update(); ZX->Draw();
851
852#endif //DRAW
853
854 return;
855}