]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtrackerV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug in ZP positioning corrected
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerV1.cxx
CommitLineData
e4f2f73d 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// Track finder //
21// //
22// Authors: //
23// Alex Bercuci <A.Bercuci@gsi.de> //
24// Markus Fasel <M.Fasel@gsi.de> //
25// //
26///////////////////////////////////////////////////////////////////////////////
27
28#include <Riostream.h>
29#include <stdio.h>
30#include <string.h>
31
32#include <TBranch.h>
33#include <TFile.h>
34#include <TGraph.h>
35#include <TH1D.h>
36#include <TH2D.h>
37#include <TLinearFitter.h>
38#include <TObjArray.h>
39#include <TROOT.h>
40#include <TTree.h>
41#include <TClonesArray.h>
42#include <TRandom.h>
43#include <TTreeStream.h>
44
45#include "AliLog.h"
46#include "AliESDEvent.h"
47#include "AliAlignObj.h"
48#include "AliRieman.h"
49#include "AliTrackPointArray.h"
50
51#include "AliTRDtracker.h"
52#include "AliTRDtrackerV1.h"
53#include "AliTRDgeometry.h"
54#include "AliTRDpadPlane.h"
55#include "AliTRDgeometry.h"
56#include "AliTRDcluster.h"
57#include "AliTRDtrack.h"
58#include "AliTRDseed.h"
59#include "AliTRDcalibDB.h"
60#include "AliTRDCommonParam.h"
61#include "AliTRDReconstructor.h"
62#include "AliTRDCalibraFillHisto.h"
63#include "AliTRDtrackerFitter.h"
64#include "AliTRDstackLayer.h"
65#include "AliTRDrecoParam.h"
66#include "AliTRDseedV1.h"
0906e73e 67#include "AliTRDtrackV1.h"
68#include "Cal/AliTRDCalDet.h"
e4f2f73d 69
70#define DEBUG
71
72ClassImp(AliTRDtrackerV1)
d76231c8 73Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
e4f2f73d 74 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
75 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
76 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
77};
78
79//____________________________________________________________________
80AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p)
81 :AliTRDtracker()
82 ,fSieveSeeding(0)
0906e73e 83 ,fTracklets(0x0)
fbb2ea06 84 ,fRecoParam(p)
e4f2f73d 85 ,fFitter(0x0)
e4f2f73d 86{
87 //
88 // Default constructor. Nothing is initialized.
89 //
90
91}
92
93//____________________________________________________________________
94AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p)
95 :AliTRDtracker(in)
96 ,fSieveSeeding(0)
0906e73e 97 ,fTracklets(0x0)
e4f2f73d 98 ,fRecoParam(p)
99 ,fFitter(0x0)
e4f2f73d 100{
101 //
102 // Standard constructor.
103 // Setting of the geometry file, debug output (if enabled)
104 // and initilize fitter helper.
105 //
106
107 fFitter = new AliTRDtrackerFitter();
108
109#ifdef DEBUG
0906e73e 110 fFitter->SetDebugStream(fDebugStreamer);
e4f2f73d 111#endif
112
113}
114
115//____________________________________________________________________
116AliTRDtrackerV1::~AliTRDtrackerV1()
117{
118 //
119 // Destructor
120 //
121
e4f2f73d 122 if(fFitter) delete fFitter;
123 if(fRecoParam) delete fRecoParam;
0906e73e 124 if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
e4f2f73d 125}
126
127//____________________________________________________________________
128Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
129{
130 //
131 // Steering stand alone tracking for full TRD detector
132 //
133 // Parameters :
134 // esd : The ESD event. On output it contains
135 // the ESD tracks found in TRD.
136 //
137 // Output :
138 // Number of tracks found in the TRD detector.
139 //
140 // Detailed description
141 // 1. Launch individual SM trackers.
142 // See AliTRDtrackerV1::Clusters2TracksSM() for details.
143 //
144
145 if(!fRecoParam){
146 AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
147 return 0;
148 }
149
150 //AliInfo("Start Track Finder ...");
151 Int_t ntracks = 0;
152 for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
153 //AliInfo(Form("Processing supermodule %i ...", ism));
154 ntracks += Clusters2TracksSM(fTrSec[ism], esd);
155 }
156 AliInfo(Form("Found %d TRD tracks.", ntracks));
157 return ntracks;
158}
159
0906e73e 160
161//_____________________________________________________________________________
162Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t /*index*/, AliTrackPoint &/*p*/) const
163{
164 //AliInfo(Form("Asking for tracklet %d", index));
165
166 if(index<0) return kFALSE;
167 //AliTRDseedV1 *tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index);
168 // etc
169 return kTRUE;
170}
171
172
173//_____________________________________________________________________________
174Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
175{
176 //
177 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
178 // backpropagated by the TPC tracker. Each seed is first propagated
179 // to the TRD, and then its prolongation is searched in the TRD.
180 // If sufficiently long continuation of the track is found in the TRD
181 // the track is updated, otherwise it's stored as originaly defined
182 // by the TPC tracker.
183 //
184
185 Int_t found = 0; // number of tracks found
186 Float_t foundMin = 20.0;
187
188 AliTRDseed::SetNTimeBins(fTimeBinsPerPlane);
189 Int_t nSeed = event->GetNumberOfTracks();
190 if(!nSeed){
191 // run stand alone tracking
192 if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
193 return 0;
194 }
195
196 Float_t *quality = new Float_t[nSeed];
197 Int_t *index = new Int_t[nSeed];
198 for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
199 AliESDtrack *seed = event->GetTrack(iSeed);
200 Double_t covariance[15];
201 seed->GetExternalCovariance(covariance);
202 quality[iSeed] = covariance[0] + covariance[2];
203 }
204 // Sort tracks according to covariance of local Y and Z
205 TMath::Sort(nSeed,quality,index,kFALSE);
206
207 // Backpropagate all seeds
208 for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
209
210 // Get the seeds in sorted sequence
211 AliESDtrack *seed = event->GetTrack(index[iSeed]);
212
213 // Check the seed status
214 ULong_t status = seed->GetStatus();
215 if ((status & AliESDtrack::kTPCout) == 0) continue;
216 if ((status & AliESDtrack::kTRDout) != 0) continue;
217
218 // Do the back prolongation
219 Int_t lbl = seed->GetLabel();
220 AliTRDtrackV1 *track = new AliTRDtrackV1(*seed);
221 //track->Print();
222 track->SetSeedLabel(lbl);
223 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup
224 fNseeds++;
225 Float_t p4 = track->GetC();
226 Int_t expectedClr = FollowBackProlongation(*track);
227 //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2()));
228 //track->Print();
229
230 //Double_t cov[15];
231 //seed->GetExternalCovariance(cov);
232 //AliInfo(Form("track %d cov[%f %f] 0", index[iSeed], cov[0], cov[2]));
233
234 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
235 (track->Pt() > 0.8)) {
236 //
237 // Make backup for back propagation
238 //
239 Int_t foundClr = track->GetNumberOfClusters();
240 if (foundClr >= foundMin) {
241 //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
242 track->CookdEdx();
243 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
244 CookLabel(track,1 - fgkLabelFraction);
245 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
246
247
248 //seed->GetExternalCovariance(cov);
249 //AliInfo(Form("track %d cov[%f %f] 0 test", index[iSeed], cov[0], cov[2]));
250
251 // Sign only gold tracks
252 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
253 if ((seed->GetKinkIndex(0) == 0) &&
254 (track->Pt() < 1.5)) UseClusters(track);
255 }
256 Bool_t isGold = kFALSE;
257
258 // Full gold track
259 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
260 if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
261
262 isGold = kTRUE;
263 }
264 //seed->GetExternalCovariance(cov);
265 //AliInfo(Form("track %d cov[%f %f] 00", index[iSeed], cov[0], cov[2]));
266
267 // Almost gold track
268 if ((!isGold) && (track->GetNCross() == 0) &&
269 (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
270 //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
271 if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
272
273 isGold = kTRUE;
274 }
275 //seed->GetExternalCovariance(cov);
276 //AliInfo(Form("track %d cov[%f %f] 01", index[iSeed], cov[0], cov[2]));
277
278 if ((!isGold) && (track->GetBackupTrack())) {
279 if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
280 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
281 isGold = kTRUE;
282 }
283 }
284 //seed->GetExternalCovariance(cov);
285 //AliInfo(Form("track %d cov[%f %f] 02", index[iSeed], cov[0], cov[2]));
286
287 //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
288 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
289 //}
290 }
291 }
292 /**/
293
294 /**/
295 // Debug part of tracking
296/* TTreeSRedirector &cstream = *fDebugStreamer;
297 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.
298 if (AliTRDReconstructor::StreamLevel() > 0) {
299 if (track->GetBackupTrack()) {
300 cstream << "Tracks"
301 << "EventNrInFile=" << eventNrInFile
302 << "ESD.=" << seed
303 << "trd.=" << track
304 << "trdback.=" << track->GetBackupTrack()
305 << "\n";
306 }
307 else {
308 cstream << "Tracks"
309 << "EventNrInFile=" << eventNrInFile
310 << "ESD.=" << seed
311 << "trd.=" << track
312 << "trdback.=" << track
313 << "\n";
314 }
315 }*/
316 /**/
317
318 //seed->GetExternalCovariance(cov);
319 //AliInfo(Form("track %d cov[%f %f] 1", index[iSeed], cov[0], cov[2]));
320
321 // Propagation to the TOF (I.Belikov)
322 if (track->GetStop() == kFALSE) {
323 //AliInfo("Track not stopped in TRD ...");
324 Double_t xtof = 371.0;
325 Double_t xTOF0 = 370.0;
326
327 Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
328 if (TMath::Abs(c2) >= 0.99) {
329 delete track;
330 continue;
331 }
332
333 PropagateToX(*track,xTOF0,fgkMaxStep);
334
335 // Energy losses taken to the account - check one more time
336 c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
337 if (TMath::Abs(c2) >= 0.99) {
338 delete track;
339 continue;
340 }
341
342 //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
343 // fHBackfit->Fill(7);
344 //delete track;
345 // continue;
346 //}
347
348 Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
349 Double_t y;
350 track->GetYAt(xtof,GetBz(),y);
351 if (y > ymax) {
352 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
353 delete track;
354 continue;
355 }
356 }else if (y < -ymax) {
357 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
358 delete track;
359 continue;
360 }
361 }
362
363 if (track->PropagateTo(xtof)) {
364 //AliInfo("set kTRDout");
365 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
366
367 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
368 for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
369 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
370 }
371 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
372 }
373 //seed->SetTRDtrack(new AliTRDtrack(*track));
374 if (track->GetNumberOfClusters() > foundMin) found++;
375 }
376 } else {
377 //AliInfo("Track stopped in TRD ...");
378
379 if ((track->GetNumberOfClusters() > 15) &&
380 (track->GetNumberOfClusters() > 0.5*expectedClr)) {
381 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
382
383 //seed->SetStatus(AliESDtrack::kTRDStop);
384 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
385 for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
386 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
387 }
388 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
389 }
390 //seed->SetTRDtrack(new AliTRDtrack(*track));
391 found++;
392 }
393 }
394
395 //if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 )
396
397 seed->SetTRDQuality(track->StatusForTOF());
398 seed->SetTRDBudget(track->GetBudget(0));
399
400// if ((seed->GetStatus()&AliESDtrack::kTRDin)!=0 ) printf("TRDin ");
401// if ((seed->GetStatus()&AliESDtrack::kTRDbackup)!=0 ) printf("TRDbackup ");
402// if ((seed->GetStatus()&AliESDtrack::kTRDStop)!=0 ) printf("TRDstop ");
403// if ((seed->GetStatus()&AliESDtrack::kTRDout)!=0 ) printf("TRDout ");
404// printf("\n");
405 delete track;
406
407 //seed->GetExternalCovariance(cov);
408 //AliInfo(Form("track %d cov[%f %f] 2", index[iSeed], cov[0], cov[2]));
409 }
410
411
412 AliInfo(Form("Number of seeds: %d",fNseeds));
413 AliInfo(Form("Number of back propagated TRD tracks: %d",found));
414
415 //fSeeds->Clear();
416 fNseeds = 0;
417
418 delete [] index;
419 delete [] quality;
420
421 return 0;
422}
423
424
425//____________________________________________________________________
426Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event)
427{
428 //
429 // Refits tracks within the TRD. The ESD event is expected to contain seeds
430 // at the outer part of the TRD.
431 // The tracks are propagated to the innermost time bin
432 // of the TRD and the ESD event is updated
433 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
434 //
435
436 Int_t nseed = 0; // contor for loaded seeds
437 Int_t found = 0; // contor for updated TRD tracks
438 AliTRDtrackV1 track;
439 for (Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++) {
440 AliESDtrack *seed = event->GetTrack(itrack);
441 new(&track) AliTRDtrackV1(*seed);
442
443 if (track.GetX() < 270.0) {
444 seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
445 //AliInfo(Form("Remove for X = %7.3f [270.]\n", track.GetX()));
446 continue;
447 }
448
449 ULong_t status = seed->GetStatus();
450 if((status & AliESDtrack::kTRDout) == 0) continue;
451 if((status & AliESDtrack::kTRDin) != 0) continue;
452 nseed++;
453
454 track.ResetCovariance(50.0);
455
456 // do the propagation and processing
457 FollowProlongation(track);
458 // computes PID for track
459 track.CookPID();
460 // update calibration references using this track
461 //track.Calibrate();
462
463 // Prolongate to TPC
464 Double_t xTPC = 250.0;
465 if (PropagateToX(track, xTPC, fgkMaxStep)) { // -with update
466 seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit);
467 track.UpdateESDtrack(seed);
468 // Add TRD track to ESDfriendTrack
469 if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
470 AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
471 calibTrack->SetOwner();
472 seed->AddCalibObject(calibTrack);
473 }
474 found++;
475 } else { // - without update
476 AliTRDtrackV1 tt(*seed);
477 if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit);
478 }
479 }
480 AliInfo(Form("Number of loaded seeds: %d",nseed));
481 AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
482
483 return 0;
484}
485
486
487//____________________________________________________________________
488Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
489{
490// Extrapolates the TRD track in the TPC direction.
491//
492// Parameters
493// t : the TRD track which has to be extrapolated
494//
495// Output
496// number of clusters attached to the track
497//
498// Detailed description
499//
500// Starting from current radial position of track <t> this function
501// extrapolates the track through the 6 TRD layers. The following steps
502// are being performed for each plane:
503// 1. prepare track:
504// a. get plane limits in the local x direction
505// b. check crossing sectors
506// c. check track inclination
507// 2. search tracklet in the tracker list (see GetTracklet() for details)
508// 3. evaluate material budget using the geo manager
509// 4. propagate and update track using the tracklet information.
510//
511// Debug level 2
512//
513
514 //AliInfo("");
515 Int_t nClustersExpected = 0;
516 Int_t lastplane = 5; //GetLastPlane(&t);
517 for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
518 //AliInfo(Form("plane %d", iplane));
519 Int_t row1 = GetGlobalTimeBin(0, iplane, 0); // to be modified to the true time bin in the geometrical acceptance
520 //AliInfo(Form("row1 %d", row1));
521
522 // Propagate track close to the plane if neccessary
523 AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row1);
524 Double_t currentx = layer->GetX();
525 if (currentx < (-fgkMaxStep + t.GetX()))
526 if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
527
528 if (!AdjustSector(&t)) break;
529
530 Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
531 //AliInfo(Form("row0 %d", row0));
532
533 // Start global position
534 Double_t xyz0[3];
535 t.GetXYZ(xyz0);
536
537 // End global position
538 Double_t x = fTrSec[0]->GetLayer(row0)->GetX(), y, z;
539 if (!t.GetProlongation(x,y,z)) break;
540 Double_t xyz1[3];
541 xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
542 xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
543 xyz1[2] = z;
544
545 // Get material budget
546 Double_t param[7];
547 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
548 Double_t xrho= param[0]*param[4];
549 Double_t xx0 = param[1]; // Get mean propagation parameters
550
551 // Propagate and update
552 //Int_t sector = t.GetSector();
553 Int_t index = 0;
554 //AliInfo(Form("sector %d", sector));
555 AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
556 //AliInfo(Form("tracklet 0x%x @ %d", tracklet, index));
557 if(!tracklet) continue;
558
559 t.PropagateTo(tracklet->GetX0(), xx0, xrho); // not correct
560 if (!AdjustSector(&t)) break;
561
562 Double_t maxChi2 = t.GetPredictedChi2(tracklet);
563 if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){
564 nClustersExpected += tracklet->GetN();
565 }
566 }
567
568#ifdef DEBUG
569 if(AliTRDReconstructor::StreamLevel() > 1){
570 Int_t index;
571 for(int iplane=0; iplane<6; iplane++){
572 AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
573 if(!tracklet) continue;
574 t.SetTracklet(tracklet, iplane, index);
575 }
576
577 TTreeSRedirector &cstreamer = *fDebugStreamer;
578 cstreamer << "FollowProlongation"
579 << "ncl=" << nClustersExpected
580 << "track.=" << &t
581 << "\n";
582 }
583#endif
584
585 return nClustersExpected;
586
587}
588
589//_____________________________________________________________________________
590Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
591{
592// Extrapolates the TRD track in the TOF direction.
593//
594// Parameters
595// t : the TRD track which has to be extrapolated
596//
597// Output
598// number of clusters attached to the track
599//
600// Detailed description
601//
602// Starting from current radial position of track <t> this function
603// extrapolates the track through the 6 TRD layers. The following steps
604// are being performed for each plane:
605// 1. prepare track:
606// a. get plane limits in the local x direction
607// b. check crossing sectors
608// c. check track inclination
609// 2. build tracklet (see AliTRDseed::AttachClusters() for details)
610// 3. evaluate material budget using the geo manager
611// 4. propagate and update track using the tracklet information.
612//
613// Debug level 2
614//
615
616 Int_t nClustersExpected = 0;
617
618 // Loop through the TRD planes
619 for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
620 //AliInfo(Form("Processing plane %d ...", iplane));
621 // Get the global time bin for the first local time bin of the given plane
622 Int_t row0 = GetGlobalTimeBin(0, iplane, fTimeBinsPerPlane-1);
623
624 // Retrive first propagation layer in the chamber
625 AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row0);
626
627 // Get the X coordinates of the propagation layer for the first time bin
628 Double_t currentx = layer->GetX(); // what if X is not defined ???
629 if (currentx < t.GetX()) continue;
630
631 // Get the global time bin for the last local time bin of the given plane
632 Int_t row1 = GetGlobalTimeBin(0, iplane, 0);
633
634 // Propagate closer to the current chamber if neccessary
635 if (currentx > (fgkMaxStep + t.GetX()) && !PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
636
637 // Rotate track to adjacent sector if neccessary
638 if (!AdjustSector(&t)) break;
639 Int_t sector = Int_t(TMath::Abs(t.GetAlpha()/AliTRDgeometry::GetAlpha()));
640 if(t.GetAlpha() < 0) sector = AliTRDgeometry::Nsect() - sector-1;
641
642 //AliInfo(Form("sector %d [%f]", sector, t.GetAlpha()));
643
644 // Check whether azimuthal angle is getting too large
645 if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
646
647 //Calculate global entry and exit positions of the track in chamber (only track prolongation)
648 Double_t xyz0[3]; // entry point
649 t.GetXYZ(xyz0);
650 //printf("Entry global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz0[0], xyz0[1], xyz0[2]);
651
652 // Get local Y and Z at the X-position of the end of the chamber
653 Double_t x0 = fTrSec[sector]->GetLayer(row1)->GetX(), y, z;
654 if (!t.GetProlongation(x0, y, z)) break;
655 //printf("Exit local x[%7.3f] y[%7.3f] z[%7.3f]\n", x0, y, z);
656
657 Double_t xyz1[3]; // exit point
658 xyz1[0] = x0 * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
659 xyz1[1] = +x0 * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
660 xyz1[2] = z;
661
662 //printf("Exit global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz1[0], xyz1[1], xyz1[2]);
663 // Find tracklet along the path inside the chamber
664 AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
665 // if the track is not already build (e.g. stand alone tracker) we build it now.
666 if(!tracklet.GetN()){ // a better check has to be implemented TODO!!!!!!!
667
668 //AliInfo(Form("Building tracklet for plane %d ...", iplane));
669 // check if we are inside detection volume
670 Int_t ichmb = fGeom->GetChamber(xyz0[2], iplane);
671 if(ichmb<0) ichmb = fGeom->GetChamber(xyz1[2], iplane);
672 if(ichmb<0){
673 // here we should decide what to do with the track. The space between the pads in 2 chambers is 4cm+. Is it making sense to continue building the tracklet here TODO????
674 AliWarning(Form("Track prolongated in the interspace between TRD detectors in plane %d. Skip plane. To be fixed !", iplane));
675 continue;
676 }
677
678 // temporary until the functionalities of AliTRDpropagationLayer and AliTRDstackLayer are merged TODO
679 AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, ichmb);
680 Int_t nrows = pp->GetNrows();
681 Double_t stacklength = pp->GetRow0ROC() - pp->GetRowEndROC();/*(nrows - 2) * pp->GetLengthIPad() + 2 * pp->GetLengthOPad() + (nrows - 1) * pp->GetRowSpacing();*/
682 Double_t z0 = fGeom->GetRow0(iplane, ichmb, 0);
683
684 Int_t nClustersChmb = 0;
685 AliTRDstackLayer stackLayer[35];
686 for(int itb=0; itb<fTimeBinsPerPlane; itb++){
687 const AliTRDpropagationLayer ksmLayer(*(fTrSec[sector]->GetLayer(row1 - itb)));
688 stackLayer[itb] = ksmLayer;
689#ifdef DEBUG
690 stackLayer[itb].SetDebugStream(fDebugStreamer);
691#endif
692 stackLayer[itb].SetRange(z0 - stacklength, stacklength);
693 stackLayer[itb].SetSector(sector);
694 stackLayer[itb].SetStackNr(ichmb);
695 stackLayer[itb].SetNRows(nrows);
696 stackLayer[itb].SetRecoParam(fRecoParam);
697 stackLayer[itb].BuildIndices();
698 nClustersChmb += stackLayer[itb].GetNClusters();
699 }
700 //AliInfo(Form("Detector p[%d] c[%d]. Building tracklet from %d clusters ... ", iplane, ichmb, nClustersChmb));
701
702 tracklet.SetRecoParam(fRecoParam);
703 tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
704 tracklet.SetPadLength(pp->GetLengthIPad());
705 tracklet.SetPlane(iplane);
706 Int_t tbRange = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
707 //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
708 tracklet.SetNTimeBinsRange(tbRange);
709 tracklet.SetX0(x0);
710 tracklet.Init(&t);
711 if(!tracklet.AttachClustersIter(stackLayer, 1000.)) continue;
712
713 //if(!tracklet.AttachClusters(stackLayer, kTRUE)) continue;
714 //if(!tracklet.Fit()) continue;
715 }
716 Int_t ncl = tracklet.GetN();
717 //AliInfo(Form("N clusters %d", ncl));
718
719 // Discard tracklet if bad quality.
720 //Check if this condition is not already checked during building of the tracklet
721 if(ncl < fTimeBinsPerPlane * fRecoParam->GetFindableClusters()){
722 //AliInfo(Form("Discard tracklet for %d nclusters", ncl));
723 continue;
724 }
725
726 // load tracklet to the tracker and the track
727 Int_t index = SetTracklet(&tracklet);
728 t.SetTracklet(&tracklet, iplane, index);
729
730 // Calculate the mean material budget along the path inside the chamber
731 Double_t param[7];
732 AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
733 // The mean propagation parameters
734 Double_t xrho = param[0]*param[4]; // density*length
735 Double_t xx0 = param[1]; // radiation length
736
737 // Propagate and update track
738 t.PropagateTo(tracklet.GetX0(), xx0, xrho);
739 if (!AdjustSector(&t)) break;
740 Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
741 if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){
742 nClustersExpected += ncl;
743 }
744 // Reset material budget if 2 consecutive gold
745 if(iplane>0 && ncl + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
746
747 // Make backup of the track until is gold
748 // TO DO update quality check of the track.
749 // consider comparison with fTimeBinsRange
750 Float_t ratio0 = ncl / Float_t(fTimeBinsPerPlane);
751 //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
752 //printf("tracklet.GetChi2() %f [< 18.0]\n", tracklet.GetChi2());
753 //printf("ratio0 %f [> 0.8]\n", ratio0);
754 //printf("ratio1 %f [> 0.6]\n", ratio1);
755 //printf("ratio0+ratio1 %f [> 1.5]\n", ratio0+ratio1);
756 //printf("t.GetNCross() %d [== 0]\n", t.GetNCross());
757 //printf("TMath::Abs(t.GetSnp()) %f [< 0.85]\n", TMath::Abs(t.GetSnp()));
758 //printf("t.GetNumberOfClusters() %d [> 20]\n", t.GetNumberOfClusters());
759
760 if (//(tracklet.GetChi2() < 18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update
761 (ratio0 > 0.8) &&
762 //(ratio1 > 0.6) &&
763 //(ratio0+ratio1 > 1.5) &&
764 (t.GetNCross() == 0) &&
765 (TMath::Abs(t.GetSnp()) < 0.85) &&
766 (t.GetNumberOfClusters() > 20)) t.MakeBackupTrack();
767
768 } // end planes loop
769
770#ifdef DEBUG
771 if(AliTRDReconstructor::StreamLevel() > 1){
772 TTreeSRedirector &cstreamer = *fDebugStreamer;
773 cstreamer << "FollowBackProlongation"
774 << "ncl=" << nClustersExpected
775 << "track.=" << &t
776 << "\n";
777 }
778#endif
779
780 return nClustersExpected;
781}
782
783//____________________________________________________________________
784void AliTRDtrackerV1::UnloadClusters()
785{
786 //
787 // Clears the arrays of clusters and tracks. Resets sectors and timebins
788 //
789
790 Int_t i;
791 Int_t nentr;
792
793 nentr = fClusters->GetEntriesFast();
794 //AliInfo(Form("clearing %d clusters", nentr));
795 for (i = 0; i < nentr; i++) {
796 delete fClusters->RemoveAt(i);
797 }
798 fNclusters = 0;
799
800 nentr = fTracklets->GetEntriesFast();
801 //AliInfo(Form("clearing %d tracklets", nentr));
802 for (i = 0; i < nentr; i++) {
803 delete fTracklets->RemoveAt(i);
804 }
805
806 nentr = fSeeds->GetEntriesFast();
807 //AliInfo(Form("clearing %d seeds", nentr));
808 for (i = 0; i < nentr; i++) {
809 delete fSeeds->RemoveAt(i);
810 }
811
812 nentr = fTracks->GetEntriesFast();
813 //AliInfo(Form("clearing %d tracks", nentr));
814 for (i = 0; i < nentr; i++) {
815 delete fTracks->RemoveAt(i);
816 }
817
818 Int_t nsec = AliTRDgeometry::kNsect;
819 for (i = 0; i < nsec; i++) {
820 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
821 fTrSec[i]->GetLayer(pl)->Clear();
822 }
823 }
824
825}
826
827//____________________________________________________________________
828AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx)
829{
830// Find tracklet for TRD track <track>
831// Parameters
832// - track
833// - sector
834// - plane
835// - index
836// Output
837// tracklet
838// index
839// Detailed description
840//
841 idx = track->GetTrackletIndex(p);
842 //AliInfo(Form("looking for tracklet in plane %d idx %d [%d]", p, idx, track->GetTrackletIndex(p)));
843 AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
844 //AliInfo(Form("found 0x%x @ %d", tracklet, idx));
845
846// Int_t *index = track->GetTrackletIndexes();
847// for (UInt_t i = 0; i < 6; i++) AliInfo(Form("index[%d] = %d", i, index[i]));
848//
849// for (UInt_t i = 0; i < 6/*kMaxTimeBinIndex*/; i++) {
850// if (index[i] < 0) continue;
851//
852// tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index[i]);
853// if(!tracklet) break;
854//
855// if(tracklet->GetPlane() != p) continue;
856//
857// idx = index[i];
858// }
859
860 return tracklet;
861}
862
863//____________________________________________________________________
864Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 */*tracklet*/)
865{
866// Add this tracklet to the list of tracklets stored in the tracker
867//
868// Parameters
869// - tracklet : pointer to the tracklet to be added to the list
870//
871// Output
872// - the index of the new tracklet in the tracker tracklets list
873//
874// Detailed description
875// Build the tracklets list if it is not yet created (late initialization)
876// and adds the new tracklet to the list.
877//
878 if(!fTracklets){
879 fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
880 fTracklets->SetOwner(kTRUE);
881 }
882 Int_t nentries = fTracklets->GetEntriesFast();
883 //AliTRDseedV1 *t = new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
884 //AliInfo(Form("0x%x @ %d", t, nentries));
885 return nentries;
886}
887
e4f2f73d 888//____________________________________________________________________
889Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
890 , AliESDEvent *esd)
891{
892 //
893 // Steer tracking for one SM.
894 //
895 // Parameters :
896 // sector : Array of (SM) propagation layers containing clusters
897 // esd : The current ESD event. On output it contains the also
898 // the ESD (TRD) tracks found in this SM.
899 //
900 // Output :
901 // Number of tracks found in this TRD supermodule.
902 //
903 // Detailed description
904 //
905 // 1. Unpack AliTRDpropagationLayers objects for each stack.
906 // 2. Launch stack tracking.
907 // See AliTRDtrackerV1::Clusters2TracksStack() for details.
908 // 3. Pack results in the ESD event.
909 //
910
911 AliTRDpadPlane *pp = 0x0;
912
913 // allocate space for esd tracks in this SM
914 TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
915 esdTrackList.SetOwner();
916 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
917 Int_t nTimeBins = cal->GetNumberOfTimeBins();
918 const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
919
920 Int_t ntracks = 0;
921 Int_t nClStack = 0;
922 for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
923 AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
924
925 nClStack = 0;
926 //AliInfo(Form("Processing stack %i ...",istack));
927 //AliInfo("Building stack propagation layers ...");
928 for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
929 pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
930 Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad()
931 + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
932 //Debug
933 Double_t z0 = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
0906e73e 934 const AliTRDpropagationLayer ksmLayer(*(sector->GetLayer(ilayer)));
935 stackLayer[ilayer] = ksmLayer;
e4f2f73d 936#ifdef DEBUG
0906e73e 937 stackLayer[ilayer].SetDebugStream(fDebugStreamer);
e4f2f73d 938#endif
939 stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
940 stackLayer[ilayer].SetSector(sector->GetSector());
941 stackLayer[ilayer].SetStackNr(istack);
942 stackLayer[ilayer].SetNRows(pp->GetNrows());
943 stackLayer[ilayer].SetRecoParam(fRecoParam);
944 stackLayer[ilayer].BuildIndices();
945 nClStack += stackLayer[ilayer].GetNClusters();
946 }
947 //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
948 if(nClStack < kFindable) continue;
949 ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
950 }
951 //AliInfo(Form("Found %d tracks in SM", ntracks));
952
953 for(int itrack=0; itrack<ntracks; itrack++)
954 esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
955
956 return ntracks;
957}
958
959//____________________________________________________________________
960Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
961 , TClonesArray *esdTrackList)
962{
963 //
964 // Make tracks in one TRD stack.
965 //
966 // Parameters :
967 // layer : Array of stack propagation layers containing clusters
968 // esdTrackList : Array of ESD tracks found by the stand alone tracker.
969 // On exit the tracks found in this stack are appended.
970 //
971 // Output :
972 // Number of tracks found in this stack.
973 //
974 // Detailed description
975 //
976 // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
977 // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations.
978 // See AliTRDtrackerV1::MakeSeeds() for more details.
979 // 3. Arrange track candidates in decreasing order of their quality
980 // 4. Classify tracks in 5 categories according to:
981 // a) number of layers crossed
982 // b) track quality
983 // 5. Sign clusters by tracks in decreasing order of track quality
984 // 6. Build AliTRDtrack out of seeding tracklets
985 // 7. Cook MC label
986 // 8. Build ESD track and register it to the output list
987 //
988
e4f2f73d 989 AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
0906e73e 990 Int_t pars[4]; // MakeSeeds parameters
e4f2f73d 991
992 //Double_t alpha = AliTRDgeometry::GetAlpha();
993 //Double_t shift = .5 * alpha;
994 Int_t configs[kNConfigs];
995
996 // Build initial seeding configurations
997 Double_t quality = BuildSeedingConfigs(layer, configs);
998#ifdef DEBUG
999 if(AliTRDReconstructor::StreamLevel() > 1)
1000 AliInfo(Form("Plane config %d %d %d Quality %f"
1001 , configs[0], configs[1], configs[2], quality));
1002#endif
1003
1004 // Initialize contors
1005 Int_t ntracks, // number of TRD track candidates
1006 ntracks1, // number of registered TRD tracks/iter
1007 ntracks2 = 0; // number of all registered TRD tracks in stack
1008 fSieveSeeding = 0;
1009 do{
1010 // Loop over seeding configurations
1011 ntracks = 0; ntracks1 = 0;
1012 for (Int_t iconf = 0; iconf<3; iconf++) {
1013 pars[0] = configs[iconf];
1014 pars[1] = layer->GetStackNr();
1015 pars[2] = ntracks;
1016 ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
1017 if(ntracks == kMaxTracksStack) break;
1018 }
1019#ifdef DEBUG
0906e73e 1020 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in stack %d iteration %d.", ntracks, pars[1], fSieveSeeding));
e4f2f73d 1021#endif
1022 if(!ntracks) break;
1023
1024 // Sort the seeds according to their quality
1025 Int_t sort[kMaxTracksStack];
1026 TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
1027
1028 // Initialize number of tracks so far and logic switches
1029 Int_t ntracks0 = esdTrackList->GetEntriesFast();
1030 Bool_t signedTrack[kMaxTracksStack];
1031 Bool_t fakeTrack[kMaxTracksStack];
1032 for (Int_t i=0; i<ntracks; i++){
1033 signedTrack[i] = kFALSE;
1034 fakeTrack[i] = kFALSE;
1035 }
1036 //AliInfo("Selecting track candidates ...");
1037
1038 // Sieve clusters in decreasing order of track quality
1039 Double_t trackParams[7];
1040// AliTRDseedV1 *lseed = 0x0;
1041 Int_t jSieve = 0, candidates;
1042 do{
1043 //AliInfo(Form("\t\tITER = %i ", jSieve));
1044
1045 // Check track candidates
1046 candidates = 0;
1047 for (Int_t itrack = 0; itrack < ntracks; itrack++) {
1048 Int_t trackIndex = sort[itrack];
1049 if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
1050
1051
1052 // Calculate track parameters from tracklets seeds
1053 Int_t labelsall[1000];
1054 Int_t nlabelsall = 0;
1055 Int_t naccepted = 0;
1056 Int_t ncl = 0;
1057 Int_t nused = 0;
1058 Int_t nlayers = 0;
1059 Int_t findable = 0;
1060 for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
1061 Int_t jseed = kNPlanes*trackIndex+jLayer;
1062 if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15)
1063 findable++;
1064
1065 if(!sseed[jseed].IsOK()) continue;
1066 sseed[jseed].UpdateUsed();
1067 ncl += sseed[jseed].GetN2();
1068 nused += sseed[jseed].GetNUsed();
1069 nlayers++;
1070
1071 // Cooking label
0906e73e 1072 for (Int_t itime = 0; itime < fTimeBinsPerPlane; itime++) {
e4f2f73d 1073 if(!sseed[jseed].IsUsable(itime)) continue;
1074 naccepted++;
1075 Int_t tindex = 0, ilab = 0;
1076 while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
1077 labelsall[nlabelsall++] = tindex;
1078 ilab++;
1079 }
1080 }
1081 }
1082 // Filter duplicated tracks
1083 if (nused > 30){
0906e73e 1084 printf("Skip %d nused %d\n", trackIndex, nused);
e4f2f73d 1085 fakeTrack[trackIndex] = kTRUE;
1086 continue;
1087 }
1088 if (Float_t(nused)/ncl >= .25){
0906e73e 1089 printf("Skip %d nused/ncl >= .25\n", trackIndex);
e4f2f73d 1090 fakeTrack[trackIndex] = kTRUE;
1091 continue;
1092 }
1093
1094 // Classify tracks
1095 Bool_t skip = kFALSE;
1096 switch(jSieve){
1097 case 0:
1098 if(nlayers < 6) {skip = kTRUE; break;}
1099 if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1100 break;
1101
1102 case 1:
1103 if(nlayers < findable){skip = kTRUE; break;}
1104 if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
1105 break;
1106
1107 case 2:
1108 if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
1109 if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
1110 break;
1111
1112 case 3:
1113 if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1114 break;
1115
1116 case 4:
1117 if (nlayers == 3){skip = kTRUE; break;}
0906e73e 1118 //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
e4f2f73d 1119 break;
1120 }
1121 if(skip){
1122 candidates++;
0906e73e 1123 printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
e4f2f73d 1124 continue;
1125 }
1126 signedTrack[trackIndex] = kTRUE;
1127
1128
1129 // Build track label - what happens if measured data ???
1130 Int_t labels[1000];
1131 Int_t outlab[1000];
1132 Int_t nlab = 0;
1133 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1134 Int_t jseed = kNPlanes*trackIndex+iLayer;
1135 if(!sseed[jseed].IsOK()) continue;
1136 for(int ilab=0; ilab<2; ilab++){
1137 if(sseed[jseed].GetLabels(ilab) < 0) continue;
1138 labels[nlab] = sseed[jseed].GetLabels(ilab);
1139 nlab++;
1140 }
1141 }
1142 Freq(nlab,labels,outlab,kFALSE);
1143 Int_t label = outlab[0];
1144 Int_t frequency = outlab[1];
1145 Freq(nlabelsall,labelsall,outlab,kFALSE);
1146 Int_t label1 = outlab[0];
1147 Int_t label2 = outlab[2];
1148 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
1149
1150
1151 // Sign clusters
1152 AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
1153 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1154 Int_t jseed = kNPlanes*trackIndex+jLayer;
1155 if(!sseed[jseed].IsOK()) continue;
1156 if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
1157 sseed[jseed].UseClusters();
1158 if(!cl){
1159 Int_t ic = 0;
1160 while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
1161 clusterIndex = sseed[jseed].GetIndexes(ic);
1162 }
1163 }
1164 if(!cl) continue;
1165
1166
1167 // Build track parameters
1168 AliTRDseedV1 *lseed =&sseed[trackIndex*6];
1169 Int_t idx = 0;
1170 while(idx<3 && !lseed->IsOK()) {
1171 idx++;
1172 lseed++;
1173 }
1174 Double_t cR = lseed->GetC();
1175 trackParams[1] = lseed->GetYref(0);
1176 trackParams[2] = lseed->GetZref(0);
1177 trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
1178 trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
1179 trackParams[5] = cR;
1180 trackParams[0] = lseed->GetX0();
1181 trackParams[6] = layer[0].GetSector();/* *alpha+shift; // Supermodule*/
1182
1183#ifdef DEBUG
1184 if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
1185
1186 if(AliTRDReconstructor::StreamLevel() >= 1){
1187 Int_t sector = layer[0].GetSector();
1188 Int_t nclusters = 0;
1189 AliTRDseedV1 *dseed[6];
1190 for(int is=0; is<6; is++){
0906e73e 1191 dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
1192 dseed[is]->SetOwner();
e4f2f73d 1193 nclusters += sseed[is].GetN2();
1194 //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
1195 }
1196 //Int_t eventNrInFile = esd->GetEventNumberInFile();
1197 //AliInfo(Form("Number of clusters %d.", nclusters));
0906e73e 1198 TTreeSRedirector &cstreamer = *fDebugStreamer;
e4f2f73d 1199 cstreamer << "Clusters2TracksStack"
1200 << "Iter=" << fSieveSeeding
1201 << "Like=" << fTrackQuality[trackIndex]
1202 << "S0.=" << dseed[0]
1203 << "S1.=" << dseed[1]
1204 << "S2.=" << dseed[2]
1205 << "S3.=" << dseed[3]
1206 << "S4.=" << dseed[4]
1207 << "S5.=" << dseed[5]
1208 << "p0=" << trackParams[0]
1209 << "p1=" << trackParams[1]
1210 << "p2=" << trackParams[2]
1211 << "p3=" << trackParams[3]
1212 << "p4=" << trackParams[4]
1213 << "p5=" << trackParams[5]
1214 << "p6=" << trackParams[6]
1215 << "Sector=" << sector
1216 << "Stack=" << pars[1]
1217 << "Label=" << label
1218 << "Label1=" << label1
1219 << "Label2=" << label2
1220 << "FakeRatio=" << fakeratio
1221 << "Freq=" << frequency
1222 << "Ncl=" << ncl
1223 << "NLayers=" << nlayers
1224 << "Findable=" << findable
1225 << "NUsed=" << nused
1226 << "\n";
1227 //???for(int is=0; is<6; is++) delete dseed[is];
1228 }
1229#endif
1230
0906e73e 1231 AliTRDtrackV1 *track = AliTRDtrackerV1::MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
e4f2f73d 1232 if(!track){
1233 AliWarning("Fail to build a TRD Track.");
1234 continue;
1235 }
0906e73e 1236 AliInfo("End of MakeTrack()");
e4f2f73d 1237 AliESDtrack esdTrack;
1238 esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
1239 esdTrack.SetLabel(track->GetLabel());
1240 new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
1241 ntracks1++;
1242 }
1243
1244 jSieve++;
1245 } while(jSieve<5 && candidates); // end track candidates sieve
1246 if(!ntracks1) break;
1247
1248 // increment counters
1249 ntracks2 += ntracks1;
1250 fSieveSeeding++;
1251
1252 // Rebuild plane configurations and indices taking only unused clusters into account
1253 quality = BuildSeedingConfigs(layer, configs);
1254 //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
1255
0906e73e 1256 for(Int_t il = 0; il < kNPlanes * fTimeBinsPerPlane; il++) layer[il].BuildIndices(fSieveSeeding);
e4f2f73d 1257
1258#ifdef DEBUG
1259 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
1260#endif
1261 } while(fSieveSeeding<10); // end stack clusters sieve
1262
1263
1264
1265 //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
1266
1267 return ntracks2;
1268}
1269
1270//___________________________________________________________________
1271Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
1272 , Int_t *configs)
1273{
1274 //
1275 // Assign probabilities to chambers according to their
1276 // capability of producing seeds.
1277 //
1278 // Parameters :
1279 //
1280 // layers : Array of stack propagation layers for all 6 chambers in one stack
1281 // configs : On exit array of configuration indexes (see GetSeedingConfig()
1282 // for details) in the decreasing order of their seeding probabilities.
1283 //
1284 // Output :
1285 //
1286 // Return top configuration quality
1287 //
1288 // Detailed description:
1289 //
1290 // To each chamber seeding configuration (see GetSeedingConfig() for
1291 // the list of all configurations) one defines 2 quality factors:
1292 // - an apriori topological quality (see GetSeedingConfig() for details) and
1293 // - a data quality based on the uniformity of the distribution of
1294 // clusters over the x range (time bins population). See CookChamberQA() for details.
1295 // The overall chamber quality is given by the product of this 2 contributions.
1296 //
1297
e4f2f73d 1298 Double_t chamberQA[kNPlanes];
1299 for(int iplane=0; iplane<kNPlanes; iplane++){
0906e73e 1300 chamberQA[iplane] = MakeSeedingPlanes(&layers[iplane*fTimeBinsPerPlane]);
e4f2f73d 1301 //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
1302 }
1303
1304 Double_t tconfig[kNConfigs];
1305 Int_t planes[4];
1306 for(int iconf=0; iconf<kNConfigs; iconf++){
1307 GetSeedingConfig(iconf, planes);
d76231c8 1308 tconfig[iconf] = fgTopologicQA[iconf];
e4f2f73d 1309 for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]];
1310 }
1311
1312 TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
1313 return tconfig[configs[0]];
1314}
1315
1316//____________________________________________________________________
1317Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
1318 , AliTRDseedV1 *sseed
1319 , Int_t *ipar)
1320{
1321 //
1322 // Make tracklet seeds in the TRD stack.
1323 //
1324 // Parameters :
1325 // layers : Array of stack propagation layers containing clusters
1326 // sseed : Array of empty tracklet seeds. On exit they are filled.
1327 // ipar : Control parameters:
1328 // ipar[0] -> seeding chambers configuration
1329 // ipar[1] -> stack index
1330 // ipar[2] -> number of track candidates found so far
1331 //
1332 // Output :
1333 // Number of tracks candidates found.
1334 //
1335 // Detailed description
1336 //
1337 // The following steps are performed:
1338 // 1. Select seeding layers from seeding chambers
1339 // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
1340 // The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
1341 // this order. The parameters controling the range of accepted clusters in
1342 // layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
1343 // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
1344 // 4. Initialize seeding tracklets in the seeding chambers.
1345 // 5. Filter 0.
1346 // Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
1347 // Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
1348 // 6. Attach clusters to seeding tracklets and find linear approximation of
1349 // the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
1350 // clusters used by current seeds should not exceed ... (25).
1351 // 7. Filter 1.
1352 // All 4 seeding tracklets should be correctly constructed (see
1353 // AliTRDseedV1::AttachClustersIter())
1354 // 8. Helix fit of the seeding tracklets
1355 // 9. Filter 2.
1356 // Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
1357 // 10. Extrapolation of the helix fit to the other 2 chambers:
1358 // a) Initialization of extrapolation tracklet with fit parameters
1359 // b) Helix fit of tracklets
1360 // c) Attach clusters and linear interpolation to extrapolated tracklets
1361 // d) Helix fit of tracklets
1362 // 11. Improve seeding tracklets quality by reassigning clusters.
1363 // See AliTRDtrackerV1::ImproveSeedQuality() for details.
1364 // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
1365 // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
1366 // 14. Cooking labels for tracklets. Should be done only for MC
1367 // 15. Register seeds.
1368 //
1369
e4f2f73d 1370 AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
1371 AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
1372 Int_t ncl, mcl; // working variable for looping over clusters
1373 Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
1374 // chi2 storage
1375 // chi2[0] = tracklet chi2 on the Z direction
1376 // chi2[1] = tracklet chi2 on the R direction
1377 Double_t chi2[4];
1378
1379
1380 // this should be data member of AliTRDtrack
1381 Double_t seedQuality[kMaxTracksStack];
1382
1383 // unpack control parameters
1384 Int_t config = ipar[0];
1385 Int_t istack = ipar[1];
1386 Int_t ntracks = ipar[2];
1387 Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);
1388#ifdef DEBUG
1389 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
1390#endif
1391
1392 // Init chambers geometry
0906e73e 1393 Int_t det, tbRange[6]; // time bins inside the detector geometry
e4f2f73d 1394 Double_t hL[kNPlanes]; // Tilting angle
1395 Float_t padlength[kNPlanes]; // pad lenghts
1396 AliTRDpadPlane *pp;
1397 for(int il=0; il<kNPlanes; il++){
1398 pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
1399 hL[il] = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
0906e73e 1400 padlength[il] = pp->GetLengthIPad();
1401 det = il; // to be fixed !!!!!
1402 tbRange[il] = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
1403 //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
e4f2f73d 1404 }
1405
1406 Double_t cond0[4], cond1[4], cond2[4];
1407 // make seeding layers (to be moved in Clusters2TracksStack)
1408 AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
0906e73e 1409 for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * fTimeBinsPerPlane], planes[isl]);
e4f2f73d 1410
1411
1412 // Start finding seeds
1413 Int_t icl = 0;
1414 while((c[3] = (*layer[3])[icl++])){
1415 if(!c[3]) continue;
1416 layer[0]->BuildCond(c[3], cond0, 0);
1417 layer[0]->GetClusters(cond0, index, ncl);
1418 Int_t jcl = 0;
1419 while(jcl<ncl) {
1420 c[0] = (*layer[0])[index[jcl++]];
1421 if(!c[0]) continue;
1422 Double_t dx = c[3]->GetX() - c[0]->GetX();
1423 Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
1424 Double_t phi = (c[3]->GetY() - c[0]->GetY())/dx;
1425 layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
1426 layer[1]->GetClusters(cond1, jndex, mcl);
1427
1428 Int_t kcl = 0;
1429 while(kcl<mcl) {
1430 c[1] = (*layer[1])[jndex[kcl++]];
1431 if(!c[1]) continue;
1432 layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
1433 c[2] = layer[2]->GetNearestCluster(cond2);
1434 if(!c[2]) continue;
1435
1436 //AliInfo("Seeding clusters found. Building seeds ...");
1437 //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
1438 for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
1439
1440 fFitter->Reset();
1441
1442 fFitter->FitRieman(c, kNSeedPlanes);
1443
1444 chi2[0] = 0.; chi2[1] = 0.;
1445 AliTRDseedV1 *tseed = 0x0;
1446 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
0906e73e 1447 Int_t jLayer = planes[iLayer];
1448 tseed = &cseed[jLayer];
e4f2f73d 1449 tseed->SetRecoParam(fRecoParam);
0906e73e 1450 tseed->SetPlane(jLayer);
1451 tseed->SetTilt(hL[jLayer]);
1452 tseed->SetPadLength(padlength[jLayer]);
1453 tseed->SetNTimeBinsRange(tbRange[jLayer]);
1454 tseed->SetX0(layer[iLayer]->GetX());//layers[jLayer*fTimeBinsPerPlane].GetX());
1455
1456 tseed->Init(fFitter->GetRiemanFitter());
1457 // temporary until new AttachClusters()
1458 tseed->SetX0(layers[(jLayer+1)*fTimeBinsPerPlane-1].GetX());
e4f2f73d 1459 chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
1460 chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
1461 }
1462
1463 Bool_t isFake = kFALSE;
1464 if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1465 if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1466 if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1467#ifdef DEBUG
1468 if(AliTRDReconstructor::StreamLevel() >= 2){
1469 Float_t yref[4], ycluster[4];
1470 for(int il=0; il<4; il++){
1471 tseed = &cseed[planes[il]];
1472 yref[il] = tseed->GetYref(0);
1473 ycluster[il] = c[il]->GetY();
1474 }
1475 Float_t threshold = .5;//1./(3. - sLayer);
1476 Int_t ll = c[3]->GetLabel(0);
0906e73e 1477 TTreeSRedirector &cs0 = *fDebugStreamer;
e4f2f73d 1478 cs0 << "MakeSeeds0"
1479 <<"isFake=" << isFake
1480 <<"label=" << ll
1481 <<"threshold=" << threshold
1482 <<"chi2=" << chi2[1]
1483 <<"yref0="<<yref[0]
1484 <<"yref1="<<yref[1]
1485 <<"yref2="<<yref[2]
1486 <<"yref3="<<yref[3]
1487 <<"ycluster0="<<ycluster[0]
1488 <<"ycluster1="<<ycluster[1]
1489 <<"ycluster2="<<ycluster[2]
1490 <<"ycluster3="<<ycluster[3]
1491 <<"\n";
1492 }
1493#endif
1494
1495 if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
1496 //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
1497 continue;
1498 }
1499 if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
1500 //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
1501 continue;
1502 }
1503 //AliInfo("Passed chi2 filter.");
1504
1505#ifdef DEBUG
1506 if(AliTRDReconstructor::StreamLevel() >= 2){
1507 Float_t minmax[2] = { -100.0, 100.0 };
1508 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1509 Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
1510 if (max < minmax[1]) minmax[1] = max;
1511 Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
1512 if (min > minmax[0]) minmax[0] = min;
1513 }
1514 Double_t xpos[4];
1515 for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
0906e73e 1516 TTreeSRedirector &cstreamer = *fDebugStreamer;
e4f2f73d 1517 cstreamer << "MakeSeeds1"
1518 << "isFake=" << isFake
1519 << "config=" << config
1520 << "Cl0.=" << c[0]
1521 << "Cl1.=" << c[1]
1522 << "Cl2.=" << c[2]
1523 << "Cl3.=" << c[3]
1524 << "X0=" << xpos[0] //layer[sLayer]->GetX()
1525 << "X1=" << xpos[1] //layer[sLayer + 1]->GetX()
1526 << "X2=" << xpos[2] //layer[sLayer + 2]->GetX()
1527 << "X3=" << xpos[3] //layer[sLayer + 3]->GetX()
1528 << "Y2exp=" << cond2[0]
1529 << "Z2exp=" << cond2[1]
1530 << "Chi2R=" << chi2[0]
1531 << "Chi2Z=" << chi2[1]
1532 << "Seed0.=" << &cseed[planes[0]]
1533 << "Seed1.=" << &cseed[planes[1]]
1534 << "Seed2.=" << &cseed[planes[2]]
1535 << "Seed3.=" << &cseed[planes[3]]
1536 << "Zmin=" << minmax[0]
1537 << "Zmax=" << minmax[1]
1538 << "\n" ;
1539 }
1540#endif
1541 // try attaching clusters to tracklets
1542 Int_t nUsedCl = 0;
1543 Int_t nlayers = 0;
1544 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
0906e73e 1545 Int_t jLayer = planes[iLayer];
1546 if(!cseed[jLayer].AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 5., kFALSE, c[iLayer])) continue;
1547 nUsedCl += cseed[jLayer].GetNUsed();
e4f2f73d 1548 if(nUsedCl > 25) break;
1549 nlayers++;
1550 }
1551 if(nlayers < kNSeedPlanes){
1552 //AliInfo("Failed updating all seeds.");
1553 continue;
1554 }
1555 // fit tracklets and cook likelihood
1556 chi2[0] = 0.; chi2[1] = 0.;
1557 fFitter->FitRieman(&cseed[0], &planes[0]);
1558 AliRieman *rim = fFitter->GetRiemanFitter();
1559 for(int iLayer=0; iLayer<4; iLayer++){
0906e73e 1560 cseed[planes[iLayer]].Init(rim);
e4f2f73d 1561 chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
1562 chi2[1] += cseed[planes[iLayer]].GetChi2Y();
1563 }
1564 Double_t chi2r = chi2[1], chi2z = chi2[0];
1565 Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
1566 if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
1567 //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
1568 continue;
1569 }
1570 //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
1571
1572
1573 // book preliminary results
1574 seedQuality[ntracks] = like;
1575 fSeedLayer[ntracks] = config;/*sLayer;*/
1576
1577 // attach clusters to the extrapolation seeds
1578 Int_t lextrap[2];
1579 GetExtrapolationConfig(config, lextrap);
1580 Int_t nusedf = 0; // debug value
1581 for(int iLayer=0; iLayer<2; iLayer++){
1582 Int_t jLayer = lextrap[iLayer];
1583
1584 // prepare extrapolated seed
1585 cseed[jLayer].Reset();
1586 cseed[jLayer].SetRecoParam(fRecoParam);
0906e73e 1587 cseed[jLayer].SetPlane(jLayer);
e4f2f73d 1588 cseed[jLayer].SetTilt(hL[jLayer]);
0906e73e 1589 cseed[jLayer].SetX0(layers[(jLayer +1) * fTimeBinsPerPlane-1].GetX());
1590 cseed[jLayer].SetPadLength(padlength[jLayer]);
1591 cseed[jLayer].SetNTimeBinsRange(tbRange[jLayer]);
1592 cseed[jLayer].Init(rim);
1593// AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*fTimeBinsPerPlane], &cseed[jLayer]);
1594// if(cd == 0x0) continue;
e4f2f73d 1595
1596 // fit extrapolated seed
1597 AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
1598 if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
1599 if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
1600 AliTRDseedV1 tseed = cseed[jLayer];
0906e73e 1601 if(!tseed.AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 1000.)) continue;
e4f2f73d 1602 cseed[jLayer] = tseed;
1603 nusedf += cseed[jLayer].GetNUsed(); // debug value
1604 AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
1605 }
1606 //AliInfo("Extrapolation done.");
1607
1608 ImproveSeedQuality(layers, cseed);
1609 //AliInfo("Improve seed quality done.");
1610
1611 nlayers = 0;
1612 Int_t nclusters = 0;
1613 Int_t findable = 0;
1614 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1615 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
1616 if (!cseed[iLayer].IsOK()) continue;
1617 nclusters += cseed[iLayer].GetN2();
1618 nlayers++;
1619 }
1620 if (nlayers < 3){
1621 //AliInfo("Failed quality check on seeds.");
1622 continue;
1623 }
1624
1625 // fit full track and cook likelihoods
1626 fFitter->FitRieman(&cseed[0]);
1627 Double_t chi2ZF = 0., chi2RF = 0.;
1628 for(int ilayer=0; ilayer<6; ilayer++){
0906e73e 1629 cseed[ilayer].Init(fFitter->GetRiemanFitter());
e4f2f73d 1630 if (!cseed[ilayer].IsOK()) continue;
1631 //tchi2 = cseed[ilayer].GetChi2Z();
1632 //printf("layer %d chi2 %e\n", ilayer, tchi2);
1633 chi2ZF += cseed[ilayer].GetChi2Z();
1634 chi2RF += cseed[ilayer].GetChi2Y();
1635 }
1636 chi2ZF /= TMath::Max((nlayers - 3.), 1.);
1637 chi2RF /= TMath::Max((nlayers - 3.), 1.);
1638
1639 // do the final track fitting
1640 fFitter->SetLayers(nlayers);
1641#ifdef DEBUG
0906e73e 1642 fFitter->SetDebugStream(fDebugStreamer);
e4f2f73d 1643#endif
1644 fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
1645 Double_t param[3];
1646 Double_t chi2[2];
1647 fFitter->GetHyperplaneFitResults(param);
1648 fFitter->GetHyperplaneFitChi2(chi2);
1649 //AliInfo("Hyperplane fit done\n");
1650
1651 // finalize tracklets
1652 Int_t labels[12];
1653 Int_t outlab[24];
1654 Int_t nlab = 0;
1655 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1656 if (!cseed[iLayer].IsOK()) continue;
1657
1658 if (cseed[iLayer].GetLabels(0) >= 0) {
1659 labels[nlab] = cseed[iLayer].GetLabels(0);
1660 nlab++;
1661 }
1662
1663 if (cseed[iLayer].GetLabels(1) >= 0) {
1664 labels[nlab] = cseed[iLayer].GetLabels(1);
1665 nlab++;
1666 }
1667 }
1668 Freq(nlab,labels,outlab,kFALSE);
1669 Int_t label = outlab[0];
1670 Int_t frequency = outlab[1];
1671 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1672 cseed[iLayer].SetFreq(frequency);
1673 cseed[iLayer].SetC(param[1]/*cR*/);
1674 cseed[iLayer].SetCC(param[0]/*cC*/);
1675 cseed[iLayer].SetChi2(chi2[0]);
1676 cseed[iLayer].SetChi2Z(chi2ZF);
1677 }
1678
1679#ifdef DEBUG
1680 if(AliTRDReconstructor::StreamLevel() >= 2){
1681 Double_t curv = (fFitter->GetRiemanFitter())->GetC();
0906e73e 1682 TTreeSRedirector &cstreamer = *fDebugStreamer;
e4f2f73d 1683 cstreamer << "MakeSeeds2"
1684 << "C=" << curv
1685 << "Chi2R=" << chi2r
1686 << "Chi2Z=" << chi2z
1687 << "Chi2TR=" << chi2[0]
1688 << "Chi2TC=" << chi2[1]
1689 << "Chi2RF=" << chi2RF
1690 << "Chi2ZF=" << chi2ZF
1691 << "Ncl=" << nclusters
1692 << "Nlayers=" << nlayers
1693 << "NUsedS=" << nUsedCl
1694 << "NUsed=" << nusedf
1695 << "Findable" << findable
1696 << "Like=" << like
1697 << "S0.=" << &cseed[0]
1698 << "S1.=" << &cseed[1]
1699 << "S2.=" << &cseed[2]
1700 << "S3.=" << &cseed[3]
1701 << "S4.=" << &cseed[4]
1702 << "S5.=" << &cseed[5]
1703 << "Label=" << label
1704 << "Freq=" << frequency
1705 << "\n";
1706 }
1707#endif
1708
1709 ntracks++;
1710 if(ntracks == kMaxTracksStack){
1711 AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
1712 return ntracks;
1713 }
1714 cseed += 6;
1715 }
1716 }
1717 }
1718 for(int isl=0; isl<4; isl++) delete layer[isl];
1719
1720 return ntracks;
1721}
1722
1723//_____________________________________________________________________________
0906e73e 1724AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
e4f2f73d 1725{
1726 //
1727 // Build a TRD track out of tracklet candidates
1728 //
1729 // Parameters :
1730 // seeds : array of tracklets
1731 // params : track parameters (see MakeSeeds() function body for a detailed description)
1732 //
1733 // Output :
1734 // The TRD track.
1735 //
1736 // Detailed description
1737 //
1738 // To be discussed with Marian !!
1739 //
1740
e4f2f73d 1741 Double_t alpha = AliTRDgeometry::GetAlpha();
1742 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
1743 Double_t c[15];
1744
1745 c[ 0] = 0.2;
1746 c[ 1] = 0.0; c[ 2] = 2.0;
1747 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
1748 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
1749 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
1750
0906e73e 1751 AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, &params[1], c, params[0], params[6]*alpha+shift);
e4f2f73d 1752 track->PropagateTo(params[0]-5.0);
1753 track->ResetCovariance(1);
0906e73e 1754 Int_t nc = FollowBackProlongation(*track);
1755 AliInfo(Form("N clusters for track %d", nc));
1756 if (nc < 30) {
e4f2f73d 1757 delete track;
0906e73e 1758 track = 0x0;
1759 } else {
1760// track->CookdEdx();
1761// track->CookdEdxTimBin(-1);
1762// CookLabel(track, 0.9);
e4f2f73d 1763 }
1764
1765 return track;
1766}
1767
0906e73e 1768//____________________________________________________________________
1769void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
1770{
1771 // to be implemented, preferably at the level of TRD tracklet. !!!!!!!
1772}
1773
e4f2f73d 1774//____________________________________________________________________
1775void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
1776 , AliTRDseedV1 *cseed)
1777{
1778 //
1779 // Sort tracklets according to "quality" and try to "improve" the first 4 worst
1780 //
1781 // Parameters :
1782 // layers : Array of propagation layers for a stack/supermodule
1783 // cseed : Array of 6 seeding tracklets which has to be improved
1784 //
1785 // Output :
1786 // cssed : Improved seeds
1787 //
1788 // Detailed description
1789 //
1790 // Iterative procedure in which new clusters are searched for each
1791 // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
1792 // can be maximized. If some optimization is found the old seeds are replaced.
1793 //
1794
1795 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1796 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1797
1798 // make a local working copy
1799 AliTRDseedV1 bseed[6];
1800 for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
1801
1802
1803 Float_t lastquality = 10000.0;
1804 Float_t lastchi2 = 10000.0;
1805 Float_t chi2 = 1000.0;
1806
1807 for (Int_t iter = 0; iter < 4; iter++) {
1808 Float_t sumquality = 0.0;
1809 Float_t squality[6];
1810 Int_t sortindexes[6];
1811
1812 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1813 squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
1814 sumquality +=squality[jLayer];
1815 }
1816 if ((sumquality >= lastquality) || (chi2 > lastchi2)) break;
1817
1818
1819 lastquality = sumquality;
1820 lastchi2 = chi2;
1821 if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
1822
1823
1824 TMath::Sort(6, squality, sortindexes, kFALSE);
1825 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1826 Int_t bLayer = sortindexes[jLayer];
1827 bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
1828 }
1829
1830 chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
1831 } // Loop: iter
1832}
1833
1834//____________________________________________________________________
0906e73e 1835Double_t AliTRDtrackerV1::MakeSeedingPlanes(AliTRDstackLayer *layers)
e4f2f73d 1836{
1837 //
1838 // Calculate plane quality for seeding.
1839 //
1840 //
1841 // Parameters :
1842 // layers : Array of propagation layers for this plane.
1843 //
1844 // Output :
1845 // plane quality factor for seeding
1846 //
1847 // Detailed description
1848 //
1849 // The quality of the plane for seeding is higher if:
1850 // 1. the average timebin population is closer to an integer number
1851 // 2. the distribution of clusters/timebin is closer to a uniform distribution.
1852 // - the slope of the first derivative of a parabolic fit is small or
1853 // - the slope of a linear fit is small
1854 //
1855
1856 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1857 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1858
1859// Double_t x;
1860// TLinearFitter fitter(1, "pol1");
1861// fitter.ClearPoints();
1862 Int_t ncl = 0;
1863 Int_t nused = 0;
1864 Int_t nClLayer;
1865 for(int itb=0; itb<nTimeBins; itb++){
1866 //x = layer[itb].GetX();
1867 //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
1868 //if(!layer[itb].GetNClusters()) continue;
1869 //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
1870 nClLayer = layers[itb].GetNClusters();
1871 ncl += nClLayer;
1872 for(Int_t incl = 0; incl < nClLayer; incl++)
1873 if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
1874 }
1875
1876 // calculate the deviation of the mean number of clusters from the
1877 // closest integer values
d76231c8 1878 Float_t nclMed = float(ncl-nused)/nTimeBins;
1879 Int_t ncli = Int_t(nclMed);
1880 Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
1881 nclDev -= (nclDev>.5) && ncli ? .5 : 0.;
1882 /*Double_t quality = */ return TMath::Exp(-2.*nclDev);
e4f2f73d 1883
1884// // get slope of the derivative
1885// if(!fitter.Eval()) return quality;
1886// fitter.PrintResults(3);
1887// Double_t a = fitter.GetParameter(1);
1888//
0906e73e 1889// printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);
e4f2f73d 1890// return quality*TMath::Exp(-a);
1891}
1892
1893//____________________________________________________________________
1894Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
1895 , Int_t planes[4]
1896 , Double_t *chi2)
1897{
1898 //
1899 // Calculate the probability of this track candidate.
1900 //
1901 // Parameters :
1902 // cseeds : array of candidate tracklets
1903 // planes : array of seeding planes (see seeding configuration)
1904 // chi2 : chi2 values (on the Z and Y direction) from the rieman fit of the track.
1905 //
1906 // Output :
1907 // likelihood value
1908 //
1909 // Detailed description
1910 //
1911 // The track quality is estimated based on the following 4 criteria:
1912 // 1. precision of the rieman fit on the Y direction (likea)
1913 // 2. chi2 on the Y direction (likechi2y)
1914 // 3. chi2 on the Z direction (likechi2z)
1915 // 4. number of attached clusters compared to a reference value
1916 // (see AliTRDrecoParam::fkFindable) (likeN)
1917 //
1918 // The distributions for each type of probabilities are given below as of
1919 // (date). They have to be checked to assure consistency of estimation.
1920 //
1921
1922 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1923 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1924 // ratio of the total number of clusters/track which are expected to be found by the tracker.
1925 Float_t fgFindable = fRecoParam->GetFindableClusters();
1926
1927
1928 Int_t nclusters = 0;
1929 Double_t sumda = 0.;
1930 for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
1931 Int_t jlayer = planes[ilayer];
1932 nclusters += cseed[jlayer].GetN2();
1933 sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
1934 }
1935 Double_t likea = TMath::Exp(-sumda*10.6);
1936 Double_t likechi2y = 0.0000000001;
1937 if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
1938 Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
1939 Int_t enc = Int_t(fgFindable*4.*nTimeBins); // Expected Number Of Clusters, normally 72
1940 Double_t likeN = TMath::Exp(-(enc - nclusters) * 0.19);
1941
1942 Double_t like = likea * likechi2y * likechi2z * likeN;
1943
1944#ifdef DEBUG
1945 //AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
1946 if(AliTRDReconstructor::StreamLevel() >= 2){
0906e73e 1947 TTreeSRedirector &cstreamer = *fDebugStreamer;
e4f2f73d 1948 cstreamer << "CookLikelihood"
1949 << "sumda=" << sumda
1950 << "chi0=" << chi2[0]
1951 << "chi1=" << chi2[1]
1952 << "likea=" << likea
1953 << "likechi2y=" << likechi2y
1954 << "likechi2z=" << likechi2z
1955 << "nclusters=" << nclusters
1956 << "likeN=" << likeN
1957 << "like=" << like
1958 << "\n";
1959 }
1960#endif
1961
1962 return like;
1963}
1964
1965//___________________________________________________________________
1966void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
1967 , Int_t *planes
1968 , Double_t *params)
1969{
1970 //
1971 // Determines the Mean number of clusters per layer.
1972 // Needed to determine good Seeding Layers
1973 //
1974 // Parameters:
1975 // - Array of AliTRDstackLayers
1976 // - Container for the params
1977 //
1978 // Detailed description
1979 //
1980 // Two Iterations:
1981 // In the first Iteration the mean is calculted using all layers.
1982 // After this, all layers outside the 1-sigma-region are rejected.
1983 // Then the mean value and the standard-deviation are calculted a second
1984 // time in order to select all layers in the 1-sigma-region as good-candidates.
1985 //
1986
1987 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1988 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1989
1990 Float_t mean = 0, stdev = 0;
1991 Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
1992 Int_t position = 0;
1993 memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
1994 memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
1995 Int_t nused = 0;
1996 for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
1997 for(Int_t ils = 0; ils < nTimeBins; ils++){
1998 position = planes[ipl]*nTimeBins + ils;
1999 ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
2000 nused = 0;
2001 for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
2002 if((layers[position].GetCluster(icl))->IsUsed()) nused++;
2003 ncl[ipl * nTimeBins + ils] -= nused;
2004 }
2005 }
2006 // Declaration of quartils:
2007 //Double_t qvals[3] = {0.0, 0.0, 0.0};
2008 //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
2009 // Iterations
2010 Int_t counter;
2011 Double_t *array;
2012 Int_t *limit;
2013 Int_t nLayers = nTimeBins * kNSeedPlanes;
2014 for(Int_t iter = 0; iter < 2; iter++){
2015 array = (iter == 0) ? &ncl[0] : &mcl[0];
2016 limit = (iter == 0) ? &nLayers : &counter;
2017 counter = 0;
2018 if(iter == 1){
2019 for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
2020 if((ncl[i] > mean + stdev) || (ncl[i] < mean - stdev)) continue; // Outside 1-sigma region
2021// if((ncl[i] > qvals[2]) || (ncl[i] < qvals[0])) continue; // Outside 1-sigma region
2022 if(ncl[i] == 0) continue; // 0-Layers also rejected
2023 mcl[counter] = ncl[i];
2024 counter++;
2025 }
2026 }
2027 if(*limit == 0) break;
2028 printf("Limit = %d\n", *limit);
2029 //using quartils instead of mean and RMS
2030// TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
2031 mean = TMath::Median(*limit, array, 0x0);
2032 stdev = TMath::RMS(*limit, array);
2033 }
2034// printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
2035// memcpy(params,qvals,sizeof(Double_t)*3);
2036 params[1] = (Double_t)TMath::Nint(mean);
2037 params[0] = (Double_t)TMath::Nint(mean - stdev);
2038 params[2] = (Double_t)TMath::Nint(mean + stdev);
2039
2040}
2041
2042//___________________________________________________________________
2043Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
2044 , Double_t *params)
2045{
2046 //
2047 // Algorithm to find optimal seeding layer
2048 // Layers inside one sigma region (given by Quantiles) are sorted
2049 // according to their difference.
2050 // All layers outside are sorted according t
2051 //
2052 // Parameters:
2053 // - Array of AliTRDstackLayers (in the current plane !!!)
2054 // - Container for the Indices of the seeding Layer candidates
2055 //
2056 // Output:
2057 // - Number of Layers inside the 1-sigma-region
2058 //
2059 // The optimal seeding layer should contain the mean number of
2060 // custers in the layers in one chamber.
2061 //
2062
2063 //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
2064 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2065 const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
2066 Int_t nTimeBins = cal->GetNumberOfTimeBins();
2067 Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
2068 memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
2069 memset(indices, 0, sizeof(Int_t)*kNTimeBins);
2070 memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
2071 Int_t nused = 0;
2072 for(Int_t ils = 0; ils < nTimeBins; ils++){
2073 ncl[ils] = layers[ils].GetNClusters();
2074 nused = 0;
2075 for(Int_t icl = 0; icl < ncl[ils]; icl++)
2076 if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
2077 ncl[ils] -= nused;
2078 }
2079
2080 Float_t mean = params[1];
2081 for(Int_t ils = 0; ils < nTimeBins; ils++){
2082 memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
2083 indices[bins[ncl[ils]+1]] = ils;
2084 for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
2085 bins[i]++;
2086 }
2087
2088 //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
2089 Int_t sbin = -1;
2090 Int_t nElements;
2091 Int_t position = 0;
2092 TRandom *r = new TRandom();
2093 Int_t iter = 0;
2094 while(1){
2095 while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
2096 // Randomly selecting one bin
2097 sbin = (Int_t)r->Poisson(mean);
2098 }
2099 printf("Bin = %d\n",sbin);
2100 //Randomly selecting one Layer in the bin
2101 nElements = bins[sbin + 1] - bins[sbin];
2102 printf("nElements = %d\n", nElements);
2103 if(iter == 5){
2104 position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
2105 break;
2106 }
2107 else if(nElements==0){
2108 iter++;
2109 continue;
2110 }
2111 position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
2112 break;
2113 }
2114 delete r;
2115 return indices[position];
2116}
2117
2118//____________________________________________________________________
2119AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
2120 , AliTRDseedV1/*AliRieman*/ *reference)
2121{
2122 //
2123 // Finds a seeding Cluster for the extrapolation chamber.
2124 //
2125 // The seeding cluster should be as close as possible to the assumed
2126 // track which is represented by a Rieman fit.
2127 // Therefore the selecting criterion is the minimum distance between
2128 // the best fitting cluster and the Reference which is derived from
2129 // the AliTRDseed. Because all layers are assumed to be equally good
2130 // a linear search is performed.
2131 //
2132 // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
2133 // - sfit: the reference
2134 //
2135 // Output: - the best seeding cluster
2136 //
2137
2138 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2139 Int_t nTimeBins = cal->GetNumberOfTimeBins();
2140
2141 // distances as squared distances
2142 Int_t index = 0;
d76231c8 2143 Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0;
e4f2f73d 2144 ypos = reference->GetYref(0);
2145 zpos = reference->GetZref(0);
2146 AliTRDcluster *currentBest = 0x0, *temp = 0x0;
2147 for(Int_t ils = 0; ils < nTimeBins; ils++){
2148 // Reference positions
2149// ypos = reference->GetYat(layers[ils].GetX());
2150// zpos = reference->GetZat(layers[ils].GetX());
2151 index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
2152 if(index == -1) continue;
2153 temp = layers[ils].GetCluster(index);
2154 if(!temp) continue;
2155 distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
d76231c8 2156 if(distance < nearestDistance){
2157 nearestDistance = distance;
e4f2f73d 2158 currentBest = temp;
2159 }
2160 }
2161 return currentBest;
2162}
2163
2164//____________________________________________________________________
2165AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
2166 , Int_t plane)
2167{
2168 //
2169 // Creates a seeding layer
2170 //
2171
2172 // constants
2173 const Int_t kMaxRows = 16;
2174 const Int_t kMaxCols = 144;
2175 const Int_t kMaxPads = 2304;
2176
2177 // Get the calculation
2178 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2179 Int_t nTimeBins = cal->GetNumberOfTimeBins();
2180
2181 // Get the geometrical data of the chamber
2182 AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
2183 Int_t nCols = pp->GetNcols();
2184 Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
2185 Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
2186 Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
2187 Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
2188 Int_t nRows = pp->GetNrows();
2189 Float_t binlength = (ymax - ymin)/nCols;
2190 //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
2191
2192 // Fill the histogram
2193 Int_t arrpos;
2194 Float_t ypos;
2195 Int_t irow, nClusters;
2196 Int_t *histogram[kMaxRows]; // 2D-Histogram
2197 Int_t hvals[kMaxPads]; memset(hvals, 0, sizeof(Int_t)*kMaxPads);
2198 Float_t *sigmas[kMaxRows];
2199 Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads);
2200 AliTRDcluster *c = 0x0;
2201 for(Int_t irs = 0; irs < kMaxRows; irs++){
2202 histogram[irs] = &hvals[irs*kMaxCols];
2203 sigmas[irs] = &svals[irs*kMaxCols];
2204 }
2205 for(Int_t iTime = 0; iTime < nTimeBins; iTime++){
2206 nClusters = layers[iTime].GetNClusters();
2207 for(Int_t incl = 0; incl < nClusters; incl++){
2208 c = layers[iTime].GetCluster(incl);
2209 ypos = c->GetY();
2210 if(ypos > ymax && ypos < ymin) continue;
2211 irow = pp->GetPadRowNumber(c->GetZ()); // Zbin
2212 if(irow < 0)continue;
2213 arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
2214 if(ypos == ymax) arrpos = nCols - 1;
2215 histogram[irow][arrpos]++;
2216 sigmas[irow][arrpos] += c->GetSigmaZ2();
2217 }
2218 }
2219
2220// Now I have everything in the histogram, do the selection
2221// printf("Starting the analysis\n");
2222 //Int_t nPads = nCols * nRows;
2223 // This is what we are interested in: The center of gravity of the best candidates
2224 Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
2225 Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
2226 Float_t *cogy[kMaxRows];
2227 Float_t *cogz[kMaxRows];
2228 // Lookup-Table storing coordinates according ti the bins
2229 Float_t yLengths[kMaxCols];
2230 Float_t zLengths[kMaxRows];
2231 for(Int_t icnt = 0; icnt < nCols; icnt++){
2232 yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
2233 }
2234 for(Int_t icnt = 0; icnt < nRows; icnt++){
2235 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
2236 }
2237
2238 // A bitfield is used to mask the pads as usable
2239 Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
2240 for(UChar_t icount = 0; icount < nRows; icount++){
2241 cogy[icount] = &cogyvals[icount*kMaxCols];
2242 cogz[icount] = &cogzvals[icount*kMaxCols];
2243 }
2244 // In this array the array position of the best candidates will be stored
2245 Int_t cand[kMaxTracksStack];
2246 Float_t sigcands[kMaxTracksStack];
2247
2248 // helper variables
2249 Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
2250 Int_t nCandidates = 0;
2251 Float_t norm, cogv;
2252 // histogram filled -> Select best bins
2253 TMath::Sort(kMaxPads, hvals, indices); // bins storing a 0 should not matter
2254 // Set Threshold
2255 Int_t maximum = hvals[indices[0]]; // best
2256 Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
2257 Int_t col, row, lower, lower1, upper, upper1;
2258 for(Int_t ib = 0; ib < kMaxPads; ib++){
2259 if(nCandidates >= kMaxTracksStack){
2260 AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack));
2261 break;
2262 }
2263 // Positions
2264 row = indices[ib]/nCols;
2265 col = indices[ib]%nCols;
2266 // here will be the threshold condition:
2267 if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
2268 if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
2269 break; // number of clusters below threshold: break;
2270 }
2271 // passing: Mark the neighbors
2272 lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
2273 lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
2274 for(Int_t ic = lower; ic < upper; ++ic)
2275 for(Int_t ir = lower1; ir < upper1; ++ir){
2276 if(ic == col && ir == row) continue;
2277 mask[ic] |= (1 << ir);
2278 }
2279 // Storing the position in an array
2280 // testing for neigboring
2281 cogv = 0;
2282 norm = 0;
2283 lower = TMath::Max(col - 1,0);
2284 upper = TMath::Min(col + 2, nCols);
2285 for(Int_t inb = lower; inb < upper; ++inb){
2286 cogv += yLengths[inb] * histogram[row][inb];
2287 norm += histogram[row][inb];
2288 }
2289 cogy[row][col] = cogv / norm;
2290 cogv = 0; norm = 0;
2291 lower = TMath::Max(row - 1, 0);
2292 upper = TMath::Min(row + 2, nRows);
2293 for(Int_t inb = lower; inb < upper; ++inb){
2294 cogv += zLengths[inb] * histogram[inb][col];
2295 norm += histogram[inb][col];
2296 }
2297 cogz[row][col] = cogv / norm;
2298 // passed the filter
2299 cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array
2300 sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
2301 // Analysis output
2302 nCandidates++;
2303 }
2304 AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
2305 fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
2306 fakeLayer->SetSector(layers[0].GetSector());
2307 AliTRDcluster **fakeClusters = 0x0;
2308 UInt_t *fakeIndices = 0x0;
2309 if(nCandidates){
2310 fakeClusters = new AliTRDcluster*[nCandidates];
2311 fakeIndices = new UInt_t[nCandidates];
2312 UInt_t fakeIndex = 0;
2313 for(Int_t ican = 0; ican < nCandidates; ican++){
2314 fakeClusters[ican] = new AliTRDcluster();
2315 fakeClusters[ican]->SetX(fakeLayer->GetX());
2316 fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
2317 fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
2318 fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
2319 fakeIndices[ican] = fakeIndex++;// fantasy number
2320 }
2321 }
2322 fakeLayer->SetRecoParam(fRecoParam);
2323 fakeLayer->SetClustersArray(fakeClusters, nCandidates);
2324 fakeLayer->SetIndexArray(fakeIndices);
2325 fakeLayer->SetNRows(nRows);
2326 fakeLayer->BuildIndices();
2327 //fakeLayer->PrintClusters();
2328
2329#ifdef DEBUG
2330 if(AliTRDReconstructor::StreamLevel() >= 3){
0906e73e 2331 TMatrixD hist(nRows, nCols);
e4f2f73d 2332 for(Int_t i = 0; i < nRows; i++)
2333 for(Int_t j = 0; j < nCols; j++)
2334 hist(i,j) = histogram[i][j];
0906e73e 2335 TTreeSRedirector &cstreamer = *fDebugStreamer;
e4f2f73d 2336 cstreamer << "MakeSeedingLayer"
2337 << "Iteration=" << fSieveSeeding
2338 << "plane=" << plane
2339 << "ymin=" << ymin
2340 << "ymax=" << ymax
2341 << "zmin=" << zmin
2342 << "zmax=" << zmax
2343 << "L.=" << fakeLayer
2344 << "Histogram.=" << &hist
2345 << "\n";
2346 }
2347#endif
2348 return fakeLayer;
2349}
2350
2351//____________________________________________________________________
0906e73e 2352void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
e4f2f73d 2353{
2354 //
2355 // Map seeding configurations to detector planes.
2356 //
2357 // Parameters :
2358 // iconfig : configuration index
2359 // planes : member planes of this configuration. On input empty.
2360 //
2361 // Output :
2362 // planes : contains the planes which are defining the configuration
2363 //
2364 // Detailed description
2365 //
2366 // Here is the list of seeding planes configurations together with
2367 // their topological classification:
2368 //
2369 // 0 - 5432 TQ 0
2370 // 1 - 4321 TQ 0
2371 // 2 - 3210 TQ 0
2372 // 3 - 5321 TQ 1
2373 // 4 - 4210 TQ 1
2374 // 5 - 5431 TQ 1
2375 // 6 - 4320 TQ 1
2376 // 7 - 5430 TQ 2
2377 // 8 - 5210 TQ 2
2378 // 9 - 5421 TQ 3
2379 // 10 - 4310 TQ 3
2380 // 11 - 5410 TQ 4
2381 // 12 - 5420 TQ 5
2382 // 13 - 5320 TQ 5
2383 // 14 - 5310 TQ 5
2384 //
2385 // The topologic quality is modeled as follows:
2386 // 1. The general model is define by the equation:
2387 // p(conf) = exp(-conf/2)
2388 // 2. According to the topologic classification, configurations from the same
2389 // class are assigned the agerage value over the model values.
2390 // 3. Quality values are normalized.
2391 //
2392 // The topologic quality distribution as function of configuration is given below:
2393 //Begin_Html
2394 // <img src="gif/topologicQA.gif">
2395 //End_Html
2396 //
2397
2398 switch(iconfig){
2399 case 0: // 5432 TQ 0
2400 planes[0] = 2;
2401 planes[1] = 3;
2402 planes[2] = 4;
2403 planes[3] = 5;
2404 break;
2405 case 1: // 4321 TQ 0
2406 planes[0] = 1;
2407 planes[1] = 2;
2408 planes[2] = 3;
2409 planes[3] = 4;
2410 break;
2411 case 2: // 3210 TQ 0
2412 planes[0] = 0;
2413 planes[1] = 1;
2414 planes[2] = 2;
2415 planes[3] = 3;
2416 break;
2417 case 3: // 5321 TQ 1
2418 planes[0] = 1;
2419 planes[1] = 2;
2420 planes[2] = 3;
2421 planes[3] = 5;
2422 break;
2423 case 4: // 4210 TQ 1
2424 planes[0] = 0;
2425 planes[1] = 1;
2426 planes[2] = 2;
2427 planes[3] = 4;
2428 break;
2429 case 5: // 5431 TQ 1
2430 planes[0] = 1;
2431 planes[1] = 3;
2432 planes[2] = 4;
2433 planes[3] = 5;
2434 break;
2435 case 6: // 4320 TQ 1
2436 planes[0] = 0;
2437 planes[1] = 2;
2438 planes[2] = 3;
2439 planes[3] = 4;
2440 break;
2441 case 7: // 5430 TQ 2
2442 planes[0] = 0;
2443 planes[1] = 3;
2444 planes[2] = 4;
2445 planes[3] = 5;
2446 break;
2447 case 8: // 5210 TQ 2
2448 planes[0] = 0;
2449 planes[1] = 1;
2450 planes[2] = 2;
2451 planes[3] = 5;
2452 break;
2453 case 9: // 5421 TQ 3
2454 planes[0] = 1;
2455 planes[1] = 2;
2456 planes[2] = 4;
2457 planes[3] = 5;
2458 break;
2459 case 10: // 4310 TQ 3
2460 planes[0] = 0;
2461 planes[1] = 1;
2462 planes[2] = 3;
2463 planes[3] = 4;
2464 break;
2465 case 11: // 5410 TQ 4
2466 planes[0] = 0;
2467 planes[1] = 1;
2468 planes[2] = 4;
2469 planes[3] = 5;
2470 break;
2471 case 12: // 5420 TQ 5
2472 planes[0] = 0;
2473 planes[1] = 2;
2474 planes[2] = 4;
2475 planes[3] = 5;
2476 break;
2477 case 13: // 5320 TQ 5
2478 planes[0] = 0;
2479 planes[1] = 2;
2480 planes[2] = 3;
2481 planes[3] = 5;
2482 break;
2483 case 14: // 5310 TQ 5
2484 planes[0] = 0;
2485 planes[1] = 1;
2486 planes[2] = 3;
2487 planes[3] = 5;
2488 break;
2489 }
2490}
2491
2492//____________________________________________________________________
0906e73e 2493void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
e4f2f73d 2494{
2495 //
2496 // Returns the extrapolation planes for a seeding configuration.
2497 //
2498 // Parameters :
2499 // iconfig : configuration index
2500 // planes : planes which are not in this configuration. On input empty.
2501 //
2502 // Output :
2503 // planes : contains the planes which are not in the configuration
2504 //
2505 // Detailed description
2506 //
2507
2508 switch(iconfig){
2509 case 0: // 5432 TQ 0
2510 planes[0] = 1;
2511 planes[1] = 0;
2512 break;
2513 case 1: // 4321 TQ 0
2514 planes[0] = 5;
2515 planes[1] = 0;
2516 break;
2517 case 2: // 3210 TQ 0
2518 planes[0] = 4;
2519 planes[1] = 5;
2520 break;
2521 case 3: // 5321 TQ 1
2522 planes[0] = 4;
2523 planes[1] = 0;
2524 break;
2525 case 4: // 4210 TQ 1
2526 planes[0] = 5;
2527 planes[1] = 3;
2528 break;
2529 case 5: // 5431 TQ 1
2530 planes[0] = 2;
2531 planes[1] = 0;
2532 break;
2533 case 6: // 4320 TQ 1
2534 planes[0] = 5;
2535 planes[1] = 1;
2536 break;
2537 case 7: // 5430 TQ 2
2538 planes[0] = 2;
2539 planes[1] = 1;
2540 break;
2541 case 8: // 5210 TQ 2
2542 planes[0] = 4;
2543 planes[1] = 3;
2544 break;
2545 case 9: // 5421 TQ 3
2546 planes[0] = 3;
2547 planes[1] = 0;
2548 break;
2549 case 10: // 4310 TQ 3
2550 planes[0] = 5;
2551 planes[1] = 2;
2552 break;
2553 case 11: // 5410 TQ 4
2554 planes[0] = 3;
2555 planes[1] = 2;
2556 break;
2557 case 12: // 5420 TQ 5
2558 planes[0] = 3;
2559 planes[1] = 1;
2560 break;
2561 case 13: // 5320 TQ 5
2562 planes[0] = 4;
2563 planes[1] = 1;
2564 break;
2565 case 14: // 5310 TQ 5
2566 planes[0] = 4;
2567 planes[1] = 2;
2568 break;
2569 }
2570}