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