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