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