]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TRD/AliTRDtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Add(const ) function for Theo (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
... / ...
CommitLineData
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/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// The standard TRD tracker //
21// Based on Kalman filltering approach //
22// //
23// Authors: //
24// M. Ivanov (Marian.Ivanov@cern.ch) //
25// //
26///////////////////////////////////////////////////////////////////////////////
27
28#include <TBranch.h>
29#include <TFile.h>
30#include <TGraph.h>
31#include <TLinearFitter.h>
32#include <TObjArray.h>
33#include <TROOT.h>
34#include <TTree.h>
35#include <TTreeStream.h>
36
37#include "AliESDEvent.h"
38#include "AliESDtrack.h"
39//#include "AliAlignObj.h"
40#include "AliGeomManager.h"
41#include "AliRieman.h"
42#include "AliTrackPointArray.h"
43
44#include "AliTRDgeometry.h"
45#include "AliTRDpadPlane.h"
46#include "AliTRDgeometry.h"
47#include "AliTRDcluster.h"
48#include "AliTRDtrack.h"
49#include "AliTRDseed.h"
50#include "AliTRDcalibDB.h"
51#include "AliTRDCommonParam.h"
52#include "AliTRDtracker.h"
53#include "AliTRDReconstructor.h"
54#include "AliTRDrecoParam.h"
55#include "AliTRDCalibraFillHisto.h"
56
57ClassImp(AliTRDtracker)
58
59const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; //
60const Float_t AliTRDtracker::fgkLabelFraction = 0.8; //
61const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; //
62const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Maximum local sine of the azimuthal angle
63const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation
64
65//_____________________________________________________________________________
66AliTRDtracker::AliTRDtracker(AliTRDReconstructor *rec)
67 :AliTracker()
68 ,fReconstructor(rec)
69 ,fGeom(0)
70 ,fNclusters(0)
71 ,fClusters(0)
72 ,fNseeds(0)
73 ,fSeeds(0)
74 ,fNtracks(0)
75 ,fTracks(0)
76 ,fTimeBinsPerPlane(0)
77 ,fAddTRDseeds(kFALSE)
78 ,fNoTilt(kFALSE)
79 ,fHBackfit(NULL)
80 ,fHClSearch(NULL)
81 ,fHRefit(NULL)
82 ,fHX(NULL)
83 ,fHNCl(NULL)
84 ,fHNClTrack(NULL)
85 ,fHMinYPos(NULL)
86 ,fHMinYNeg(NULL)
87 ,fHMinZ(NULL)
88 ,fHMinD(NULL)
89 ,fHDeltaX(NULL)
90 ,fHXCl(NULL)
91 ,fDebugStreamer(0)
92{
93 //
94 // Default constructor
95 //
96
97 for (Int_t i = 0; i < kTrackingSectors; i++) {
98 fTrSec[i] = 0;
99 }
100
101 InitLogHists();
102
103}
104
105//_____________________________________________________________________________
106AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
107 :AliTracker(t)
108 ,fReconstructor(t.fReconstructor)
109 ,fGeom(0)
110 ,fNclusters(0)
111 ,fClusters(0)
112 ,fNseeds(0)
113 ,fSeeds(0)
114 ,fNtracks(0)
115 ,fTracks(0)
116 ,fTimeBinsPerPlane(0)
117 ,fAddTRDseeds(kFALSE)
118 ,fNoTilt(kFALSE)
119 ,fHBackfit(NULL)
120 ,fHClSearch(NULL)
121 ,fHRefit(NULL)
122 ,fHX(NULL)
123 ,fHNCl(NULL)
124 ,fHNClTrack(NULL)
125 ,fHMinYPos(NULL)
126 ,fHMinYNeg(NULL)
127 ,fHMinZ(NULL)
128 ,fHMinD(NULL)
129 ,fHDeltaX(NULL)
130 ,fHXCl(NULL)
131 ,fDebugStreamer(0)
132{
133 //
134 // Copy constructor
135 //
136
137}
138
139//_____________________________________________________________________________
140AliTRDtracker::AliTRDtracker(const TFile */*geomfile*/, AliTRDReconstructor *rec)
141 :AliTracker()
142 ,fReconstructor(rec)
143 ,fGeom(0)
144 ,fNclusters(0)
145 ,fClusters(new TObjArray(2000))
146 ,fNseeds(0)
147 ,fSeeds(new TObjArray(2000))
148 ,fNtracks(0)
149 ,fTracks(new TObjArray(1000))
150 ,fTimeBinsPerPlane(0)
151 ,fAddTRDseeds(kFALSE)
152 ,fNoTilt(kFALSE)
153 ,fHBackfit(NULL)
154 ,fHClSearch(NULL)
155 ,fHRefit(NULL)
156 ,fHX(NULL)
157 ,fHNCl(NULL)
158 ,fHNClTrack(NULL)
159 ,fHMinYPos(NULL)
160 ,fHMinYNeg(NULL)
161 ,fHMinZ(NULL)
162 ,fHMinD(NULL)
163 ,fHDeltaX(NULL)
164 ,fHXCl(NULL)
165 ,fDebugStreamer(0)
166{
167 //
168 // Main constructor
169 //
170
171 TDirectory *savedir = gDirectory;
172
173 fGeom = new AliTRDgeometry();
174
175 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
176 Int_t trS = geomS;
177 fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
178 }
179
180 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
181 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
182 if (tiltAngle < 0.1) {
183 fNoTilt = kTRUE;
184 }
185
186 if (!AliTRDcalibDB::Instance()) {
187 AliFatal("Could not get calibration object");
188 }
189 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
190
191 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
192
193 savedir->cd();
194
195 InitLogHists();
196
197}
198
199//_____________________________________________________________________________
200AliTRDtracker::~AliTRDtracker()
201{
202 //
203 // Destructor of AliTRDtracker
204 //
205
206 if (fClusters) {
207 fClusters->Delete();
208 delete fClusters;
209 }
210
211 if (fTracks) {
212 fTracks->Delete();
213 delete fTracks;
214 }
215
216 if (fSeeds) {
217 fSeeds->Delete();
218 delete fSeeds;
219 }
220
221 if (fGeom) {
222 delete fGeom;
223 }
224
225 for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
226 delete fTrSec[geomS];
227 }
228
229 if (fDebugStreamer) {
230 delete fDebugStreamer;
231 }
232
233}
234
235//_____________________________________________________________________________
236Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
237{
238 //
239 // Transform internal TRD ID to global detector ID
240 //
241
242 Int_t isector = fGeom->GetSector(lid);
243 Int_t istack = fGeom->GetStack(lid);
244 Int_t ilayer = fGeom->GetLayer(lid);
245
246 AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
247 switch (ilayer) {
248 case 0:
249 iLayer = AliGeomManager::kTRD1;
250 break;
251 case 1:
252 iLayer = AliGeomManager::kTRD2;
253 break;
254 case 2:
255 iLayer = AliGeomManager::kTRD3;
256 break;
257 case 3:
258 iLayer = AliGeomManager::kTRD4;
259 break;
260 case 4:
261 iLayer = AliGeomManager::kTRD5;
262 break;
263 case 5:
264 iLayer = AliGeomManager::kTRD6;
265 break;
266 };
267
268 Int_t modId = isector * fGeom->Nstack() + istack;
269 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
270
271 return volid;
272
273}
274
275//_____________________________________________________________________________
276Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
277{
278 //
279 // Transform global detector ID to local detector ID
280 //
281
282 Int_t modId = 0;
283 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
284
285 Int_t isector = modId / fGeom->Nstack();
286 Int_t istack = modId % fGeom->Nstack();
287 Int_t iLayer = -1;
288
289 switch (layerId) {
290 case AliGeomManager::kTRD1:
291 iLayer = 0;
292 break;
293 case AliGeomManager::kTRD2:
294 iLayer = 1;
295 break;
296 case AliGeomManager::kTRD3:
297 iLayer = 2;
298 break;
299 case AliGeomManager::kTRD4:
300 iLayer = 3;
301 break;
302 case AliGeomManager::kTRD5:
303 iLayer = 4;
304 break;
305 case AliGeomManager::kTRD6:
306 iLayer = 5;
307 break;
308 default:
309 iLayer =-1;
310 }
311
312 if (iLayer < 0) {
313 return -1;
314 }
315
316 Int_t lid = fGeom->GetDetector(iLayer,istack,isector);
317
318 return lid;
319
320}
321
322//_____________________________________________________________________________
323Bool_t AliTRDtracker::AdjustSector(AliTRDtrack* const track) const
324{
325 //
326 // Rotates the track when necessary
327 //
328
329 Double_t alpha = AliTRDgeometry::GetAlpha();
330 Double_t y = track->GetY();
331 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
332
333 if (y > ymax) {
334 if (!track->Rotate( alpha)) {
335 return kFALSE;
336 }
337 }
338 else if (y < -ymax) {
339 if (!track->Rotate(-alpha)) {
340 return kFALSE;
341 }
342 }
343
344 return kTRUE;
345
346}
347
348//_____________________________________________________________________________
349AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack * const track, Int_t plane
350 , Int_t timebin, UInt_t &index)
351{
352 //
353 // Try to find cluster in the backup list
354 //
355
356 AliTRDcluster *cl =0;
357 const Int_t *indexes = track->GetBackupIndexes();
358
359 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
360 if (indexes[i] == 0) {
361 break;
362 }
363 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
364 if (!cli) {
365 break;
366 }
367 if (cli->GetLocalTimeBin() != timebin) {
368 continue;
369 }
370 Int_t ilayer = fGeom->GetLayer(cli->GetDetector());
371 if (ilayer == plane) {
372 cl = cli;
373 index = indexes[i];
374 break;
375 }
376 }
377
378 return cl;
379
380}
381
382//_____________________________________________________________________________
383Int_t AliTRDtracker::GetLastPlane(AliTRDtrack * const track)
384{
385 //
386 // Return last updated plane
387 //
388
389 Int_t lastplane = 0;
390 const Int_t *indexes = track->GetBackupIndexes();
391
392 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
393 AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
394 if (!cli) {
395 break;
396 }
397 Int_t ilayer = fGeom->GetLayer(cli->GetDetector());
398 if (ilayer > lastplane) {
399 lastplane = ilayer;
400 }
401 }
402
403 return lastplane;
404
405}
406
407//_____________________________________________________________________________
408Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
409{
410 //
411 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
412 // backpropagated by the TPC tracker. Each seed is first propagated
413 // to the TRD, and then its prolongation is searched in the TRD.
414 // If sufficiently long continuation of the track is found in the TRD
415 // the track is updated, otherwise it's stored as originaly defined
416 // by the TPC tracker.
417 //
418
419 Int_t found = 0; // number of tracks found
420 Float_t foundMin = 20.0;
421
422 Int_t nSeed = event->GetNumberOfTracks();
423 if(!nSeed){
424 // run stand alone tracking
425 if (fReconstructor->IsSeeding()) Clusters2Tracks(event);
426 return 0;
427 }
428
429 Float_t *quality = new Float_t[nSeed];
430 Int_t *index = new Int_t[nSeed];
431 for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
432 AliESDtrack *seed = event->GetTrack(iSeed);
433 Double_t covariance[15];
434 seed->GetExternalCovariance(covariance);
435 quality[iSeed] = covariance[0] + covariance[2];
436 }
437 // Sort tracks according to covariance of local Y and Z
438 TMath::Sort(nSeed,quality,index,kFALSE);
439
440 // Backpropagate all seeds
441 for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
442
443 // Get the seeds in sorted sequence
444 AliESDtrack *seed = event->GetTrack(index[iSeed]);
445 fHBackfit->Fill(0); // All seeds
446
447 // Check the seed status
448 ULong_t status = seed->GetStatus();
449 if ((status & AliESDtrack::kTPCout) == 0) {
450 fHBackfit->Fill(1); // TPC outer edge reached
451 continue;
452 }
453 if ((status & AliESDtrack::kTRDout) != 0) {
454 fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
455 continue;
456 }
457
458 // Do the back prolongation
459 Int_t lbl = seed->GetLabel();
460 AliTRDtrack *track = new AliTRDtrack(*seed);
461 track->SetSeedLabel(lbl);
462 seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
463 fNseeds++;
464 Float_t p4 = track->GetC();
465 Int_t expectedClr = FollowBackProlongation(*track);
466 fHBackfit->Fill(3); // Back prolongation done
467 fHX->Fill(track->GetX());
468
469 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
470 (track->Pt() > 0.8)) {
471
472 fHBackfit->Fill(4);
473
474 //
475 // Make backup for back propagation
476 //
477
478 Int_t foundClr = track->GetNumberOfClusters();
479 if (foundClr >= foundMin) {
480 track->CookdEdx();
481 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
482 CookLabel(track,1 - fgkLabelFraction);
483 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
484
485
486 // Sign only gold tracks
487 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
488 if ((seed->GetKinkIndex(0) == 0) &&
489 (track->Pt() < 1.5)) {
490 UseClusters(track);
491 }
492 }
493 Bool_t isGold = kFALSE;
494
495 // Full gold track
496 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
497 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
498 if (track->GetBackupTrack()) {
499 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
500 }
501 isGold = kTRUE;
502 //fHBackfit->Fill()
503 }
504
505 // Almost gold track
506 if ((!isGold) && (track->GetNCross() == 0) &&
507 (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
508 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
509 if (track->GetBackupTrack()) {
510 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
511 }
512 isGold = kTRUE;
513 }
514
515 if ((!isGold) && (track->GetBackupTrack())) {
516 if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
517 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
518 isGold = kTRUE;
519 }
520 }
521
522 if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
523 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
524 }
525 }
526 }
527 /**/
528
529 /**/
530 // Debug part of tracking
531 TTreeSRedirector &cstream = *fDebugStreamer;
532 Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
533 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
534 if (track->GetBackupTrack()) {
535 cstream << "Tracks"
536 << "EventNrInFile=" << eventNrInFile
537 << "ESD.=" << seed
538 << "trd.=" << track
539 << "trdback.=" << track->GetBackupTrack()
540 << "\n";
541 }
542 else {
543 cstream << "Tracks"
544 << "EventNrInFile=" << eventNrInFile
545 << "ESD.=" << seed
546 << "trd.=" << track
547 << "trdback.=" << track
548 << "\n";
549 }
550 }
551 /**/
552
553 // Propagation to the TOF (I.Belikov)
554 if (track->GetStop() == kFALSE) {
555 fHBackfit->Fill(5);
556
557 Double_t xtof = 371.0;
558 Double_t xTOF0 = 370.0;
559
560 Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
561 if (TMath::Abs(c2) >= 0.99) {
562 fHBackfit->Fill(6);
563 delete track;
564 continue;
565 }
566
567 PropagateToX(*track,xTOF0,fgkMaxStep);
568
569 // Energy losses taken to the account - check one more time
570 c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
571 if (TMath::Abs(c2) >= 0.99) {
572 fHBackfit->Fill(7);
573 delete track;
574 continue;
575 }
576
577 //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
578 // fHBackfit->Fill(7);
579 //delete track;
580 // continue;
581 //}
582
583 Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
584 Double_t y;
585 track->GetYAt(xtof,GetBz(),y);
586 if (y > ymax) {
587 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
588 fHBackfit->Fill(8);
589 delete track;
590 continue;
591 }
592 }
593 else if (y < -ymax) {
594 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
595 fHBackfit->Fill(9);
596 delete track;
597 continue;
598 }
599 }
600
601 if (track->PropagateTo(xtof)) {
602 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
603 fHBackfit->Fill(10);
604
605 seed->SetNumberOfTRDslices(AliTRDCalPID::kNSlicesLQ);
606 for (Int_t i = 0; i < AliTRDtrack::kNplane; i++) {
607 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
608 seed->SetTRDslice(track->GetPIDsignals(i,j),i,j);
609 }
610 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
611 }
612 //seed->SetTRDtrack(new AliTRDtrack(*track));
613 if (track->GetNumberOfClusters() > foundMin) {
614 fHBackfit->Fill(11);
615 found++;
616 }
617 }
618 }
619 else {
620 fHBackfit->Fill(12);
621
622 if ((track->GetNumberOfClusters() > 15) &&
623 (track->GetNumberOfClusters() > 0.5*expectedClr)) {
624 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
625 fHBackfit->Fill(13);
626
627 //seed->SetStatus(AliESDtrack::kTRDStop);
628
629 seed->SetNumberOfTRDslices(AliTRDCalPID::kNSlicesLQ);
630 for (Int_t i = 0; i < AliTRDtrack::kNplane; i++) {
631 for (Int_t j = 0; j <AliTRDCalPID::kNSlicesLQ; j++) {
632 seed->SetTRDslice(track->GetPIDsignals(i,j),i,j);
633 }
634 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
635 }
636 //seed->SetTRDtrack(new AliTRDtrack(*track));
637 found++;
638 }
639 }
640
641 seed->SetTRDQuality(track->StatusForTOF());
642 seed->SetTRDBudget(track->GetBudget(0));
643
644 fHBackfit->Fill(14);
645 delete track;
646 }
647
648 AliInfo(Form("Number of seeds: %d",fNseeds));
649 AliInfo(Form("Number of back propagated TRD tracks: %d",found));
650
651 fSeeds->Clear();
652 fNseeds = 0;
653
654 delete [] index;
655 delete [] quality;
656
657 SaveLogHists();
658
659 return 0;
660}
661
662//_____________________________________________________________________________
663Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
664{
665 //
666 // Refits tracks within the TRD. The ESD event is expected to contain seeds
667 // at the outer part of the TRD.
668 // The tracks are propagated to the innermost time bin
669 // of the TRD and the ESD event is updated
670 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
671 //
672
673 //Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
674 //Float_t foundMin = fgkMinClustersInTrack * timeBins;
675 Int_t nseed = 0;
676 Int_t found = 0;
677 Int_t pidQ = 0;
678 //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
679 AliTRDtrack seed2;
680
681 // Calibration fill 2D
682 AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
683 if (!calibra) {
684 AliInfo("Could not get Calibra instance\n");
685 }
686
687 Int_t n = event->GetNumberOfTracks();
688 for (Int_t i = 0; i < n; i++) {
689
690 AliESDtrack *seed = event->GetTrack(i);
691 new (&seed2) AliTRDtrack(*seed);
692 fHRefit->Fill(0);
693
694 if (seed2.GetX() < 270.0) {
695 seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
696 fHRefit->Fill(1);
697 continue;
698 }
699
700 ULong_t status = seed->GetStatus();
701 if ((status & AliESDtrack::kTRDout) == 0) {
702 fHRefit->Fill(2);
703 continue;
704 }
705 if ((status & AliESDtrack::kTRDin) != 0) {
706 fHRefit->Fill(3);
707 continue;
708 }
709
710 nseed++;
711 fHRefit->Fill(4);
712
713 seed2.ResetCovariance(50.0);
714
715 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
716 const Int_t *indexes2 = seed2.GetIndexes();
717 for (Int_t l = 0; l < AliTRDtrack::kNplane;++l) {
718 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ;j++) {
719 pt->SetPIDsignals(seed2.GetPIDsignals(l,j),l,j);
720 }
721 pt->SetPIDTimBin(seed2.GetPIDTimBin(l),l);
722 }
723
724 Int_t *indexes3 = const_cast<Int_t *>(pt->GetBackupIndexes());
725 for (Int_t l = 0; l < 200;++l) {
726 if (indexes2[l] == 0) {
727 break;
728 }
729 indexes3[l] = indexes2[l];
730 }
731
732 FollowProlongation(*pt);
733 pt->CookdEdx();
734 pt->CookdEdxTimBin(seed->GetID());
735
736 //calculate PID methods
737 pt->SetPIDMethod(AliTRDtrack::kLQ);
738 pt->CookPID(pidQ);
739 seed->SetTRDpid(pt->GetPID());
740 seed->SetTRDntracklets(pidQ<<3);
741
742 // update calibration
743 if(calibra->GetHisto2d()) calibra->UpdateHistograms(pt);
744 found++;
745
746 Double_t xTPC = 250.0;
747 if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
748 seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
749 fHRefit->Fill(5);
750
751 for (Int_t l = 0; l < AliTRDtrack::kNplane; ++l) {
752 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
753 seed->SetTRDslice(pt->GetPIDsignals(l,j),l,j);
754 }
755 seed->SetTRDTimBin(pt->GetPIDTimBin(l),l);
756 }
757 } else {
758 // If not prolongation to TPC - propagate without update
759 fHRefit->Fill(5);
760 AliTRDtrack *seed2t = new AliTRDtrack(*seed);
761 seed2t->ResetCovariance(5.0);
762 AliTRDtrack *pt2 = new AliTRDtrack(*seed2t,seed2t->GetAlpha());
763 delete seed2t;
764
765 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
766 pt2->CookdEdx();
767 pt2->CookdEdxTimBin(seed->GetID());
768 seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
769 fHRefit->Fill(6);
770
771 for (Int_t l = 0; l < AliTRDtrack::kNplane; ++l) {
772 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
773 seed->SetTRDslice(pt2->GetPIDsignals(l,j),l,j);
774 }
775 seed->SetTRDTimBin(pt2->GetPIDTimBin(l),l);
776 }
777 }
778
779 // Add TRD track to ESDfriendTrack - maybe this tracks are
780 // not useful for post-processing - TODO make decision
781 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
782 seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
783 }
784 delete pt2;
785
786 }
787
788 // Add TRD track to ESDfriendTrack
789 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
790 seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
791 }
792 delete pt;
793
794 }
795
796 AliInfo(Form("Number of loaded seeds: %d",nseed));
797 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
798
799 SaveLogHists();
800
801 return 0;
802
803}
804
805//_____________________________________________________________________________
806Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
807{
808 //
809 // Starting from current position on track=t this function tries
810 // to extrapolate the track up to timeBin=0 and to confirm prolongation
811 // if a close cluster is found. Returns the number of clusters
812 // expected to be found in sensitive layers
813 // GeoManager used to estimate mean density
814 //
815
816 Int_t sector;
817 Int_t lastplane = GetLastPlane(&t);
818 Double_t xx0 = 0.0;
819 Double_t xrho= 0.0;
820 Int_t expectedNumberOfClusters = 0;
821
822 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
823
824 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
825 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
826
827 // Propagate track close to the plane if neccessary
828 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
829 if (currentx < (-fgkMaxStep + t.GetX())) {
830 // Propagate closer to chamber - safety space fgkMaxStep
831 if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
832 break;
833 }
834 }
835
836 if (!AdjustSector(&t)) {
837 break;
838 }
839
840 // Get material budget
841 Double_t xyz0[3];
842 Double_t xyz1[3];
843 Double_t param[7];
844 Double_t x;
845 Double_t y;
846 Double_t z;
847
848 // Starting global position
849 t.GetXYZ(xyz0);
850 // End global position
851 x = fTrSec[0]->GetLayer(row0)->GetX();
852 if (!t.GetProlongation(x,y,z)) {
853 break;
854 }
855 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
856 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
857 xyz1[2] = z;
858 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
859 xrho= param[0]*param[4];
860 xx0 = param[1]; // Get mean propagation parameters
861
862 // Flags for marking the track momentum and direction. The track is
863 // marked only if it has at least 1 cluster picked up in the current
864 // chamber.
865 Bool_t kUPDATED = kFALSE;
866 Bool_t kMARKED = kFALSE;
867
868 // Propagate and update
869 sector = t.GetSector();
870 //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
871 for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
872
873 // Mark track kinematics
874 if (itime > 10 && kUPDATED && !kMARKED) {
875 t.SetTrackSegmentDirMom(iplane);
876 kMARKED = kTRUE;
877 }
878
879 Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
880 expectedNumberOfClusters++;
881 t.SetNExpected(t.GetNExpected() + 1);
882 if (t.GetX() > 345.0) {
883 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
884 }
885 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
886 AliTRDcluster *cl = 0;
887 UInt_t index = 0;
888 Double_t maxChi2 = fgkMaxChi2;
889 x = timeBin.GetX();
890
891 if (timeBin) {
892
893 AliTRDcluster *cl0 = timeBin[0];
894 if (!cl0) {
895 // No clusters in given time bin
896 continue;
897 }
898
899 Int_t layer = fGeom->GetLayer(cl0->GetDetector());
900 if (layer > lastplane) {
901 continue;
902 }
903
904 Int_t timebin = cl0->GetLocalTimeBin();
905 AliTRDcluster *cl2 = GetCluster(&t,layer,timebin,index);
906
907 if (cl2) {
908 cl = cl2;
909 //Double_t h01 = GetTiltFactor(cl); //I.B's fix
910 //maxChi2=t.GetPredictedChi2(cl,h01);
911 }
912
913 if (cl) {
914
915 //if (cl->GetNPads()<5)
916 Double_t dxsample = timeBin.GetdX();
917 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
918 Double_t h01 = GetTiltFactor(cl);
919 Int_t det = cl->GetDetector();
920 Int_t llayer = fGeom->GetLayer(det);
921
922 if (t.GetX() > 345.0) {
923 t.SetNLast(t.GetNLast() + 1);
924 t.SetChi2Last(t.GetChi2Last() + maxChi2);
925 }
926
927 Double_t xcluster = cl->GetX();
928 t.PropagateTo(xcluster,xx0,xrho);
929 if (!AdjustSector(&t)) {
930 break; //I.B's fix
931 }
932
933 maxChi2 = t.GetPredictedChi2(cl,h01);
934 if (maxChi2 < 1e+10) {
935 if (!t.UpdateMI(cl,maxChi2,index,h01,llayer)) {
936 // ????
937 }
938 else {
939 //SetCluster(cl, GetNumberOfClusters()-1);
940 kUPDATED = kTRUE;
941 }
942 }
943
944 }
945
946 }
947
948 }
949
950 }
951
952 return expectedNumberOfClusters;
953
954}
955
956//_____________________________________________________________________________
957Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
958{
959 //
960 // Starting from current radial position of track <t> this function
961 // extrapolates the track up to the outer timebin and in the sensitive
962 // layers confirms prolongation if a close cluster is found.
963 // Returns the number of clusters expected to be found in sensitive layers
964 // Uses the geomanager for material description
965 //
966 // return number of assigned clusters ?
967 //
968
969 Int_t sector;
970
971 Double_t xx0 = 0.0;
972 Double_t xrho = 0.0;
973
974 Float_t ratio0 = 0.0;
975
976 Int_t expectedNumberOfClusters = 0;
977
978 AliTRDtracklet tracklet;
979
980 const Int_t kNclusters = 1000;
981 Int_t clusters[kNclusters];
982 for (Int_t i = 0; i < kNclusters; i++) {
983 clusters[i] = -1;
984 }
985
986// // Calibration fill 2D
987// AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
988// if (!calibra) {
989// AliInfo("Could not get Calibra instance\n");
990// }
991// if (calibra->GetMITracking()) {
992// calibra->ResetTrack();
993// }
994
995 // Loop through the TRD planes
996 for (Int_t iplane = 0; iplane < AliTRDtrack::kNplane; iplane++) {
997
998 Int_t hb = iplane * 10;
999 fHClSearch->Fill(hb);
1000
1001 // Get the global time bin numbers for the first an last
1002 // local time bin of the given plane
1003 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1004 Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
1005
1006 // Get the X coordinates of the propagation layer for the first time bin
1007 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1008 if (currentx < t.GetX()) {
1009 fHClSearch->Fill(hb+1);
1010 continue;
1011 }
1012
1013 // Propagate closer to the current chamber if neccessary
1014 if (currentx > (fgkMaxStep + t.GetX())) {
1015 if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
1016 fHClSearch->Fill(hb+2);
1017 break;
1018 }
1019 }
1020
1021 // Rotate track to adjacent sector if neccessary
1022 if (!AdjustSector(&t)) {
1023 fHClSearch->Fill(hb+3);
1024 break;
1025 }
1026
1027 // Check whether azimuthal angle is getting too large
1028 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
1029 fHClSearch->Fill(hb+4);
1030 break;
1031 }
1032
1033 Double_t xyz0[3];
1034 Double_t xyz1[3];
1035 Double_t param[7];
1036 Double_t x;
1037 Double_t y;
1038 Double_t z;
1039 // Global start position (beginning of chamber)
1040 t.GetXYZ(xyz0);
1041 // X-position of the end of the chamber
1042 x = fTrSec[0]->GetLayer(rowlast)->GetX();
1043 // Get local Y and Z at the X-position of the end of the chamber
1044 if (!t.GetProlongation(x,y,z)) {
1045 fHClSearch->Fill(hb+5);
1046 break;
1047 }
1048 // Global end position (end of chamber)
1049 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1050 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1051 xyz1[2] = z;
1052
1053 // Calculate the mean material budget along the path inside the chamber
1054 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1055 // The mean propagation parameters (density*length and radiation length)
1056 xrho = param[0]*param[4];
1057 xx0 = param[1];
1058
1059 // Find the clusters and tracklet along the path inside the chamber
1060 sector = t.GetSector();
1061 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
1062 fHNCl->Fill(tracklet.GetN());
1063
1064 // Discard if in less than 1/3 of the available timebins
1065 // clusters are found
1066 if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
1067 fHClSearch->Fill(hb+6);
1068 continue;
1069 }
1070
1071 //
1072 // Propagate and update track
1073 //
1074 for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1075 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
1076 expectedNumberOfClusters++;
1077 t.SetNExpected(t.GetNExpected() + 1);
1078 if (t.GetX() > 345.0) {
1079 t.SetNExpectedLast(t.GetNExpectedLast() + 1);
1080 }
1081 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
1082 AliTRDcluster *cl = 0;
1083 UInt_t index = 0;
1084 Double_t maxChi2 = fgkMaxChi2;
1085 x = timeBin.GetX();
1086
1087 if (timeBin) {
1088 if (clusters[ilayer] > 0) {
1089 index = clusters[ilayer];
1090 cl = (AliTRDcluster *)GetCluster(index);
1091 //Double_t h01 = GetTiltFactor(cl); // I.B's fix
1092 //maxChi2=t.GetPredictedChi2(cl,h01); //
1093 }
1094
1095 if (cl) {
1096 //if (cl->GetNPads() < 5)
1097 Double_t dxsample = timeBin.GetdX();
1098 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1099 Double_t h01 = GetTiltFactor(cl);
1100 Int_t det = cl->GetDetector();
1101 Int_t layer = fGeom->GetLayer(det);
1102 if (t.GetX() > 345.0) {
1103 t.SetNLast(t.GetNLast() + 1);
1104 t.SetChi2Last(t.GetChi2Last() + maxChi2);
1105 }
1106 Double_t xcluster = cl->GetX();
1107 t.PropagateTo(xcluster,xx0,xrho);
1108 maxChi2 = t.GetPredictedChi2(cl,h01);
1109
1110 if (maxChi2<1e+10)
1111 if (!t.UpdateMI(cl,maxChi2,index,h01,layer)) {
1112 if (!t.Update(cl,maxChi2,index,h01)) {
1113 // ????
1114 }
1115 } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
1116
1117
1118// if (calibra->GetMITracking()) {
1119// calibra->UpdateHistograms(cl,&t);
1120// }
1121
1122 // Reset material budget if 2 consecutive gold
1123 if (layer > 0) {
1124 if ((t.GetTracklets(layer).GetN() + t.GetTracklets(layer-1).GetN()) > 20) {
1125 t.SetBudget(2,0.0);
1126 }
1127 }
1128
1129 }
1130
1131 }
1132
1133 }
1134
1135 ratio0 = ncl / Float_t(fTimeBinsPerPlane);
1136 Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
1137 if ((tracklet.GetChi2() < 18.0) &&
1138 (ratio0 > 0.8) &&
1139 (ratio1 > 0.6) &&
1140 (ratio0+ratio1 > 1.5) &&
1141 (t.GetNCross() == 0) &&
1142 (TMath::Abs(t.GetSnp()) < 0.85) &&
1143 (t.GetNumberOfClusters() > 20)){
1144 //if (ratio0 > 0.8) {
1145 t.MakeBackupTrack(); // Make backup of the track until is gold
1146 }
1147
1148 }
1149
1150 return expectedNumberOfClusters;
1151
1152}
1153
1154//_____________________________________________________________________________
1155Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
1156{
1157 //
1158 // Starting from current X-position of track <t> this function
1159 // extrapolates the track up to radial position <xToGo>.
1160 // Returns 1 if track reaches the plane, and 0 otherwise
1161 //
1162
1163 const Double_t kEpsilon = 0.00001;
1164
1165 // Current track X-position
1166 Double_t xpos = t.GetX();
1167
1168 // Direction: inward or outward
1169 Double_t dir = (xpos < xToGo) ? 1.0 : -1.0;
1170
1171 while (((xToGo - xpos) * dir) > kEpsilon) {
1172
1173 Double_t xyz0[3];
1174 Double_t xyz1[3];
1175 Double_t param[7];
1176 Double_t x;
1177 Double_t y;
1178 Double_t z;
1179
1180 // The next step size
1181 Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1182
1183 // Get the global position of the starting point
1184 t.GetXYZ(xyz0);
1185
1186 // X-position after next step
1187 x = xpos + step;
1188
1189 // Get local Y and Z at the X-position of the next step
1190 if (!t.GetProlongation(x,y,z)) {
1191 return 0; // No prolongation possible
1192 }
1193
1194 // The global position of the end point of this prolongation step
1195 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1196 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1197 xyz1[2] = z;
1198
1199 // Calculate the mean material budget between start and
1200 // end point of this prolongation step
1201 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1202
1203 // Propagate the track to the X-position after the next step
1204 if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
1205 return 0;
1206 }
1207
1208 // Rotate the track if necessary
1209 AdjustSector(&t);
1210
1211 // New track X-position
1212 xpos = t.GetX();
1213
1214 }
1215
1216 return 1;
1217
1218}
1219
1220//_____________________________________________________________________________
1221Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1222{
1223 //
1224 // Fills clusters into TRD tracking_sectors
1225 // Note that the numbering scheme for the TRD tracking_sectors
1226 // differs from that of TRD sectors
1227 //
1228
1229
1230 if (ReadClusters(fClusters, cTree)) {
1231 AliError("Problem with reading the clusters !");
1232 return 1;
1233 }
1234 Int_t ncl = fClusters->GetEntriesFast();
1235 fNclusters = ncl;
1236 AliInfo(Form("Sorting %d clusters",ncl));
1237
1238 UInt_t index;
1239 while (ncl--) {
1240
1241 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1242 Int_t detector = c->GetDetector();
1243 Int_t localTimeBin = c->GetLocalTimeBin();
1244 Int_t sector = fGeom->GetSector(detector);
1245 Int_t layer = fGeom->GetLayer(detector);
1246 Int_t trackingSector = sector;
1247
1248 //if (c->GetQ() > 10) {
1249 // Int_t stack = fGeom->GetStack(detector);
1250 //}
1251
1252 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(layer,localTimeBin);
1253 if (gtb < 0) {
1254 continue;
1255 }
1256 Int_t trLayer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1257
1258 index = ncl;
1259
1260 fHXCl->Fill(c->GetX());
1261
1262 fTrSec[trackingSector]->GetLayer(trLayer)->SetX(c->GetX());
1263 fTrSec[trackingSector]->GetLayer(trLayer)->InsertCluster(c,index);
1264
1265 }
1266
1267 return 0;
1268
1269}
1270
1271//_____________________________________________________________________________
1272void AliTRDtracker::UnloadClusters()
1273{
1274 //
1275 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1276 //
1277
1278 Int_t i;
1279 Int_t nentr;
1280
1281 nentr = fClusters->GetEntriesFast();
1282 for (i = 0; i < nentr; i++) {
1283 delete fClusters->RemoveAt(i);
1284 }
1285 fNclusters = 0;
1286
1287 nentr = fSeeds->GetEntriesFast();
1288 for (i = 0; i < nentr; i++) {
1289 delete fSeeds->RemoveAt(i);
1290 }
1291
1292 nentr = fTracks->GetEntriesFast();
1293 for (i = 0; i < nentr; i++) {
1294 delete fTracks->RemoveAt(i);
1295 }
1296
1297 Int_t nsec = AliTRDgeometry::kNsector;
1298 for (i = 0; i < nsec; i++) {
1299 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1300 fTrSec[i]->GetLayer(pl)->Clear();
1301 }
1302 }
1303
1304}
1305
1306//_____________________________________________________________________________
1307Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd)
1308{
1309 //
1310 // Creates seeds using clusters between position inner plane and outer plane
1311 //
1312
1313 const Double_t kMaxTheta = 1.0;
1314 const Double_t kMaxPhi = 2.0;
1315
1316 const Double_t kRoad0y = 6.0; // Road for middle cluster
1317 const Double_t kRoad0z = 8.5; // Road for middle cluster
1318
1319 const Double_t kRoad1y = 2.0; // Road in y for seeded cluster
1320 const Double_t kRoad1z = 20.0; // Road in z for seeded cluster
1321
1322 const Double_t kRoad2y = 3.0; // Road in y for extrapolated cluster
1323 const Double_t kRoad2z = 20.0; // Road in z for extrapolated cluster
1324 const Int_t kMaxSeed = 3000;
1325
1326 Int_t maxSec = AliTRDgeometry::kNsector;
1327
1328 // Linear fitters in planes
1329 TLinearFitter fitterTC(2,"hyp2"); // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1330 TLinearFitter fitterT2(4,"hyp4"); // Fitting with tilting pads - kz not fixed
1331 fitterTC.StoreData(kTRUE);
1332 fitterT2.StoreData(kTRUE);
1333 AliRieman rieman(1000); // Rieman fitter
1334 AliRieman rieman2(1000); // Rieman fitter
1335
1336 // Find the maximal and minimal layer for the planes
1337 Int_t layers[6][2];
1338 AliTRDpropagationLayer *reflayers[6];
1339 for (Int_t i = 0; i < 6; i++) {
1340 layers[i][0] = 10000;
1341 layers[i][1] = 0;
1342 }
1343 for (Int_t ns = 0; ns < maxSec; ns++) {
1344 for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1345 AliTRDpropagationLayer &propLayer = *(fTrSec[ns]->GetLayer(ilayer));
1346 if (propLayer == 0) {
1347 continue;
1348 }
1349 Int_t det = propLayer[0]->GetDetector();
1350 Int_t layer = fGeom->GetLayer(det);
1351 if (ilayer < layers[layer][0]) {
1352 layers[layer][0] = ilayer;
1353 }
1354 if (ilayer > layers[layer][1]) {
1355 layers[layer][1] = ilayer;
1356 }
1357 }
1358 }
1359
1360 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0);
1361 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
1362 Double_t hL[6]; // Tilting angle
1363 Double_t xcl[6]; // X - position of reference cluster
1364 Double_t ycl[6]; // Y - position of reference cluster
1365 Double_t zcl[6]; // Z - position of reference cluster
1366
1367 AliTRDcluster *cl[6] = { 0, 0, 0, 0, 0, 0 }; // Seeding clusters
1368 Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1369
1370 Double_t chi2R = 0.0;
1371 Double_t chi2Z = 0.0;
1372 Double_t chi2RF = 0.0;
1373 Double_t chi2ZF = 0.0;
1374
1375 Int_t nclusters; // Total number of clusters
1376 for (Int_t i = 0; i < 6; i++) {
1377 hL[i] = h01;
1378 if (i%2==1) {
1379 hL[i]*=-1.0;
1380 }
1381 }
1382
1383 // Registered seed
1384 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1385 AliTRDseed *seed[kMaxSeed];
1386 for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1387 seed[iseed]= &pseed[iseed*6];
1388 }
1389 AliTRDseed *cseed = seed[0];
1390
1391 Double_t seedquality[kMaxSeed];
1392 Double_t seedquality2[kMaxSeed];
1393 Double_t seedparams[kMaxSeed][7];
1394 Int_t seedlayer[kMaxSeed];
1395 Int_t registered = 0;
1396 Int_t sort[kMaxSeed];
1397
1398 //
1399 // Seeding part
1400 //
1401 for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1402 //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1403
1404 registered = 0; // Reset registerd seed counter
1405 cseed = seed[registered];
1406 Float_t iter = 0.0;
1407
1408 for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1409 //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1410
1411 iter += 1.0;
1412 Int_t dseed = 5 + Int_t(iter) * 3;
1413
1414 // Initialize seeding layers
1415 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
1416 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1417 xcl[ilayer] = reflayers[ilayer]->GetX();
1418 }
1419
1420 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1421 AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1422 AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1423 AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1424 AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1425
1426 Int_t maxn3 = layer3;
1427 for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
1428
1429 AliTRDcluster *cl3 = layer3[icl3];
1430 if (!cl3) {
1431 continue;
1432 }
1433 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
1434 ycl[sLayer+3] = cl3->GetY();
1435 zcl[sLayer+3] = cl3->GetZ();
1436 Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1437 Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1438 Int_t maxn0 = layer0;
1439
1440 for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1441
1442 AliTRDcluster *cl0 = layer0[icl0];
1443 if (!cl0) {
1444 continue;
1445 }
1446 if (cl3->IsUsed() && cl0->IsUsed()) {
1447 continue;
1448 }
1449 ycl[sLayer+0] = cl0->GetY();
1450 zcl[sLayer+0] = cl0->GetZ();
1451 if (ycl[sLayer+0] > yymax0) {
1452 break;
1453 }
1454 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1455 if (TMath::Abs(tanphi) > kMaxPhi) {
1456 continue;
1457 }
1458 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1459 if (TMath::Abs(tantheta) > kMaxTheta) {
1460 continue;
1461 }
1462 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1463
1464 // Expected position in 1 layer
1465 Double_t y1exp = ycl[sLayer+0] + (tanphi) * (xcl[sLayer+1]-xcl[sLayer+0]);
1466 Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1467 Float_t yymin1 = y1exp - kRoad0y - tanphi;
1468 Float_t yymax1 = y1exp + kRoad0y + tanphi;
1469 Int_t maxn1 = layer1;
1470
1471 for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1472
1473 AliTRDcluster *cl1 = layer1[icl1];
1474 if (!cl1) {
1475 continue;
1476 }
1477 Int_t nusedCl = 0;
1478 if (cl3->IsUsed()) nusedCl++;
1479 if (cl0->IsUsed()) nusedCl++;
1480 if (cl1->IsUsed()) nusedCl++;
1481 if (nusedCl > 1) {
1482 continue;
1483 }
1484 ycl[sLayer+1] = cl1->GetY();
1485 zcl[sLayer+1] = cl1->GetZ();
1486 if (ycl[sLayer+1] > yymax1) {
1487 break;
1488 }
1489 if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1490 continue;
1491 }
1492 if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1493 continue;
1494 }
1495 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1496
1497 Double_t y2exp = ycl[sLayer+0]+(tanphi) * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1498 Double_t z2exp = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1499 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1500 if (index2 <= 0) {
1501 continue;
1502 }
1503 AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1504 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1505 ycl[sLayer+2] = cl2->GetY();
1506 zcl[sLayer+2] = cl2->GetZ();
1507 if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1508 continue;
1509 }
1510
1511 rieman.Reset();
1512 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1513 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1514 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1515 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1516 rieman.Update();
1517
1518 // Reset fitter
1519 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1520 cseed[iLayer].Reset();
1521 }
1522 chi2Z = 0.0;
1523 chi2R = 0.0;
1524
1525 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1526 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1527 chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1528 * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1529 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1530 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1531 chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1532 * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1533 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1534 }
1535 if (TMath::Sqrt(chi2R) > 1.0/iter) {
1536 continue;
1537 }
1538 if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1539 continue;
1540 }
1541
1542 Float_t minmax[2] = { -100.0, 100.0 };
1543 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1544 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
1545 + 1.0 - cseed[sLayer+iLayer].GetZref(0);
1546 if (max < minmax[1]) {
1547 minmax[1] = max;
1548 }
1549 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
1550 - 1.0 - cseed[sLayer+iLayer].GetZref(0);
1551 if (min > minmax[0]) {
1552 minmax[0] = min;
1553 }
1554 }
1555
1556 Bool_t isFake = kFALSE;
1557 if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1558 isFake = kTRUE;
1559 }
1560 if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1561 isFake = kTRUE;
1562 }
1563 if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1564 isFake = kTRUE;
1565 }
1566
1567 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
1568 if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
1569 TTreeSRedirector &cstream = *fDebugStreamer;
1570 cstream << "Seeds0"
1571 << "isFake=" << isFake
1572 << "Cl0.=" << cl0
1573 << "Cl1.=" << cl1
1574 << "Cl2.=" << cl2
1575 << "Cl3.=" << cl3
1576 << "Xref=" << xref
1577 << "X0=" << xcl[sLayer+0]
1578 << "X1=" << xcl[sLayer+1]
1579 << "X2=" << xcl[sLayer+2]
1580 << "X3=" << xcl[sLayer+3]
1581 << "Y2exp=" << y2exp
1582 << "Z2exp=" << z2exp
1583 << "Chi2R=" << chi2R
1584 << "Chi2Z=" << chi2Z
1585 << "Seed0.=" << &cseed[sLayer+0]
1586 << "Seed1.=" << &cseed[sLayer+1]
1587 << "Seed2.=" << &cseed[sLayer+2]
1588 << "Seed3.=" << &cseed[sLayer+3]
1589 << "Zmin=" << minmax[0]
1590 << "Zmax=" << minmax[1]
1591 << "\n";
1592 }
1593 }
1594
1595 ////////////////////////////////////////////////////////////////////////////////////
1596 //
1597 // Fit seeding part
1598 //
1599 ////////////////////////////////////////////////////////////////////////////////////
1600
1601 cl[sLayer+0] = cl0;
1602 cl[sLayer+1] = cl1;
1603 cl[sLayer+2] = cl2;
1604 cl[sLayer+3] = cl3;
1605 Bool_t isOK = kTRUE;
1606
1607 for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1608
1609 cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1610 cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1611 cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
1612
1613 for (Int_t jter = 0; jter < 2; jter++) {
1614
1615 //
1616 // In iteration 0 we try only one pad-row
1617 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
1618 //
1619 AliTRDseed tseed = cseed[sLayer+jLayer];
1620 Float_t roadz = padlength[sLayer+jLayer] * 0.5;
1621 if (jter > 0) {
1622 roadz = padlength[sLayer+jLayer];
1623 }
1624
1625 Float_t quality = 10000.0;
1626
1627 for (Int_t iTime = 2; iTime < 20; iTime++) {
1628
1629 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1630 Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1631 Double_t zexp = cl[sLayer+jLayer]->GetZ();
1632
1633 if (jter > 0) {
1634 // Try 2 pad-rows in second iteration
1635 zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
1636 if (zexp > cl[sLayer+jLayer]->GetZ()) {
1637 zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1638 }
1639 if (zexp < cl[sLayer+jLayer]->GetZ()) {
1640 zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1641 }
1642 }
1643
1644 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
1645 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1646 if (index <= 0) {
1647 continue;
1648 }
1649 AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);
1650
1651 tseed.SetIndexes(iTime,index);
1652 tseed.SetClusters(iTime,clu); // Register cluster
1653 tseed.SetX(iTime,dxlayer); // Register cluster
1654 tseed.SetY(iTime,clu->GetY()); // Register cluster
1655 tseed.SetZ(iTime,clu->GetZ()); // Register cluster
1656
1657 }
1658
1659 tseed.Update();
1660
1661 // Count the number of clusters and distortions into quality
1662 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1663 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1664 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1665 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1666 if ((jter == 0) && tseed.IsOK()) {
1667 cseed[sLayer+jLayer] = tseed;
1668 quality = tquality;
1669 if (tquality < 5) {
1670 break;
1671 }
1672 }
1673 if (tseed.IsOK() && (tquality < quality)) {
1674 cseed[sLayer+jLayer] = tseed;
1675 }
1676
1677 } // Loop: jter
1678
1679 if (!cseed[sLayer+jLayer].IsOK()) {
1680 isOK = kFALSE;
1681 break;
1682 }
1683
1684 cseed[sLayer+jLayer].CookLabels();
1685 cseed[sLayer+jLayer].UpdateUsed();
1686 nusedCl += cseed[sLayer+jLayer].GetNUsed();
1687 if (nusedCl > 25) {
1688 isOK = kFALSE;
1689 break;
1690 }
1691
1692 } // Loop: jLayer
1693
1694 if (!isOK) {
1695 continue;
1696 }
1697 nclusters = 0;
1698 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1699 if (cseed[sLayer+iLayer].IsOK()) {
1700 nclusters += cseed[sLayer+iLayer].GetN2();
1701 }
1702 }
1703
1704 // Iteration 0
1705 rieman.Reset();
1706 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1707 rieman.AddPoint(xcl[sLayer+iLayer]
1708 ,cseed[sLayer+iLayer].GetYfitR(0)
1709 ,cseed[sLayer+iLayer].GetZProb()
1710 ,1
1711 ,10);
1712 }
1713 rieman.Update();
1714
1715 chi2R = 0.0;
1716 chi2Z = 0.0;
1717
1718 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1719 cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1720 chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1721 * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1722 cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1723 cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1724 chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1725 * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1726 cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1727 }
1728 Double_t curv = rieman.GetC();
1729
1730 //
1731 // Likelihoods
1732 //
1733 Double_t sumda = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1734 + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1735 + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1736 + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
1737 Double_t likea = TMath::Exp(-sumda*10.6);
1738 Double_t likechi2 = 0.0000000001;
1739 if (chi2R < 0.5) {
1740 likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1741 }
1742 Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1743 Double_t likeN = TMath::Exp(-(72 - nclusters) * 0.19);
1744 Double_t like = likea * likechi2 * likechi2z * likeN;
1745 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1746 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1747 - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
1748 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
1749
1750 seedquality[registered] = like;
1751 seedlayer[registered] = sLayer;
1752 if (TMath::Log(0.000000000000001 + like) < -15) {
1753 continue;
1754 }
1755 AliTRDseed seedb[6];
1756 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1757 seedb[iLayer] = cseed[iLayer];
1758 }
1759
1760 ////////////////////////////////////////////////////////////////////////////////////
1761 //
1762 // Full track fit part
1763 //
1764 ////////////////////////////////////////////////////////////////////////////////////
1765
1766 Int_t nlayers = 0;
1767 Int_t nusedf = 0;
1768 Int_t findable = 0;
1769
1770 //
1771 // Add new layers - avoid long extrapolation
1772 //
1773 Int_t tLayer[2] = { 0, 0 };
1774 if (sLayer == 2) {
1775 tLayer[0] = 1;
1776 tLayer[1] = 0;
1777 }
1778 if (sLayer == 1) {
1779 tLayer[0] = 5;
1780 tLayer[1] = 0;
1781 }
1782 if (sLayer == 0) {
1783 tLayer[0] = 4;
1784 tLayer[1] = 5;
1785 }
1786
1787 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1788 Int_t jLayer = tLayer[iLayer]; // Set tracking layer
1789 cseed[jLayer].Reset();
1790 cseed[jLayer].SetTilt(hL[jLayer]);
1791 cseed[jLayer].SetPadLength(padlength[jLayer]);
1792 cseed[jLayer].SetX0(xcl[jLayer]);
1793 // Get pad length and rough cluster
1794 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
1795 ,cseed[jLayer].GetZref(0)
1796 ,kRoad2y
1797 ,kRoad2z);
1798 if (indexdummy <= 0) {
1799 continue;
1800 }
1801 AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
1802 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
1803 }
1804 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1805
1806 for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1807
1808 Int_t jLayer = tLayer[iLayer]; // set tracking layer
1809 if ((jLayer == 0) && !(cseed[1].IsOK())) {
1810 continue; // break not allowed
1811 }
1812 if ((jLayer == 5) && !(cseed[4].IsOK())) {
1813 continue; // break not allowed
1814 }
1815 Float_t zexp = cseed[jLayer].GetZref(0);
1816 Double_t zroad = padlength[jLayer] * 0.5 + 1.0;
1817
1818 for (Int_t jter = 0; jter < 2; jter++) {
1819
1820 AliTRDseed tseed = cseed[jLayer];
1821 Float_t quality = 10000.0;
1822
1823 for (Int_t iTime = 2; iTime < 20; iTime++) {
1824 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1825 Double_t dxlayer = layer.GetX()-xcl[jLayer];
1826 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
1827 Float_t yroad = kRoad1y;
1828 Int_t index = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
1829 if (index <= 0) {
1830 continue;
1831 }
1832 AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);
1833 tseed.SetIndexes(iTime,index);
1834 tseed.SetClusters(iTime,clu); // Register cluster
1835 tseed.SetX(iTime,dxlayer); // Register cluster
1836 tseed.SetY(iTime,clu->GetY()); // Register cluster
1837 tseed.SetZ(iTime,clu->GetZ()); // Register cluster
1838 }
1839
1840 tseed.Update();
1841 if (tseed.IsOK()) {
1842 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1843 Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
1844 + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1845 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1846 if (tquality < quality) {
1847 cseed[jLayer] = tseed;
1848 quality = tquality;
1849 }
1850 }
1851
1852 zroad *= 2.0;
1853
1854 } // Loop: jter
1855
1856 if ( cseed[jLayer].IsOK()) {
1857 cseed[jLayer].CookLabels();
1858 cseed[jLayer].UpdateUsed();
1859 nusedf += cseed[jLayer].GetNUsed();
1860 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1861 }
1862
1863 } // Loop: iLayer
1864
1865 // Make copy
1866 AliTRDseed bseed[6];
1867 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1868 bseed[jLayer] = cseed[jLayer];
1869 }
1870 Float_t lastquality = 10000.0;
1871 Float_t lastchi2 = 10000.0;
1872 Float_t chi2 = 1000.0;
1873
1874 for (Int_t jter = 0; jter < 4; jter++) {
1875
1876 // Sort tracklets according "quality", try to "improve" 4 worst
1877 Float_t sumquality = 0.0;
1878 Float_t squality[6];
1879 Int_t sortindexes[6];
1880
1881 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1882 if (bseed[jLayer].IsOK()) {
1883 AliTRDseed &tseed = bseed[jLayer];
1884 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1885 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1886 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1887 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1888 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1889 squality[jLayer] = tquality;
1890 }
1891 else {
1892 squality[jLayer] = -1.0;
1893 }
1894 sumquality +=squality[jLayer];
1895 }
1896
1897 if ((sumquality >= lastquality) ||
1898 (chi2 > lastchi2)) {
1899 break;
1900 }
1901 lastquality = sumquality;
1902 lastchi2 = chi2;
1903 if (jter > 0) {
1904 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1905 cseed[jLayer] = bseed[jLayer];
1906 }
1907 }
1908 TMath::Sort(6,squality,sortindexes,kFALSE);
1909
1910 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1911
1912 Int_t bLayer = sortindexes[jLayer];
1913 AliTRDseed tseed = bseed[bLayer];
1914
1915 for (Int_t iTime = 2; iTime < 20; iTime++) {
1916
1917 AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1918 Double_t dxlayer = layer.GetX() - xcl[bLayer];
1919 Double_t zexp = tseed.GetZref(0);
1920 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1921 Float_t roadz = padlength[bLayer] + 1;
1922 if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
1923 roadz = padlength[bLayer] * 0.5;
1924 }
1925 if (tseed.GetZfit(1)*tseed.GetZref(1) < 0.0) {
1926 roadz = padlength[bLayer] * 0.5;
1927 }
1928 if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
1929 zexp = tseed.GetZProb();
1930 roadz = padlength[bLayer] * 0.5;
1931 }
1932
1933 Double_t yexp = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
1934 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1935 if (index <= 0) {
1936 continue;
1937 }
1938 AliTRDcluster *clu = (AliTRDcluster *) GetCluster(index);
1939
1940 tseed.SetIndexes(iTime,index);
1941 tseed.SetClusters(iTime,clu); // Register cluster
1942 tseed.SetX(iTime,dxlayer); // Register cluster
1943 tseed.SetY(iTime,clu->GetY()); // Register cluster
1944 tseed.SetZ(iTime,clu->GetZ()); // Register cluster
1945
1946 }
1947
1948 tseed.Update();
1949 if (tseed.IsOK()) {
1950 Float_t dangle = tseed.GetYfit(1) - tseed.GetYref(1);
1951 Double_t zcor = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1952 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0
1953 + TMath::Abs(dangle) / 0.1
1954 + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1955 + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1956 if (tquality<squality[bLayer]) {
1957 bseed[bLayer] = tseed;
1958 }
1959 }
1960
1961 } // Loop: jLayer
1962
1963 chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
1964
1965 } // Loop: jter
1966
1967 nclusters = 0;
1968 nlayers = 0;
1969 findable = 0;
1970 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1971 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
1972 findable++;
1973 }
1974 if (cseed[iLayer].IsOK()) {
1975 nclusters += cseed[iLayer].GetN2();
1976 nlayers++;
1977 }
1978 }
1979 if (nlayers < 3) {
1980 continue;
1981 }
1982 rieman.Reset();
1983 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1984 if (cseed[iLayer].IsOK()) {
1985 rieman.AddPoint(xcl[iLayer]
1986 ,cseed[iLayer].GetYfitR(0)
1987 ,cseed[iLayer].GetZProb()
1988 ,1
1989 ,10);
1990 }
1991 }
1992 rieman.Update();
1993
1994 chi2RF = 0.0;
1995 chi2ZF = 0.0;
1996 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1997 if (cseed[iLayer].IsOK()) {
1998 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
1999 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2000 * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2001 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2002 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2003 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2004 * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2005 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
2006 }
2007 }
2008 chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2009 chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2010 curv = rieman.GetC();
2011
2012 Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2013 Double_t dzmf = rieman.GetDZat(xref2);
2014 Double_t zmf = rieman.GetZat(xref2);
2015
2016 //
2017 // Fit hyperplane
2018 //
2019 Int_t npointsT = 0;
2020 fitterTC.ClearPoints();
2021 fitterT2.ClearPoints();
2022 rieman2.Reset();
2023
2024 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2025
2026 if (!cseed[iLayer].IsOK()) {
2027 continue;
2028 }
2029
2030 for (Int_t itime = 0; itime < 25; itime++) {
2031
2032 if (!cseed[iLayer].IsUsable(itime)) {
2033 continue;
2034 }
2035 // X relative to the middle chamber
2036 Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2037 Double_t y = cseed[iLayer].GetY(itime);
2038 Double_t z = cseed[iLayer].GetZ(itime);
2039 // ExB correction to the correction
2040 // Tilted rieman
2041 Double_t uvt[6];
2042 // Global x
2043 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
2044 Double_t t = 1.0 / (x2*x2 + y*y);
2045 uvt[1] = t; // t
2046 uvt[0] = 2.0 * x2 * uvt[1]; // u
2047 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2048 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2049 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2050 Double_t error = 2.0 * 0.2 * uvt[1];
2051 fitterT2.AddPoint(uvt,uvt[4],error);
2052
2053 //
2054 // Constrained rieman
2055 //
2056 z = cseed[iLayer].GetZ(itime);
2057 uvt[0] = 2.0 * x2 * t; // u
2058 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2059 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
2060 fitterTC.AddPoint(uvt,uvt[2],error);
2061 rieman2.AddPoint(x2,y,z,1,10);
2062 npointsT++;
2063
2064 }
2065
2066 } // Loop: iLayer
2067
2068 rieman2.Update();
2069 fitterTC.Eval();
2070 fitterT2.Eval();
2071 Double_t rpolz0 = fitterT2.GetParameter(3);
2072 Double_t rpolz1 = fitterT2.GetParameter(4);
2073
2074 //
2075 // Linear fitter - not possible to make boundaries
2076 // Do not accept non possible z and dzdx combinations
2077 //
2078 Bool_t acceptablez = kTRUE;
2079 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2080 if (cseed[iLayer].IsOK()) {
2081 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2082 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
2083 acceptablez = kFALSE;
2084 }
2085 }
2086 }
2087 if (!acceptablez) {
2088 fitterT2.FixParameter(3,zmf);
2089 fitterT2.FixParameter(4,dzmf);
2090 fitterT2.Eval();
2091 fitterT2.ReleaseParameter(3);
2092 fitterT2.ReleaseParameter(4);
2093 rpolz0 = fitterT2.GetParameter(3);
2094 rpolz1 = fitterT2.GetParameter(4);
2095 }
2096
2097 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2098 Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
2099 Double_t polz1c = fitterTC.GetParameter(2);
2100 Double_t polz0c = polz1c * xref2;
2101 Double_t aC = fitterTC.GetParameter(0);
2102 Double_t bC = fitterTC.GetParameter(1);
2103 Double_t cC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
2104 Double_t aR = fitterT2.GetParameter(0);
2105 Double_t bR = fitterT2.GetParameter(1);
2106 Double_t dR = fitterT2.GetParameter(2);
2107 Double_t cR = 1.0 + bR*bR - dR*aR;
2108 Double_t dca = 0.0;
2109 if (cR > 0.0) {
2110 dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2111 cR = aR / TMath::Sqrt(cR);
2112 }
2113
2114 Double_t chi2ZT2 = 0.0;
2115 Double_t chi2ZTC = 0.0;
2116 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2117 if (cseed[iLayer].IsOK()) {
2118 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2119 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
2120 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2121 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
2122 }
2123 }
2124 chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2125 chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2126
2127 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2128 Float_t sumdaf = 0.0;
2129 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2130 if (cseed[iLayer].IsOK()) {
2131 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2132 / cseed[iLayer].GetSigmaY2());
2133 }
2134 }
2135 sumdaf /= Float_t (nlayers - 2.0);
2136
2137 //
2138 // Likelihoods for full track
2139 //
2140 Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
2141 Double_t likechi2C = TMath::Exp(-chi2TC * 0.677);
2142 Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2143 Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
2144 seedquality2[registered] = likezf * likechi2TR * likeaf;
2145
2146 // Still needed ????
2147// Bool_t isGold = kFALSE;
2148//
2149// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2150// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
2151// if (isGold &&nusedf<10){
2152// for (Int_t jLayer=0;jLayer<6;jLayer++){
2153// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
2154// seed[index][jLayer].UseClusters(); //sign gold
2155// }
2156// }
2157
2158 Int_t index0 = 0;
2159 if (!cseed[0].IsOK()) {
2160 index0 = 1;
2161 if (!cseed[1].IsOK()) {
2162 index0 = 2;
2163 }
2164 }
2165 seedparams[registered][0] = cseed[index0].GetX0();
2166 seedparams[registered][1] = cseed[index0].GetYref(0);
2167 seedparams[registered][2] = cseed[index0].GetZref(0);
2168 seedparams[registered][5] = cR;
2169 seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2170 seedparams[registered][4] = cseed[index0].GetZref(1)
2171 / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
2172 seedparams[registered][6] = ns;
2173
2174 Int_t labels[12];
2175 Int_t outlab[24];
2176 Int_t nlab = 0;
2177 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2178 if (!cseed[iLayer].IsOK()) {
2179 continue;
2180 }
2181 if (cseed[iLayer].GetLabels(0) >= 0) {
2182 labels[nlab] = cseed[iLayer].GetLabels(0);
2183 nlab++;
2184 }
2185 if (cseed[iLayer].GetLabels(1) >= 0) {
2186 labels[nlab] = cseed[iLayer].GetLabels(1);
2187 nlab++;
2188 }
2189 }
2190 Freq(nlab,labels,outlab,kFALSE);
2191 Int_t label = outlab[0];
2192 Int_t frequency = outlab[1];
2193 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2194 cseed[iLayer].SetFreq(frequency);
2195 cseed[iLayer].SetC(cR);
2196 cseed[iLayer].SetCC(cC);
2197 cseed[iLayer].SetChi2(chi2TR);
2198 cseed[iLayer].SetChi2Z(chi2ZF);
2199 }
2200
2201 // Debugging print
2202 if (1 || (!isFake)) {
2203 Float_t zvertex = GetZ();
2204 TTreeSRedirector &cstream = *fDebugStreamer;
2205 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
2206 cstream << "Seeds1"
2207 << "isFake=" << isFake
2208 << "Vertex=" << zvertex
2209 << "Rieman2.=" << &rieman2
2210 << "Rieman.=" << &rieman
2211 << "Xref=" << xref
2212 << "X0=" << xcl[0]
2213 << "X1=" << xcl[1]
2214 << "X2=" << xcl[2]
2215 << "X3=" << xcl[3]
2216 << "X4=" << xcl[4]
2217 << "X5=" << xcl[5]
2218 << "Chi2R=" << chi2R
2219 << "Chi2Z=" << chi2Z
2220 << "Chi2RF=" << chi2RF // Chi2 of trackletes on full track
2221 << "Chi2ZF=" << chi2ZF // Chi2 z on tracklets on full track
2222 << "Chi2ZT2=" << chi2ZT2 // Chi2 z on tracklets on full track - rieman tilt
2223 << "Chi2ZTC=" << chi2ZTC // Chi2 z on tracklets on full track - rieman tilt const
2224 << "Chi2TR=" << chi2TR // Chi2 without vertex constrain
2225 << "Chi2TC=" << chi2TC // Chi2 with vertex constrain
2226 << "C=" << curv // Non constrained - no tilt correction
2227 << "DR=" << dR // DR parameter - tilt correction
2228 << "DCA=" << dca // DCA - tilt correction
2229 << "CR=" << cR // Non constrained curvature - tilt correction
2230 << "CC=" << cC // Constrained curvature
2231 << "Polz0=" << polz0c
2232 << "Polz1=" << polz1c
2233 << "RPolz0=" << rpolz0
2234 << "RPolz1=" << rpolz1
2235 << "Ncl=" << nclusters
2236 << "Nlayers=" << nlayers
2237 << "NUsedS=" << nusedCl
2238 << "NUsed=" << nusedf
2239 << "Findable=" << findable
2240 << "Like=" << like
2241 << "LikePrim=" << likePrim
2242 << "Likechi2C=" << likechi2C
2243 << "Likechi2TR=" << likechi2TR
2244 << "Likezf=" << likezf
2245 << "LikeF=" << seedquality2[registered]
2246 << "S0.=" << &cseed[0]
2247 << "S1.=" << &cseed[1]
2248 << "S2.=" << &cseed[2]
2249 << "S3.=" << &cseed[3]
2250 << "S4.=" << &cseed[4]
2251 << "S5.=" << &cseed[5]
2252 << "SB0.=" << &seedb[0]
2253 << "SB1.=" << &seedb[1]
2254 << "SB2.=" << &seedb[2]
2255 << "SB3.=" << &seedb[3]
2256 << "SB4.=" << &seedb[4]
2257 << "SB5.=" << &seedb[5]
2258 << "Label=" << label
2259 << "Freq=" << frequency
2260 << "sLayer=" << sLayer
2261 << "\n";
2262 }
2263 }
2264
2265 if (registered<kMaxSeed - 1) {
2266 registered++;
2267 cseed = seed[registered];
2268 }
2269
2270 } // End of loop over layer 1
2271
2272 } // End of loop over layer 0
2273
2274 } // End of loop over layer 3
2275
2276 } // End of loop over seeding time bins
2277
2278 //
2279 // Choose best
2280 //
2281
2282 TMath::Sort(registered,seedquality2,sort,kTRUE);
2283 Bool_t signedseed[kMaxSeed];
2284 for (Int_t i = 0; i < registered; i++) {
2285 signedseed[i] = kFALSE;
2286 }
2287
2288 for (Int_t jter = 0; jter < 5; jter++) {
2289
2290 for (Int_t iseed = 0; iseed < registered; iseed++) {
2291
2292 Int_t index = sort[iseed];
2293 if (signedseed[index]) {
2294 continue;
2295 }
2296 Int_t labelsall[1000];
2297 Int_t nlabelsall = 0;
2298 Int_t naccepted = 0;;
2299 Int_t sLayer = seedlayer[index];
2300 Int_t ncl = 0;
2301 Int_t nused = 0;
2302 Int_t nlayers = 0;
2303 Int_t findable = 0;
2304
2305 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2306
2307 if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
2308 findable++;
2309 }
2310 if (seed[index][jLayer].IsOK()) {
2311 seed[index][jLayer].UpdateUsed();
2312 ncl +=seed[index][jLayer].GetN2();
2313 nused +=seed[index][jLayer].GetNUsed();
2314 nlayers++;
2315 // Cooking label
2316 for (Int_t itime = 0; itime < 25; itime++) {
2317 if (seed[index][jLayer].IsUsable(itime)) {
2318 naccepted++;
2319 for (Int_t ilab = 0; ilab < 3; ilab++) {
2320 Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
2321 if (tindex >= 0) {
2322 labelsall[nlabelsall] = tindex;
2323 nlabelsall++;
2324 }
2325 }
2326 }
2327 }
2328 }
2329
2330 }
2331
2332 if (nused > 30) {
2333 continue;
2334 }
2335
2336 if (jter == 0) {
2337 if (nlayers < 6) {
2338 continue;
2339 }
2340 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2341 continue; // Gold
2342 }
2343 }
2344
2345 if (jter == 1) {
2346 if (nlayers < findable) {
2347 continue;
2348 }
2349 if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2350 continue;
2351 }
2352 }
2353
2354 if (jter == 2) {
2355 if ((nlayers == findable) ||
2356 (nlayers == 6)) {
2357 continue;
2358 }
2359 if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2360 continue;
2361 }
2362 }
2363
2364 if (jter == 3) {
2365 if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2366 continue;
2367 }
2368 }
2369
2370 if (jter == 4) {
2371 if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2372 continue;
2373 }
2374 }
2375
2376 signedseed[index] = kTRUE;
2377
2378 Int_t labels[1000];
2379 Int_t outlab[1000];
2380 Int_t nlab = 0;
2381 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2382 if (seed[index][iLayer].IsOK()) {
2383 if (seed[index][iLayer].GetLabels(0) >= 0) {
2384 labels[nlab] = seed[index][iLayer].GetLabels(0);
2385 nlab++;
2386 }
2387 if (seed[index][iLayer].GetLabels(1) >= 0) {
2388 labels[nlab] = seed[index][iLayer].GetLabels(1);
2389 nlab++;
2390 }
2391 }
2392 }
2393 Freq(nlab,labels,outlab,kFALSE);
2394 Int_t label = outlab[0];
2395 Int_t frequency = outlab[1];
2396 Freq(nlabelsall,labelsall,outlab,kFALSE);
2397 Int_t label1 = outlab[0];
2398 Int_t label2 = outlab[2];
2399 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2400 Float_t ratio = Float_t(nused) / Float_t(ncl);
2401 if (ratio < 0.25) {
2402 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2403 if ((seed[index][jLayer].IsOK()) &&
2404 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
2405 seed[index][jLayer].UseClusters(); // Sign gold
2406 }
2407 }
2408 }
2409
2410 Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
2411 TTreeSRedirector &cstream = *fDebugStreamer;
2412
2413 //
2414 // Register seed
2415 //
2416 AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2417 AliTRDtrack dummy;
2418 if (!track) {
2419 track = &dummy;
2420 }
2421 else {
2422 AliESDtrack esdtrack;
2423 esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
2424 esdtrack.SetLabel(label);
2425 esd->AddTrack(&esdtrack);
2426 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
2427 cstream << "Tracks"
2428 << "EventNrInFile=" << eventNrInFile
2429 << "ESD.=" << &esdtrack
2430 << "trd.=" << track
2431 << "trdback.=" << track
2432 << "\n";
2433 }
2434 }
2435
2436 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
2437 cstream << "Seeds2"
2438 << "Iter=" << jter
2439 << "Track.=" << track
2440 << "Like=" << seedquality[index]
2441 << "LikeF=" << seedquality2[index]
2442 << "S0.=" << &seed[index][0]
2443 << "S1.=" << &seed[index][1]
2444 << "S2.=" << &seed[index][2]
2445 << "S3.=" << &seed[index][3]
2446 << "S4.=" << &seed[index][4]
2447 << "S5.=" << &seed[index][5]
2448 << "Label=" << label
2449 << "Label1=" << label1
2450 << "Label2=" << label2
2451 << "FakeRatio=" << fakeratio
2452 << "Freq=" << frequency
2453 << "Ncl=" << ncl
2454 << "Nlayers=" << nlayers
2455 << "Findable=" << findable
2456 << "NUsed=" << nused
2457 << "sLayer=" << sLayer
2458 << "EventNrInFile=" << eventNrInFile
2459 << "\n";
2460 }
2461
2462 } // Loop: iseed
2463
2464 } // Loop: jter
2465
2466 } // End of loop over sectors
2467
2468 delete [] pseed;
2469
2470 return 0;
2471}
2472
2473//_____________________________________________________________________________
2474Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
2475{
2476 //
2477 // Reads AliTRDclusters from the file.
2478 // The names of the cluster tree and branches
2479 // should match the ones used in AliTRDclusterizer::WriteClusters()
2480 //
2481
2482 Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
2483 TObjArray *clusterArray = new TObjArray(nsize+1000);
2484
2485 TBranch *branch = clusterTree->GetBranch("TRDcluster");
2486 if (!branch) {
2487 AliError("Can't get the branch !");
2488 return 1;
2489 }
2490 branch->SetAddress(&clusterArray);
2491
2492 // Loop through all entries in the tree
2493 Int_t nEntries = (Int_t) clusterTree->GetEntries();
2494 Int_t nbytes = 0;
2495 AliTRDcluster *c = NULL;
2496 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2497
2498 // Import the tree
2499 nbytes += clusterTree->GetEvent(iEntry);
2500
2501 // Get the number of points in the detector
2502 Int_t nCluster = clusterArray->GetEntriesFast();
2503 // Loop through all TRD digits
2504 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2505 if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
2506 array->AddLast(c);
2507 //printf("Add cluster 0x%x.\n", c);
2508 clusterArray->RemoveAt(iCluster);
2509 }
2510
2511 }
2512 delete clusterArray;
2513
2514 return 0;
2515
2516}
2517
2518//_____________________________________________________________________________
2519Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
2520{
2521 //
2522 // Get track space point with index i
2523 // Origin: C.Cheshkov
2524 //
2525
2526 AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2527 Int_t idet = cl->GetDetector();
2528 Int_t isector = fGeom->GetSector(idet);
2529 Int_t istack = fGeom->GetStack(idet);
2530 Int_t ilayer = fGeom->GetLayer(idet);
2531 Double_t local[3];
2532 local[0] = GetX(isector,ilayer,cl->GetLocalTimeBin());
2533 local[1] = cl->GetY();
2534 local[2] = cl->GetZ();
2535 Double_t global[3];
2536 fGeom->RotateBack(idet,local,global);
2537 p.SetXYZ(global[0],global[1],global[2]);
2538 AliGeomManager::ELayerID iGeoLayer = AliGeomManager::kTRD1;
2539 switch (ilayer) {
2540 case 0:
2541 iGeoLayer = AliGeomManager::kTRD1;
2542 break;
2543 case 1:
2544 iGeoLayer = AliGeomManager::kTRD2;
2545 break;
2546 case 2:
2547 iGeoLayer = AliGeomManager::kTRD3;
2548 break;
2549 case 3:
2550 iGeoLayer = AliGeomManager::kTRD4;
2551 break;
2552 case 4:
2553 iGeoLayer = AliGeomManager::kTRD5;
2554 break;
2555 case 5:
2556 iGeoLayer = AliGeomManager::kTRD6;
2557 break;
2558 };
2559 Int_t modId = isector * fGeom->Nstack() + istack;
2560 UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,modId);
2561 p.SetVolumeID(volid);
2562
2563 return kTRUE;
2564
2565}
2566
2567//_____________________________________________________________________________
2568void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
2569{
2570 //
2571 // This cooks a label. Mmmmh, smells good...
2572 //
2573
2574 Int_t label = 123456789;
2575 Int_t index;
2576 Int_t i;
2577 Int_t j;
2578 Int_t ncl = pt->GetNumberOfClusters();
2579
2580 const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
2581
2582 Bool_t labelAdded;
2583
2584 Int_t **s = new Int_t* [kRange];
2585 for (i = 0; i < kRange; i++) {
2586 s[i] = new Int_t[2];
2587 }
2588 for (i = 0; i < kRange; i++) {
2589 s[i][0] = -1;
2590 s[i][1] = 0;
2591 }
2592
2593 Int_t t0;
2594 Int_t t1;
2595 Int_t t2;
2596
2597 for (i = 0; i < ncl; i++) {
2598 index = pt->GetClusterIndex(i);
2599 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2600 t0=c->GetLabel(0);
2601 t1=c->GetLabel(1);
2602 t2=c->GetLabel(2);
2603 }
2604
2605 for (i = 0; i < ncl; i++) {
2606 index = pt->GetClusterIndex(i);
2607 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2608 for (Int_t k = 0; k < 3; k++) {
2609 label = c->GetLabel(k);
2610 labelAdded = kFALSE;
2611 j = 0;
2612 if (label >= 0) {
2613 while ((!labelAdded) && (j < kRange)) {
2614 if ((s[j][0] == label) ||
2615 (s[j][1] == 0)) {
2616 s[j][0] = label;
2617 s[j][1] = s[j][1] + 1;
2618 labelAdded = kTRUE;
2619 }
2620 j++;
2621 }
2622 }
2623 }
2624 }
2625
2626 Int_t max = 0;
2627 label = -123456789;
2628
2629 for (i = 0; i < kRange; i++) {
2630 if (s[i][1] > max) {
2631 max = s[i][1];
2632 label = s[i][0];
2633 }
2634 }
2635
2636 for (i = 0; i < kRange; i++) {
2637 delete []s[i];
2638 }
2639
2640 delete []s;
2641
2642 if ((1.0 - Float_t(max)/ncl) > wrong) {
2643 label = -label;
2644 }
2645
2646 pt->SetLabel(label);
2647
2648}
2649
2650//_____________________________________________________________________________
2651void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
2652{
2653 //
2654 // Use clusters, but don't abuse them!
2655 //
2656
2657 const Float_t kmaxchi2 = 18;
2658 const Float_t kmincl = 10;
2659
2660 AliTRDtrack *track = (AliTRDtrack *) t;
2661
2662 Int_t ncl = t->GetNumberOfClusters();
2663 for (Int_t i = from; i < ncl; i++) {
2664 Int_t index = t->GetClusterIndex(i);
2665 AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2666 Int_t ilayer = fGeom->GetLayer(c->GetDetector());
2667 if (track->GetTracklets(ilayer).GetChi2() > kmaxchi2) {
2668 continue;
2669 }
2670 if (track->GetTracklets(ilayer).GetN() < kmincl) {
2671 continue;
2672 }
2673 if (!(c->IsUsed())) {
2674 c->Use();
2675 }
2676 }
2677
2678}
2679
2680//_____________________________________________________________________________
2681Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
2682{
2683 //
2684 // Parametrised "expected" error of the cluster reconstruction in Y
2685 //
2686
2687 Double_t s = 0.08 * 0.08;
2688 return s;
2689
2690}
2691
2692//_____________________________________________________________________________
2693Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
2694{
2695 //
2696 // Parametrised "expected" error of the cluster reconstruction in Z
2697 //
2698
2699 Double_t s = 9.0 * 9.0 / 12.0;
2700 return s;
2701
2702}
2703
2704//_____________________________________________________________________________
2705Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
2706{
2707 //
2708 // Returns radial position which corresponds to time bin <localTB>
2709 // in tracking sector <sector> and plane <plane>
2710 //
2711
2712 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
2713 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2714
2715 return fTrSec[sector]->GetLayer(pl)->GetX();
2716
2717}
2718
2719
2720//_____________________________________________________________________________
2721AliTRDtracker::AliTRDtrackingSector
2722 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2723 :fN(0)
2724 ,fGeom(geo)
2725 ,fGeomSector(gs)
2726{
2727 //
2728 // AliTRDtrackingSector Constructor
2729 //
2730
2731 AliTRDpadPlane *padPlane = 0;
2732 AliTRDpropagationLayer *ppl = 0;
2733
2734 // Get holes description from geometry
2735 Bool_t holes[AliTRDgeometry::kNstack];
2736 for (Int_t istack = 0; istack < AliTRDgeometry::kNstack; istack++) {
2737 holes[istack] = fGeom->IsHole(0,istack,gs);
2738 }
2739
2740 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
2741 fTimeBinIndex[i] = -1;
2742 }
2743
2744 Double_t x;
2745 Double_t dx;
2746 Double_t rho;
2747 Double_t radLength;
2748
2749 // Add layers for each of the planes
2750 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2751 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2752
2753 const Int_t kNstacks = AliTRDgeometry::Nstack();
2754 Int_t tbIndex;
2755 Double_t ymax = 0;
2756 Double_t ymaxsensitive = 0;
2757 Double_t *zc = new Double_t[kNstacks];
2758 Double_t *zmax = new Double_t[kNstacks];
2759 Double_t *zmaxsensitive = new Double_t[kNstacks];
2760
2761 for (Int_t layer = 0; layer < AliTRDgeometry::Nlayer(); layer++) {
2762
2763 ymax = fGeom->GetChamberWidth(layer) / 2.0;
2764 padPlane = fGeom->GetPadPlane(layer,0);
2765 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
2766
2767 for (Int_t st = 0; st < kNstacks; st++) {
2768 zmax[st] = fGeom->GetChamberLength(layer,st) / 2.0;
2769 Float_t pad = padPlane->GetRowSize(1);
2770 Float_t row0 = fGeom->GetRow0(layer,st,0);
2771 Int_t nPads = fGeom->GetRowMax(layer,st,0);
2772 zmaxsensitive[st] = Float_t(nPads) * pad / 2.0;
2773 zc[st] = -(pad * nPads) / 2.0 + row0;
2774 }
2775
2776 AliTRDcalibDB *fCalibration = AliTRDcalibDB::Instance();
2777 dx = fCalibration->GetVdrift(0,0,0)
2778 / AliTRDCommonParam::Instance()->GetSamplingFrequency();
2779 rho = 0.00295 * 0.85; //????
2780 radLength = 11.0;
2781
2782 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(layer);
2783 //Double_t xbottom = x0 - dxDrift;
2784 //Double_t xtop = x0 + dxAmp;
2785
2786 //temporary !! (A.Bercuci)
2787 Int_t t0 = (Int_t)fCalibration->GetT0Average(AliTRDgeometry::GetDetector(layer, 2, gs));
2788
2789 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
2790 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2791
2792 Double_t xlayer = iTime * dx - dxAmp;
2793 //if (xlayer<0) xlayer = dxAmp / 2.0;
2794 x = x0 - xlayer;
2795
2796 tbIndex = CookTimeBinIndex(layer,iTime);
2797 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,layer);
2798 ppl->SetYmax(ymax,ymaxsensitive);
2799 ppl->SetZ(zc,zmax,zmaxsensitive);
2800 ppl->SetHoles(holes);
2801 if(iTime == t0) ppl->SetT0();
2802
2803 InsertLayer(ppl);
2804
2805 }
2806
2807 }
2808
2809 MapTimeBinLayers();
2810
2811 delete [] zc;
2812 delete [] zmax;
2813 delete [] zmaxsensitive;
2814
2815}
2816
2817//_____________________________________________________________________________
2818AliTRDtracker::AliTRDtrackingSector
2819 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2820 :fN(0)
2821 ,fGeom(0)
2822 ,fGeomSector(0)
2823{
2824 //
2825 // Copy constructor
2826 //
2827
2828}
2829
2830//_____________________________________________________________________________
2831AliTRDtracker::AliTRDtrackingSector
2832 ::~AliTRDtrackingSector()
2833{
2834 //
2835 // Destructor
2836 //
2837
2838 for (Int_t i = 0; i < fN; i++) {
2839 delete fLayers[i];
2840 }
2841
2842}
2843
2844//_____________________________________________________________________________
2845Int_t AliTRDtracker::AliTRDtrackingSector
2846 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2847{
2848 //
2849 // Depending on the digitization parameters calculates global
2850 // (i.e. for a whole TRD stack of 6 planes) time bin index for
2851 // timebin <localTB> in plane <plane>
2852 //
2853
2854 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
2855 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1;
2856 if (localTB < 0) {
2857 return -1;
2858 }
2859 if (gtb < 0) {
2860 return -1;
2861 }
2862
2863 return gtb;
2864
2865}
2866
2867//_____________________________________________________________________________
2868void AliTRDtracker::AliTRDtrackingSector
2869 ::MapTimeBinLayers()
2870{
2871 //
2872 // For all sensitive time bins sets corresponding layer index
2873 // in the array fTimeBins
2874 //
2875
2876 Int_t index;
2877
2878 for (Int_t i = 0; i < fN; i++) {
2879
2880 index = fLayers[i]->GetTimeBinIndex();
2881
2882 if (index < 0) {
2883 continue;
2884 }
2885 if (index >= (Int_t) kMaxTimeBinIndex) {
2886 //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
2887 // ,index,kMaxTimeBinIndex-1));
2888 continue;
2889 }
2890
2891 fTimeBinIndex[index] = i;
2892
2893 }
2894
2895}
2896
2897//_____________________________________________________________________________
2898Int_t AliTRDtracker::AliTRDtrackingSector
2899 ::GetLayerNumber(Double_t x) const
2900{
2901 //
2902 // Returns the number of time bin which in radial position is closest to <x>
2903 //
2904
2905 if (x >= fLayers[fN-1]->GetX()) {
2906 return fN - 1;
2907 }
2908 if (x <= fLayers[ 0]->GetX()) {
2909 return 0;
2910 }
2911
2912 Int_t b = 0;
2913 Int_t e = fN - 1;
2914 Int_t m = (b + e) / 2;
2915
2916 for ( ; b < e; m = (b + e) / 2) {
2917 if (x > fLayers[m]->GetX()) {
2918 b = m + 1;
2919 }
2920 else {
2921 e = m;
2922 }
2923 }
2924
2925 if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
2926 return m + 1;
2927 }
2928 else {
2929 return m;
2930 }
2931
2932}
2933
2934//_____________________________________________________________________________
2935Int_t AliTRDtracker::AliTRDtrackingSector
2936 ::GetInnerTimeBin() const
2937{
2938 //
2939 // Returns number of the innermost SENSITIVE propagation layer
2940 //
2941
2942 return GetLayerNumber(0);
2943
2944}
2945
2946//_____________________________________________________________________________
2947Int_t AliTRDtracker::AliTRDtrackingSector
2948 ::GetOuterTimeBin() const
2949{
2950 //
2951 // Returns number of the outermost SENSITIVE time bin
2952 //
2953
2954 return GetLayerNumber(GetNumberOfTimeBins() - 1);
2955
2956}
2957
2958//_____________________________________________________________________________
2959Int_t AliTRDtracker::AliTRDtrackingSector
2960 ::GetNumberOfTimeBins() const
2961{
2962 //
2963 // Returns number of SENSITIVE time bins
2964 //
2965
2966 Int_t tb;
2967 Int_t layer;
2968
2969 for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
2970 layer = GetLayerNumber(tb);
2971 if (layer >= 0) {
2972 break;
2973 }
2974 }
2975
2976 return tb + 1;
2977
2978}
2979
2980//_____________________________________________________________________________
2981void AliTRDtracker::AliTRDtrackingSector
2982 ::InsertLayer(AliTRDpropagationLayer *pl)
2983{
2984 //
2985 // Insert layer <pl> in fLayers array.
2986 // Layers are sorted according to X coordinate.
2987 //
2988
2989 if (fN == ((Int_t) kMaxLayersPerSector)) {
2990 //AliWarning("Too many layers !\n");
2991 return;
2992 }
2993
2994 if (fN == 0) {
2995 fLayers[fN++] = pl;
2996 return;
2997 }
2998
2999 Int_t i = Find(pl->GetX());
3000
3001 memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3002
3003 fLayers[i] = pl;
3004 fN++;
3005
3006}
3007
3008//_____________________________________________________________________________
3009Int_t AliTRDtracker::AliTRDtrackingSector
3010 ::Find(Double_t x) const
3011{
3012 //
3013 // Returns index of the propagation layer nearest to X
3014 //
3015
3016 if (x <= fLayers[0]->GetX()) {
3017 return 0;
3018 }
3019
3020 if (x > fLayers[fN-1]->GetX()) {
3021 return fN;
3022 }
3023
3024 Int_t b = 0;
3025 Int_t e = fN-1;
3026 Int_t m = (b + e) / 2;
3027
3028 for (; b < e; m = (b + e) / 2) {
3029 if (x > fLayers[m]->GetX()) {
3030 b = m + 1;
3031 }
3032 else {
3033 e = m;
3034 }
3035 }
3036
3037 return m;
3038
3039}
3040
3041//_____________________________________________________________________________
3042Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
3043{
3044 //
3045 // Returns correction factor for tilted pads geometry
3046 //
3047
3048 Int_t det = c->GetDetector();
3049 Int_t layer = fGeom->GetLayer(det);
3050 AliTRDpadPlane *padPlane = fGeom->GetPadPlane(layer,0);
3051 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
3052
3053 if (fNoTilt) {
3054 h01 = 0;
3055 }
3056
3057 return h01;
3058
3059}
3060
3061//_____________________________________________________________________________
3062Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
3063 , AliTRDtrack *track
3064 , Int_t *clusters
3065 , AliTRDtracklet &tracklet)
3066{
3067 //
3068 // Try to find the nearest clusters to the track in the time bins
3069 // between <t0> and <t1>.
3070 // Also the corresponding tracklet is calculated
3071 // Correction coeficients - depend on TRD parameters - to be changed accordingly
3072 //
3073
3074 const Int_t kN1 = 100;
3075 const Int_t kN2 = 10;
3076
3077 Double_t x[kN1];
3078 Double_t yt[kN1];
3079 Double_t zt[kN1];
3080 Float_t zmean[kN1];
3081 Float_t nmean[kN1];
3082
3083 Double_t dz[kN2][kN1];
3084 Double_t dy[kN2][kN1];
3085 Int_t indexes[kN2][kN1]; // Indexes of the clusters in the road
3086 Int_t best[kN2][kN1]; // Index of best matching cluster
3087 AliTRDcluster *cl[kN2][kN1]; // Pointers to the clusters in the road
3088
3089 Double_t xmean = 0.0; // Reference x
3090 Int_t clfound = 0;
3091
3092 // Initialize the arrays
3093 for (Int_t it = 0; it < kN1; it++) {
3094
3095 x[it] = 0.0;
3096 yt[it] = 0.0;
3097 zt[it] = 0.0;
3098 clusters[it] = -2;
3099 zmean[it] = 0.0;
3100 nmean[it] = 0.0;
3101
3102 for (Int_t ih = 0; ih < kN2; ih++) {
3103 indexes[ih][it] = -2;
3104 cl[ih][it] = 0;
3105 dz[ih][it] = -100.0;
3106 dy[ih][it] = -100.0;
3107 best[ih][it] = 0;
3108 }
3109
3110 }
3111
3112 Double_t x0 = track->GetX();
3113 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3114 Int_t nall = 0;
3115 Int_t nfound = 0;
3116 Double_t h01 = 0.0;
3117 Int_t layer = -1;
3118 Int_t detector = -1;
3119 Float_t padlength = 0.0;
3120
3121 AliTRDtrack track2(* track);
3122 Float_t snpy = track->GetSnp();
3123 Float_t tany = TMath::Sqrt(snpy*snpy / ((1.-snpy)*(1.+snpy)));
3124 if (snpy < 0.0) {
3125 tany *= -1.0;
3126 }
3127
3128 Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetSignedPt());
3129 Double_t sz2 = ExpectedSigmaZ2(x0,track->GetTgl());
3130 Double_t road = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3131 if (road > 6.0) {
3132 road = 6.0;
3133 }
3134
3135 for (Int_t it = 0; it < t1-t0; it++) {
3136
3137 Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3138 AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3139 if (timeBin == 0) {
3140 continue; // No indexes1
3141 }
3142
3143 Int_t maxn = timeBin;
3144 x[it] = timeBin.GetX();
3145 track2.PropagateTo(x[it]);
3146 yt[it] = track2.GetY();
3147 zt[it] = track2.GetZ();
3148
3149 Double_t y = yt[it];
3150 Double_t z = zt[it];
3151 Double_t chi2 = 1000000.0;
3152 nall++;
3153
3154 //
3155 // Find 2 nearest cluster at given time bin
3156 //
3157 Int_t checkPoint[4] = { 0, 0, 0, 0 };
3158 Double_t minY = 123456789.0;
3159 Double_t minD[2] = { 1.0, 1.0 };
3160
3161 for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
3162
3163 AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
3164 h01 = GetTiltFactor(c);
3165 if (layer < 0) {
3166 Int_t det = c->GetDetector();
3167 layer = fGeom->GetLayer(det);
3168 padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3169 }
3170
3171 if (c->GetY() > (y + road)) {
3172 break;
3173 }
3174
3175 fHDeltaX->Fill(c->GetX() - x[it]);
3176
3177 if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) {
3178 minY = c->GetY() - y;
3179 minD[0] = c->GetY() - y;
3180 minD[1] = c->GetZ() - z;
3181 }
3182
3183 checkPoint[0]++;
3184
3185 fHMinZ->Fill(c->GetZ() - z);
3186 if ((c->GetZ() - z) * (c->GetZ() - z) > 2 * (12.0 * sz2)) {
3187 continue;
3188 }
3189 checkPoint[1]++;
3190
3191 Double_t dist = TMath::Abs(c->GetZ() - z);
3192 if (dist > (0.5 * padlength + 6.0 * sigmaz)) {
3193 continue; // 6 sigma boundary cut
3194 }
3195 checkPoint[2]++;
3196
3197 // Sigma boundary cost function
3198 Double_t cost = 0.0;
3199 if (dist> (0.5 * padlength - sigmaz)){
3200 cost = (dist - 0.5*padlength) / (2.0 * sigmaz);
3201 if (cost > -1) {
3202 cost = (cost + 1.0) * (cost + 1.0);
3203 }
3204 else {
3205 cost = 0.0;
3206 }
3207 }
3208 chi2 = track2.GetPredictedChi2(c,h01) + cost;
3209 clfound++;
3210
3211 if (chi2 > maxChi2[1]) {
3212 continue;
3213 }
3214 checkPoint[3]++;
3215
3216 // Store the clusters in the road
3217 detector = c->GetDetector();
3218 for (Int_t ih = 2; ih < 9; ih++) {
3219 if (cl[ih][it] == 0) {
3220 cl[ih][it] = c;
3221 indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
3222 break;
3223 }
3224 }
3225
3226 if (chi2 < maxChi2[0]) {
3227 maxChi2[1] = maxChi2[0];
3228 maxChi2[0] = chi2;
3229 indexes[1][it] = indexes[0][it];
3230 cl[1][it] = cl[0][it];
3231 indexes[0][it] = timeBin.GetIndex(i);
3232 cl[0][it] = c;
3233 continue;
3234 }
3235 maxChi2[1] = chi2;
3236 cl[1][it] = c;
3237 indexes[1][it] = timeBin.GetIndex(i);
3238
3239 }
3240
3241 for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++) {
3242 fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]);
3243 }
3244
3245 if (checkPoint[3]) {
3246 if (track->GetSignedPt() > 0) {
3247 fHMinYPos->Fill(minY);
3248 }
3249 else {
3250 fHMinYNeg->Fill(minY);
3251 }
3252 fHMinD->Fill(minD[0],minD[1]);
3253 }
3254
3255 if (cl[0][it]) {
3256 nfound++;
3257 xmean += x[it];
3258 }
3259
3260 }
3261
3262 if (nfound < 4) {
3263 return 0;
3264 }
3265 xmean /= Float_t(nfound); // Middle x
3266 track2.PropagateTo(xmean); // Propagate track to the center
3267
3268 //
3269 // Choose one of the variants
3270 //
3271 Float_t sumz = 0.0;
3272 Float_t sum = 0.0;
3273 Double_t sumdy = 0.0;
3274 Double_t sumdy2 = 0.0;
3275 Double_t sumx = 0.0;
3276 Double_t sumxy = 0.0;
3277 Double_t sumx2 = 0.0;
3278 Double_t mpads = 0.0;
3279
3280 Int_t changes[kN2];
3281
3282 Int_t ngood[kN2];
3283 Int_t nbad[kN2];
3284
3285 Double_t meanz[kN2];
3286 Double_t moffset[kN2]; // Mean offset
3287 Double_t mean[kN2]; // Mean value
3288 Double_t angle[kN2]; // Angle
3289
3290 Double_t smoffset[kN2]; // Sigma of mean offset
3291 Double_t smean[kN2]; // Sigma of mean value
3292 Double_t sangle[kN2]; // Sigma of angle
3293 Double_t smeanangle[kN2]; // Correlation
3294
3295 Double_t sigmas[kN2];
3296 Double_t tchi2s[kN2]; // Chi2s for tracklet
3297
3298 for (Int_t it = 0; it < kN2; it++) {
3299
3300 ngood[it] = 0;
3301 nbad[it] = 0;
3302
3303 meanz[it] = 0.0;
3304 moffset[it] = 0.0; // Mean offset
3305 mean[it] = 0.0; // Mean value
3306 angle[it] = 0.0; // Angle
3307
3308 smoffset[it] = 1.0e5; // Sigma of mean offset
3309 smean[it] = 1.0e5; // Sigma of mean value
3310 sangle[it] = 1.0e5; // Sigma of angle
3311 smeanangle[it] = 0.0; // Correlation
3312
3313 sigmas[it] = 1.0e5;
3314 tchi2s[it] = 1.0e5; // Chi2s for tracklet
3315
3316 }
3317
3318 //
3319 // Calculate zmean
3320 //
3321 for (Int_t it = 0; it < t1 - t0; it++) {
3322 if (!cl[0][it]) {
3323 continue;
3324 }
3325 for (Int_t dt = -3; dt <= 3; dt++) {
3326 if (it+dt < 0) {
3327 continue;
3328 }
3329 if (it+dt > t1-t0) {
3330 continue;
3331 }
3332 if (!cl[0][it+dt]) {
3333 continue;
3334 }
3335 zmean[it] += cl[0][it+dt]->GetZ();
3336 nmean[it] += 1.0;
3337 }
3338 zmean[it] /= nmean[it];
3339 }
3340
3341 for (Int_t it = 0; it < t1 - t0; it++) {
3342
3343 best[0][it] = 0;
3344
3345 for (Int_t ih = 0; ih < 10; ih++) {
3346 dz[ih][it] = -100.0;
3347 dy[ih][it] = -100.0;
3348 if (!cl[ih][it]) {
3349 continue;
3350 }
3351 Double_t xcluster = cl[ih][it]->GetX();
3352 Double_t ytrack;
3353 Double_t ztrack;
3354 track2.GetProlongation(xcluster,ytrack,ztrack );
3355 dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // Calculate distance from track in z
3356 dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01 - ytrack; // and in y
3357 }
3358
3359 // Minimize changes
3360 if (!cl[0][it]) {
3361 continue;
3362 }
3363 if ((TMath::Abs(cl[0][it]->GetZ()-zmean[it]) > padlength * 0.8) &&
3364 (cl[1][it])) {
3365 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it]) < padlength * 0.5) {
3366 best[0][it] = 1;
3367 }
3368 }
3369
3370 }
3371
3372 //
3373 // Iterative choice of "best path"
3374 //
3375 Int_t label = TMath::Abs(track->GetLabel());
3376 Int_t bestiter = 0;
3377
3378 for (Int_t iter = 0; iter < 9; iter++) {
3379
3380 changes[iter] = 0;
3381 sumz = 0;
3382 sum = 0;
3383 sumdy = 0;
3384 sumdy2 = 0;
3385 sumx = 0;
3386 sumx2 = 0;
3387 sumxy = 0;
3388 mpads = 0;
3389 ngood[iter] = 0;
3390 nbad[iter] = 0;
3391
3392 // Linear fit
3393 for (Int_t it = 0; it < t1 - t0; it++) {
3394
3395 if (!cl[best[iter][it]][it]) {
3396 continue;
3397 }
3398
3399 // Calculates pad-row changes
3400 Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3401 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3402 for (Int_t itd = it - 1; itd >= 0; itd--) {
3403 if (cl[best[iter][itd]][itd]) {
3404 zbefore = cl[best[iter][itd]][itd]->GetZ();
3405 break;
3406 }
3407 }
3408 for (Int_t itd = it + 1; itd < t1 - t0; itd++) {
3409 if (cl[best[iter][itd]][itd]) {
3410 zafter = cl[best[iter][itd]][itd]->GetZ();
3411 break;
3412 }
3413 }
3414 if ((TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore) > 0.1) &&
3415 (TMath::Abs(cl[best[iter][it]][it]->GetZ()- zafter) > 0.1)) {
3416 changes[iter]++;
3417 }
3418
3419 // Distance to reference x
3420 Double_t dx = x[it]-xmean;
3421 sumz += cl[best[iter][it]][it]->GetZ();
3422 sum++;
3423 sumdy += dy[best[iter][it]][it];
3424 sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3425 sumx += dx;
3426 sumx2 += dx*dx;
3427 sumxy += dx*dy[best[iter][it]][it];
3428 mpads += cl[best[iter][it]][it]->GetNPads();
3429 if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3430 (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3431 (cl[best[iter][it]][it]->GetLabel(2) == label)) {
3432 ngood[iter]++;
3433 }
3434 else {
3435 nbad[iter]++;
3436 }
3437
3438 }
3439
3440 //
3441 // Calculates line parameters
3442 //
3443 Double_t det = sum*sumx2 - sumx*sumx;
3444 angle[iter] = (sum*sumxy - sumx*sumdy) / det;
3445 mean[iter] = (sumx2*sumdy - sumx*sumxy) / det;
3446 meanz[iter] = sumz / sum;
3447 moffset[iter] = sumdy / sum;
3448 mpads /= sum; // Mean number of pads
3449
3450 Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3451 Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3452
3453 for (Int_t it = 0; it < t1 - t0; it++) {
3454 if (!cl[best[iter][it]][it]) {
3455 continue;
3456 }
3457 Double_t dx = x[it] - xmean;
3458 Double_t ytr = mean[iter] + angle[iter] * dx;
3459 sigma2 += (dy[best[iter][it]][it] - ytr)
3460 * (dy[best[iter][it]][it] - ytr);
3461 sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3462 * (dy[best[iter][it]][it] - moffset[iter]);
3463 sum++;
3464 }
3465 sigma2 /= (sum - 2); // Normalized residuals
3466 sigma1 /= (sum - 1); // Normalized residuals
3467 smean[iter] = sigma2 * (sumx2 / det); // Estimated error2 of mean
3468 sangle[iter] = sigma2 * ( sum / det); // Estimated error2 of angle
3469 smeanangle[iter] = sigma2 * (-sumx / det); // Correlation
3470 sigmas[iter] = TMath::Sqrt(sigma1);
3471 smoffset[iter] = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3472
3473 //
3474 // Iterative choice of "better path"
3475 //
3476 for (Int_t it = 0; it < t1 - t0; it++) {
3477
3478 if (!cl[best[iter][it]][it]) {
3479 continue;
3480 }
3481
3482 // Add unisochronity + angular effect contribution
3483 Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3484 Double_t sweight = 1.0/sigmatr2 + 1.0/track->GetSigmaY2();
3485 Double_t weighty = (moffset[iter] / sigmatr2) / sweight; // Weighted mean
3486 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3487 Double_t mindist = 100000.0;
3488 Int_t ihbest = 0;
3489
3490 for (Int_t ih = 0; ih < kN2; ih++) {
3491 if (!cl[ih][it]) {
3492 break;
3493 }
3494 Double_t dist2 = (dy[ih][it] - weighty) / sigmacl;
3495 dist2 *= dist2; // Chi2 distance
3496 if (dist2 < mindist) {
3497 mindist = dist2;
3498 ihbest = ih;
3499 }
3500 }
3501
3502 best[iter+1][it] = ihbest;
3503
3504 }
3505
3506 //
3507 // Update best hypothesy if better chi2 according tracklet position and angle
3508 //
3509 sy2 = smean[iter] + track->GetSigmaY2();
3510 Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee;
3511 Double_t say = track->GetSigmaSnpY(); // track->fCey;
3512
3513 Double_t detchi = sy2*sa2 - say*say;
3514 Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix
3515
3516 Double_t chi20 = mean[bestiter] * mean[bestiter] * invers[0]
3517 + angle[bestiter] * angle[bestiter] * invers[1]
3518 + 2.0 * mean[bestiter] * angle[bestiter] * invers[2];
3519 Double_t chi21 = mean[iter] * mean[iter] * invers[0]
3520 + angle[iter] * angle[iter] * invers[1]
3521 + 2.0 * mean[iter] * angle[iter] * invers[2];
3522 tchi2s[iter] = chi21;
3523
3524 if ((changes[iter] <= changes[bestiter]) &&
3525 (chi21 < chi20)) {
3526 bestiter = iter;
3527 }
3528
3529 }
3530
3531 //
3532 // Set clusters
3533 //
3534 Double_t sigma2 = sigmas[0]; // Choose as sigma from 0 iteration
3535 Short_t maxpos = -1;
3536 Float_t maxcharge = 0.0;
3537 Short_t maxpos4 = -1;
3538 Float_t maxcharge4 = 0.0;
3539 Short_t maxpos5 = -1;
3540 Float_t maxcharge5 = 0.0;
3541
3542 Double_t exB = AliTRDCommonParam::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0));
3543 Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
3544 if (mpads > 3.5) {
3545 expectederr += (mpads - 3.5) * 0.04;
3546 }
3547 if (changes[bestiter] > 1) {
3548 expectederr += changes[bestiter] * 0.01;
3549 }
3550 expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
3551
3552 for (Int_t it = 0; it < t1 - t0; it++) {
3553
3554 if (!cl[best[bestiter][it]][it]) {
3555 continue;
3556 }
3557
3558 // Set cluster error
3559 ((AliCluster*)cl[best[bestiter][it]][it])->SetSigmaY2(expectederr);
3560 if (!cl[best[bestiter][it]][it]->IsUsed()) {
3561 cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
3562 }
3563
3564 // Time bins with maximal charge
3565 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
3566 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3567 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3568 }
3569
3570 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3571 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
3572 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3573 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3574 }
3575 }
3576
3577 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3578 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
3579 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3580 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3581 }
3582 }
3583
3584 // Time bins with maximal charge
3585 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
3586 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3587 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3588 }
3589
3590 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3591 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
3592 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3593 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3594 }
3595 }
3596
3597 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3598 if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
3599 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3600 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3601 }
3602 }
3603
3604 clusters[it+t0] = indexes[best[bestiter][it]][it];
3605
3606 }
3607
3608 // Set tracklet parameters
3609 Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
3610 if (mpads > 3.5) {
3611 trackleterr2 += (mpads - 3.5) * 0.04;
3612 }
3613 trackleterr2 += changes[bestiter] * 0.01;
3614 trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
3615 trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
3616
3617 tracklet.Set(xmean
3618 ,track2.GetY() + moffset[bestiter]
3619 ,meanz[bestiter]
3620 ,track2.GetAlpha()
3621 ,trackleterr2);
3622 tracklet.SetTilt(h01);
3623 tracklet.SetP0(mean[bestiter]);
3624 tracklet.SetP1(angle[bestiter]);
3625 tracklet.SetN(nfound);
3626 tracklet.SetNCross(changes[bestiter]);
3627 tracklet.SetPlane(layer);
3628 tracklet.SetSigma2(expectederr);
3629 tracklet.SetChi2(tchi2s[bestiter]);
3630 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
3631 track->SetTracklets(layer,tracklet);
3632 track->SetNWrong(track->GetNWrong() + nbad[0]);
3633
3634 //
3635 // Debuging part
3636 //
3637 TClonesArray array0("AliTRDcluster");
3638 TClonesArray array1("AliTRDcluster");
3639 array0.ExpandCreateFast(t1 - t0 + 1);
3640 array1.ExpandCreateFast(t1 - t0 + 1);
3641 TTreeSRedirector &cstream = *fDebugStreamer;
3642 AliTRDcluster dummy;
3643 Double_t dy0[100];
3644 Double_t dyb[100];
3645
3646 for (Int_t it = 0; it < t1 - t0; it++) {
3647 dy0[it] = dy[0][it];
3648 dyb[it] = dy[best[bestiter][it]][it];
3649 if (cl[0][it]) {
3650 new(array0[it]) AliTRDcluster(*cl[0][it]);
3651 }
3652 else {
3653 new(array0[it]) AliTRDcluster(dummy);
3654 }
3655 if(cl[best[bestiter][it]][it]) {
3656 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3657 }
3658 else{
3659 new(array1[it]) AliTRDcluster(dummy);
3660 }
3661 }
3662
3663 TGraph graph0(t1-t0,x,dy0);
3664 TGraph graph1(t1-t0,x,dyb);
3665 TGraph graphy(t1-t0,x,yt);
3666 TGraph graphz(t1-t0,x,zt);
3667
3668 if (fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 0) {
3669 cstream << "tracklet"
3670 << "track.=" << track // Track parameters
3671 << "tany=" << tany // Tangent of the local track angle
3672 << "xmean=" << xmean // Xmean - reference x of tracklet
3673 << "tilt=" << h01 // Tilt angle
3674 << "nall=" << nall // Number of foundable clusters
3675 << "nfound=" << nfound // Number of found clusters
3676 << "clfound=" << clfound // Total number of found clusters in road
3677 << "mpads=" << mpads // Mean number of pads per cluster
3678 << "layer=" << layer // Layer number
3679 << "detector=" << detector // Detector number
3680 << "road=" << road // The width of the used road
3681 << "graph0.=" << &graph0 // x - y = dy for closest cluster
3682 << "graph1.=" << &graph1 // x - y = dy for second closest cluster
3683 << "graphy.=" << &graphy // y position of the track
3684 << "graphz.=" << &graphz // z position of the track
3685 //<< "fCl.=" << &array0 // closest cluster
3686 //<< "fCl2.=" << &array1 // second closest cluster
3687 << "maxpos=" << maxpos // Maximal charge postion
3688 << "maxcharge=" << maxcharge // Maximal charge
3689 << "maxpos4=" << maxpos4 // Maximal charge postion - after bin 4
3690 << "maxcharge4=" << maxcharge4 // Maximal charge - after bin 4
3691 << "maxpos5=" << maxpos5 // Maximal charge postion - after bin 5
3692 << "maxcharge5=" << maxcharge5 // Maximal charge - after bin 5
3693 << "bestiter=" << bestiter // Best iteration number
3694 << "tracklet.=" << &tracklet // Corrspond to the best iteration
3695 << "tchi20=" << tchi2s[0] // Chi2 of cluster in the 0 iteration
3696 << "tchi2b=" << tchi2s[bestiter] // Chi2 of cluster in the best iteration
3697 << "sigmas0=" << sigmas[0] // Residuals sigma
3698 << "sigmasb=" << sigmas[bestiter] // Residulas sigma
3699 << "ngood0=" << ngood[0] // Number of good clusters in 0 iteration
3700 << "nbad0=" << nbad[0] // Number of bad clusters in 0 iteration
3701 << "ngoodb=" << ngood[bestiter] // in best iteration
3702 << "nbadb=" << nbad[bestiter] // in best iteration
3703 << "changes0=" << changes[0] // Changes of pardrows in iteration number 0
3704 << "changesb=" << changes[bestiter] // Changes of pardrows in best iteration
3705 << "moffset0=" << moffset[0] // Offset fixing angle in iter=0
3706 << "smoffset0=" << smoffset[0] // Sigma of offset fixing angle in iter=0
3707 << "moffsetb=" << moffset[bestiter] // Offset fixing angle in iter=best
3708 << "smoffsetb=" << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
3709 << "mean0=" << mean[0] // Mean dy in iter=0;
3710 << "smean0=" << smean[0] // Sigma of mean dy in iter=0
3711 << "meanb=" << mean[bestiter] // Mean dy in iter=best
3712 << "smeanb=" << smean[bestiter] // Sigma of mean dy in iter=best
3713 << "angle0=" << angle[0] // Angle deviation in the iteration number 0
3714 << "sangle0=" << sangle[0] // Sigma of angular deviation in iteration number 0
3715 << "angleb=" << angle[bestiter] // Angle deviation in the best iteration
3716 << "sangleb=" << sangle[bestiter] // Sigma of angle deviation in the best iteration
3717 << "expectederr=" << expectederr // Expected error of cluster position
3718 << "\n";
3719 }
3720
3721 return nfound;
3722
3723}
3724
3725//_____________________________________________________________________________
3726Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
3727 , Int_t *outlist, Bool_t down)
3728{
3729 //
3730 // Sort eleements according occurancy
3731 // The size of output array has is 2*n
3732 //
3733
3734 if (n <= 0) {
3735 return 0;
3736 }
3737
3738 Int_t *sindexS = new Int_t[n]; // Temporary array for sorting
3739 Int_t *sindexF = new Int_t[2*n];
3740 for (Int_t i = 0; i < n; i++) {
3741 sindexF[i] = 0;
3742 }
3743
3744 TMath::Sort(n,inlist,sindexS,down);
3745
3746 Int_t last = inlist[sindexS[0]];
3747 Int_t val = last;
3748 sindexF[0] = 1;
3749 sindexF[0+n] = last;
3750 Int_t countPos = 0;
3751
3752 // Find frequency
3753 for (Int_t i = 1; i < n; i++) {
3754 val = inlist[sindexS[i]];
3755 if (last == val) {
3756 sindexF[countPos]++;
3757 }
3758 else {
3759 countPos++;
3760 sindexF[countPos+n] = val;
3761 sindexF[countPos]++;
3762 last = val;
3763 }
3764 }
3765 if (last == val) {
3766 countPos++;
3767 }
3768
3769 // Sort according frequency
3770 TMath::Sort(countPos,sindexF,sindexS,kTRUE);
3771
3772 for (Int_t i = 0; i < countPos; i++) {
3773 outlist[2*i ] = sindexF[sindexS[i]+n];
3774 outlist[2*i+1] = sindexF[sindexS[i]];
3775 }
3776
3777 delete [] sindexS;
3778 delete [] sindexF;
3779
3780 return countPos;
3781
3782}
3783
3784//_____________________________________________________________________________
3785AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed * const seeds, Double_t *params)
3786{
3787 //
3788 // Build a TRD track out of tracklet candidates
3789 //
3790 // Parameters :
3791 // seeds : array of tracklets
3792 // params : track parameters (see MakeSeeds() function body for a detailed description)
3793 //
3794 // Output :
3795 // The TRD track.
3796 //
3797 // Detailed description
3798 //
3799 // To be discussed with Marian !!
3800 //
3801
3802 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3803 Int_t nTimeBins = cal->GetNumberOfTimeBinsDCS();
3804
3805
3806 Double_t alpha = AliTRDgeometry::GetAlpha();
3807 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
3808 Double_t c[15];
3809
3810 c[ 0] = 0.2;
3811 c[ 1] = 0.0; c[ 2] = 2.0;
3812 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
3813 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
3814 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
3815
3816 Int_t index = 0;
3817 AliTRDcluster *cl = 0;
3818
3819 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
3820
3821 if (seeds[ilayer].IsOK()) {
3822 for (Int_t itime = nTimeBins-1; itime > 0; itime--) {
3823 if (seeds[ilayer].GetIndexes(itime) > 0) {
3824 index = seeds[ilayer].GetIndexes(itime);
3825 cl = seeds[ilayer].GetClusters(itime);
3826 //printf("l[%d] index[%d] tb[%d] cptr[%p]\n", ilayer, index, itime, cl);
3827 break;
3828 }
3829 }
3830 }
3831 if (index > 0) {
3832 break;
3833 }
3834 }
3835 if (cl == 0) {
3836 return 0;
3837 }
3838
3839 AliTRDtrack *track = new AliTRDtrack(cl
3840 ,index
3841 ,&params[1]
3842 ,c
3843 ,params[0]
3844 ,params[6]*alpha+shift);
3845 // SetCluster(cl, 0); // A. Bercuci
3846 track->PropagateTo(params[0]-5.0);
3847 track->ResetCovariance(1);
3848
3849 Int_t rc = FollowBackProlongation(*track);
3850 if (rc < 30) {
3851 delete track;
3852 track = 0;
3853 }
3854 else {
3855 track->CookdEdx();
3856 track->CookdEdxTimBin(-1);
3857 CookLabel(track,0.9);
3858 }
3859
3860 return track;
3861
3862}
3863
3864//_____________________________________________________________________________
3865void AliTRDtracker::InitLogHists()
3866{
3867 //
3868 // Create the log histograms
3869 //
3870
3871 fHBackfit = new TH1D("logTRD_backfit" ,""
3872 , 40,-0.5, 39.5);
3873 fHRefit = new TH1D("logTRD_refit" ,""
3874 , 40,-0.5, 39.5);
3875 fHClSearch = new TH1D("logTRD_clSearch",""
3876 , 60,-0.5, 59.5);
3877
3878 fHX = new TH1D("logTRD_X" ,";x (cm)"
3879 , 200, 50, 400);
3880 fHNCl = new TH1D("logTRD_ncl" ,""
3881 , 40,-0.5, 39.5);
3882 fHNClTrack = new TH1D("logTRD_nclTrack",""
3883 , 180,-0.5,179.5);
3884
3885 fHMinYPos = new TH1D("logTRD_minYPos" ,";#delta Y (cm)"
3886 , 400, -6, 6);
3887 fHMinYNeg = new TH1D("logTRD_minYNeg" ,";#delta Y (cm)"
3888 , 400, -6, 6);
3889 fHMinZ = new TH1D("logTRD_minZ" ,";#delta Z (cm)"
3890 , 400, -20, 20);
3891 fHMinD = new TH2D("logTRD_minD" ,";#delta Y (cm);#delta Z (cm)"
3892 , 100, -6, 6
3893 , 100, -50, 50);
3894
3895 fHDeltaX = new TH1D("logTRD_deltaX" ,";#delta X (cm)"
3896 , 100, -5, 5);
3897 fHXCl = new TH1D("logTRD_xCl" ,";cluster x position (cm)"
3898 ,1000, 280, 380);
3899
3900 const Char_t *nameFindCl[4] = { "logTRD_clY"
3901 , "logTRD_clZ"
3902 , "logTRD_clB"
3903 , "logTRD_clG" };
3904
3905 for (Int_t i = 0; i < 4; i++) {
3906 fHFindCl[i] = new TH1D(nameFindCl[i],"",30,-0.5,29.5);
3907 }
3908
3909}
3910
3911//_____________________________________________________________________________
3912void AliTRDtracker::SaveLogHists()
3913{
3914 //
3915 // Save the log histograms in AliESDs.root
3916 //
3917
3918 TDirectory *sav = gDirectory;
3919 TFile *logFile = 0;
3920
3921 TSeqCollection *col = gROOT->GetListOfFiles();
3922 Int_t nn = col->GetEntries();
3923 for (Int_t i = 0; i < nn; i++) {
3924 logFile = (TFile *) col->At(i);
3925 if (strstr(logFile->GetName(),"AliESDs.root")) {
3926 break;
3927 }
3928 }
3929
3930 logFile->cd();
3931
3932 fHBackfit->Write(fHBackfit->GetName(),TObject::kOverwrite);
3933 fHRefit->Write(fHRefit->GetName(),TObject::kOverwrite);
3934 fHClSearch->Write(fHClSearch->GetName(),TObject::kOverwrite);
3935 fHX->Write(fHX->GetName(),TObject::kOverwrite);
3936 fHNCl->Write(fHNCl->GetName(),TObject::kOverwrite);
3937 fHNClTrack->Write(fHNClTrack->GetName(),TObject::kOverwrite);
3938
3939 fHMinYPos->Write(fHMinYPos->GetName(),TObject::kOverwrite);
3940 fHMinYNeg->Write(fHMinYNeg->GetName(),TObject::kOverwrite);
3941 fHMinD->Write(fHMinD->GetName(),TObject::kOverwrite);
3942 fHMinZ->Write(fHMinZ->GetName(),TObject::kOverwrite);
3943
3944 fHDeltaX->Write(fHDeltaX->GetName(),TObject::kOverwrite);
3945 fHXCl->Write(fHXCl->GetName(),TObject::kOverwrite);
3946
3947 for (Int_t i = 0; i < 4; i++) {
3948 fHFindCl[i]->Write(fHFindCl[i]->GetName(),TObject::kOverwrite);
3949 }
3950
3951 logFile->Flush();
3952
3953 sav->cd();
3954
3955}