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