git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
ALIfied Pythia version for PYQUEN.
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerV1.cxx
CommitLineData
e4f2f73d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// Track finder //
21// //
22// Authors: //
23// Alex Bercuci <A.Bercuci@gsi.de> //
24// Markus Fasel <M.Fasel@gsi.de> //
25// //
26///////////////////////////////////////////////////////////////////////////////
27
28#include <Riostream.h>
29#include <stdio.h>
30#include <string.h>
31
32#include <TBranch.h>
33#include <TFile.h>
34#include <TGraph.h>
35#include <TH1D.h>
36#include <TH2D.h>
37#include <TLinearFitter.h>
38#include <TObjArray.h>
39#include <TROOT.h>
40#include <TTree.h>
41#include <TClonesArray.h>
42#include <TRandom.h>
43#include <TTreeStream.h>
44
45#include "AliLog.h"
46#include "AliESDEvent.h"
47#include "AliAlignObj.h"
48#include "AliRieman.h"
49#include "AliTrackPointArray.h"
50
51#include "AliTRDtracker.h"
52#include "AliTRDtrackerV1.h"
53#include "AliTRDgeometry.h"
54#include "AliTRDpadPlane.h"
55#include "AliTRDgeometry.h"
56#include "AliTRDcluster.h"
57#include "AliTRDtrack.h"
58#include "AliTRDseed.h"
59#include "AliTRDcalibDB.h"
60#include "AliTRDCommonParam.h"
61#include "AliTRDReconstructor.h"
62#include "AliTRDCalibraFillHisto.h"
63#include "AliTRDtrackerFitter.h"
64#include "AliTRDstackLayer.h"
65#include "AliTRDrecoParam.h"
66#include "AliTRDseedV1.h"
67
68#define DEBUG
69
70ClassImp(AliTRDtrackerV1)
d76231c8 71Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
e4f2f73d 72 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
73 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
74 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
75};
76
77//____________________________________________________________________
78AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p)
79 :AliTRDtracker()
80 ,fSieveSeeding(0)
fbb2ea06 81 ,fRecoParam(p)
e4f2f73d 82 ,fFitter(0x0)
fbb2ea06 83 ,fDebugStreamerV1(0x0)
e4f2f73d 84{
85 //
86 // Default constructor. Nothing is initialized.
87 //
88
89}
90
91//____________________________________________________________________
92AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p)
93 :AliTRDtracker(in)
94 ,fSieveSeeding(0)
95 ,fRecoParam(p)
96 ,fFitter(0x0)
fbb2ea06 97 ,fDebugStreamerV1(0x0)
e4f2f73d 98{
99 //
100 // Standard constructor.
101 // Setting of the geometry file, debug output (if enabled)
102 // and initilize fitter helper.
103 //
104
105 fFitter = new AliTRDtrackerFitter();
106
107#ifdef DEBUG
fbb2ea06 108 fDebugStreamerV1 = new TTreeSRedirector("TRDdebug.root");
109 fFitter->SetDebugStream(fDebugStreamerV1);
e4f2f73d 110#endif
111
112}
113
114//____________________________________________________________________
115AliTRDtrackerV1::~AliTRDtrackerV1()
116{
117 //
118 // Destructor
119 //
120
fbb2ea06 121 if(fDebugStreamerV1) delete fDebugStreamerV1;
e4f2f73d 122 if(fFitter) delete fFitter;
123 if(fRecoParam) delete fRecoParam;
fbb2ea06 124
e4f2f73d 125}
126
127//____________________________________________________________________
128Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
129{
130 //
131 // Steering stand alone tracking for full TRD detector
132 //
133 // Parameters :
134 // esd : The ESD event. On output it contains
135 // the ESD tracks found in TRD.
136 //
137 // Output :
138 // Number of tracks found in the TRD detector.
139 //
140 // Detailed description
141 // 1. Launch individual SM trackers.
142 // See AliTRDtrackerV1::Clusters2TracksSM() for details.
143 //
144
145 if(!fRecoParam){
146 AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
147 return 0;
148 }
149
150 //AliInfo("Start Track Finder ...");
151 Int_t ntracks = 0;
152 for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
153 //AliInfo(Form("Processing supermodule %i ...", ism));
154 ntracks += Clusters2TracksSM(fTrSec[ism], esd);
155 }
156 AliInfo(Form("Found %d TRD tracks.", ntracks));
157 return ntracks;
158}
159
160//____________________________________________________________________
161Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
162 , AliESDEvent *esd)
163{
164 //
165 // Steer tracking for one SM.
166 //
167 // Parameters :
168 // sector : Array of (SM) propagation layers containing clusters
169 // esd : The current ESD event. On output it contains the also
170 // the ESD (TRD) tracks found in this SM.
171 //
172 // Output :
173 // Number of tracks found in this TRD supermodule.
174 //
175 // Detailed description
176 //
177 // 1. Unpack AliTRDpropagationLayers objects for each stack.
178 // 2. Launch stack tracking.
179 // See AliTRDtrackerV1::Clusters2TracksStack() for details.
180 // 3. Pack results in the ESD event.
181 //
182
183 AliTRDpadPlane *pp = 0x0;
184
185 // allocate space for esd tracks in this SM
186 TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
187 esdTrackList.SetOwner();
188 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
189 Int_t nTimeBins = cal->GetNumberOfTimeBins();
190 const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
191
192 Int_t ntracks = 0;
193 Int_t nClStack = 0;
194 for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
195 AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
196
197 nClStack = 0;
198 //AliInfo(Form("Processing stack %i ...",istack));
199 //AliInfo("Building stack propagation layers ...");
200 for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
201 pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
202 Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad()
203 + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
204 //Debug
205 Double_t z0 = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
fbb2ea06 206 const AliTRDpropagationLayer kSMlayer(*(sector->GetLayer(ilayer)));
207 stackLayer[ilayer] = kSMlayer;
e4f2f73d 208#ifdef DEBUG
fbb2ea06 209 stackLayer[ilayer].SetDebugStream(fDebugStreamerV1);
e4f2f73d 210#endif
211 stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
212 stackLayer[ilayer].SetSector(sector->GetSector());
213 stackLayer[ilayer].SetStackNr(istack);
214 stackLayer[ilayer].SetNRows(pp->GetNrows());
215 stackLayer[ilayer].SetRecoParam(fRecoParam);
216 stackLayer[ilayer].BuildIndices();
217 nClStack += stackLayer[ilayer].GetNClusters();
218 }
219 //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
220 if(nClStack < kFindable) continue;
221 ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
222 }
223 //AliInfo(Form("Found %d tracks in SM", ntracks));
224
225 for(int itrack=0; itrack<ntracks; itrack++)
226 esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
227
228 return ntracks;
229}
230
231//____________________________________________________________________
232Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
233 , TClonesArray *esdTrackList)
234{
235 //
236 // Make tracks in one TRD stack.
237 //
238 // Parameters :
239 // layer : Array of stack propagation layers containing clusters
240 // esdTrackList : Array of ESD tracks found by the stand alone tracker.
241 // On exit the tracks found in this stack are appended.
242 //
243 // Output :
244 // Number of tracks found in this stack.
245 //
246 // Detailed description
247 //
248 // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
249 // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations.
250 // See AliTRDtrackerV1::MakeSeeds() for more details.
251 // 3. Arrange track candidates in decreasing order of their quality
252 // 4. Classify tracks in 5 categories according to:
253 // a) number of layers crossed
254 // b) track quality
255 // 5. Sign clusters by tracks in decreasing order of track quality
256 // 6. Build AliTRDtrack out of seeding tracklets
257 // 7. Cook MC label
258 // 8. Build ESD track and register it to the output list
259 //
260
fbb2ea06 261 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
262 Int_t nTimeBins = cal->GetNumberOfTimeBins();
e4f2f73d 263 AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
fbb2ea06 264 Int_t pars[3]; // MakeSeeds parameters
e4f2f73d 265
266 //Double_t alpha = AliTRDgeometry::GetAlpha();
267 //Double_t shift = .5 * alpha;
268 Int_t configs[kNConfigs];
269
270 // Build initial seeding configurations
271 Double_t quality = BuildSeedingConfigs(layer, configs);
272#ifdef DEBUG
273 if(AliTRDReconstructor::StreamLevel() > 1)
274 AliInfo(Form("Plane config %d %d %d Quality %f"
275 , configs[0], configs[1], configs[2], quality));
276#endif
277
278 // Initialize contors
279 Int_t ntracks, // number of TRD track candidates
280 ntracks1, // number of registered TRD tracks/iter
281 ntracks2 = 0; // number of all registered TRD tracks in stack
282 fSieveSeeding = 0;
283 do{
284 // Loop over seeding configurations
285 ntracks = 0; ntracks1 = 0;
286 for (Int_t iconf = 0; iconf<3; iconf++) {
287 pars[0] = configs[iconf];
288 pars[1] = layer->GetStackNr();
289 pars[2] = ntracks;
290 ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
291 if(ntracks == kMaxTracksStack) break;
292 }
293#ifdef DEBUG
fbb2ea06 294 if(AliTRDReconstructor::StreamLevel() > 1)
295 AliInfo(Form("Candidate TRD tracks %d in stack %d.", ntracks, pars[1]));
e4f2f73d 296#endif
297 if(!ntracks) break;
298
299 // Sort the seeds according to their quality
300 Int_t sort[kMaxTracksStack];
301 TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
302
303 // Initialize number of tracks so far and logic switches
304 Int_t ntracks0 = esdTrackList->GetEntriesFast();
305 Bool_t signedTrack[kMaxTracksStack];
306 Bool_t fakeTrack[kMaxTracksStack];
307 for (Int_t i=0; i<ntracks; i++){
308 signedTrack[i] = kFALSE;
309 fakeTrack[i] = kFALSE;
310 }
311 //AliInfo("Selecting track candidates ...");
312
313 // Sieve clusters in decreasing order of track quality
314 Double_t trackParams[7];
315// AliTRDseedV1 *lseed = 0x0;
316 Int_t jSieve = 0, candidates;
317 do{
318 //AliInfo(Form("\t\tITER = %i ", jSieve));
319
320 // Check track candidates
321 candidates = 0;
322 for (Int_t itrack = 0; itrack < ntracks; itrack++) {
323 Int_t trackIndex = sort[itrack];
324 if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
325
326
327 // Calculate track parameters from tracklets seeds
328 Int_t labelsall[1000];
329 Int_t nlabelsall = 0;
330 Int_t naccepted = 0;
331 Int_t ncl = 0;
332 Int_t nused = 0;
333 Int_t nlayers = 0;
334 Int_t findable = 0;
335 for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
336 Int_t jseed = kNPlanes*trackIndex+jLayer;
337 if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15)
338 findable++;
339
340 if(!sseed[jseed].IsOK()) continue;
341 sseed[jseed].UpdateUsed();
342 ncl += sseed[jseed].GetN2();
343 nused += sseed[jseed].GetNUsed();
344 nlayers++;
345
346 // Cooking label
fbb2ea06 347 for (Int_t itime = 0; itime < nTimeBins; itime++) {
e4f2f73d 348 if(!sseed[jseed].IsUsable(itime)) continue;
349 naccepted++;
350 Int_t tindex = 0, ilab = 0;
351 while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
352 labelsall[nlabelsall++] = tindex;
353 ilab++;
354 }
355 }
356 }
357 // Filter duplicated tracks
358 if (nused > 30){
fbb2ea06 359 //printf("Skip nused %d\n", nused);
e4f2f73d 360 fakeTrack[trackIndex] = kTRUE;
361 continue;
362 }
363 if (Float_t(nused)/ncl >= .25){
fbb2ea06 364 //printf("Skip nused/ncl >= .25\n");
e4f2f73d 365 fakeTrack[trackIndex] = kTRUE;
366 continue;
367 }
368
369 // Classify tracks
370 Bool_t skip = kFALSE;
371 switch(jSieve){
372 case 0:
373 if(nlayers < 6) {skip = kTRUE; break;}
374 if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
375 break;
376
377 case 1:
378 if(nlayers < findable){skip = kTRUE; break;}
379 if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
380 break;
381
382 case 2:
383 if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
384 if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
385 break;
386
387 case 3:
388 if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
389 break;
390
391 case 4:
392 if (nlayers == 3){skip = kTRUE; break;}
fbb2ea06 393 if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
e4f2f73d 394 break;
395 }
396 if(skip){
397 candidates++;
fbb2ea06 398 //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
e4f2f73d 399 continue;
400 }
401 signedTrack[trackIndex] = kTRUE;
402
403
404 // Build track label - what happens if measured data ???
405 Int_t labels[1000];
406 Int_t outlab[1000];
407 Int_t nlab = 0;
408 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
409 Int_t jseed = kNPlanes*trackIndex+iLayer;
410 if(!sseed[jseed].IsOK()) continue;
411 for(int ilab=0; ilab<2; ilab++){
412 if(sseed[jseed].GetLabels(ilab) < 0) continue;
413 labels[nlab] = sseed[jseed].GetLabels(ilab);
414 nlab++;
415 }
416 }
417 Freq(nlab,labels,outlab,kFALSE);
418 Int_t label = outlab[0];
419 Int_t frequency = outlab[1];
420 Freq(nlabelsall,labelsall,outlab,kFALSE);
421 Int_t label1 = outlab[0];
422 Int_t label2 = outlab[2];
423 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
424
425
426 // Sign clusters
427 AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
428 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
429 Int_t jseed = kNPlanes*trackIndex+jLayer;
430 if(!sseed[jseed].IsOK()) continue;
431 if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
432 sseed[jseed].UseClusters();
433 if(!cl){
434 Int_t ic = 0;
435 while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
436 clusterIndex = sseed[jseed].GetIndexes(ic);
437 }
438 }
439 if(!cl) continue;
440
441
442 // Build track parameters
443 AliTRDseedV1 *lseed =&sseed[trackIndex*6];
444 Int_t idx = 0;
445 while(idx<3 && !lseed->IsOK()) {
446 idx++;
447 lseed++;
448 }
449 Double_t cR = lseed->GetC();
450 trackParams[1] = lseed->GetYref(0);
451 trackParams[2] = lseed->GetZref(0);
452 trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
453 trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
454 trackParams[5] = cR;
455 trackParams[0] = lseed->GetX0();
456 trackParams[6] = layer[0].GetSector();/* *alpha+shift; // Supermodule*/
457
458#ifdef DEBUG
459 if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
460
461 if(AliTRDReconstructor::StreamLevel() >= 1){
462 Int_t sector = layer[0].GetSector();
463 Int_t nclusters = 0;
464 AliTRDseedV1 *dseed[6];
465 for(int is=0; is<6; is++){
fbb2ea06 466 dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is], kTRUE);
e4f2f73d 467 nclusters += sseed[is].GetN2();
468 //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
469 }
470 //Int_t eventNrInFile = esd->GetEventNumberInFile();
471 //AliInfo(Form("Number of clusters %d.", nclusters));
fbb2ea06 472 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
e4f2f73d 473 cstreamer << "Clusters2TracksStack"
474 << "Iter=" << fSieveSeeding
475 << "Like=" << fTrackQuality[trackIndex]
476 << "S0.=" << dseed[0]
477 << "S1.=" << dseed[1]
478 << "S2.=" << dseed[2]
479 << "S3.=" << dseed[3]
480 << "S4.=" << dseed[4]
481 << "S5.=" << dseed[5]
482 << "p0=" << trackParams[0]
483 << "p1=" << trackParams[1]
484 << "p2=" << trackParams[2]
485 << "p3=" << trackParams[3]
486 << "p4=" << trackParams[4]
487 << "p5=" << trackParams[5]
488 << "p6=" << trackParams[6]
489 << "Sector=" << sector
490 << "Stack=" << pars[1]
491 << "Label=" << label
492 << "Label1=" << label1
493 << "Label2=" << label2
494 << "FakeRatio=" << fakeratio
495 << "Freq=" << frequency
496 << "Ncl=" << ncl
497 << "NLayers=" << nlayers
498 << "Findable=" << findable
499 << "NUsed=" << nused
500 << "\n";
501 //???for(int is=0; is<6; is++) delete dseed[is];
502 }
503#endif
504
fbb2ea06 505 AliTRDtrack *track = AliTRDtrackerV1::RegisterSeed(&sseed[trackIndex*kNPlanes], trackParams);
e4f2f73d 506 if(!track){
507 AliWarning("Fail to build a TRD Track.");
508 continue;
509 }
510 AliESDtrack esdTrack;
511 esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
512 esdTrack.SetLabel(track->GetLabel());
513 new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
514 ntracks1++;
515 }
516
517 jSieve++;
518 } while(jSieve<5 && candidates); // end track candidates sieve
519 if(!ntracks1) break;
520
521 // increment counters
522 ntracks2 += ntracks1;
523 fSieveSeeding++;
524
525 // Rebuild plane configurations and indices taking only unused clusters into account
526 quality = BuildSeedingConfigs(layer, configs);
527 //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
528
fbb2ea06 529 for(Int_t il = 0; il < kNPlanes * nTimeBins; il++) layer[il].BuildIndices(fSieveSeeding);
e4f2f73d 530
531#ifdef DEBUG
532 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
533#endif
534 } while(fSieveSeeding<10); // end stack clusters sieve
535
536
537
538 //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
539
540 return ntracks2;
541}
542
543//___________________________________________________________________
544Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
545 , Int_t *configs)
546{
547 //
548 // Assign probabilities to chambers according to their
549 // capability of producing seeds.
550 //
551 // Parameters :
552 //
553 // layers : Array of stack propagation layers for all 6 chambers in one stack
554 // configs : On exit array of configuration indexes (see GetSeedingConfig()
555 // for details) in the decreasing order of their seeding probabilities.
556 //
557 // Output :
558 //
559 // Return top configuration quality
560 //
561 // Detailed description:
562 //
563 // To each chamber seeding configuration (see GetSeedingConfig() for
564 // the list of all configurations) one defines 2 quality factors:
565 // - an apriori topological quality (see GetSeedingConfig() for details) and
566 // - a data quality based on the uniformity of the distribution of
567 // clusters over the x range (time bins population). See CookChamberQA() for details.
568 // The overall chamber quality is given by the product of this 2 contributions.
569 //
570
fbb2ea06 571 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
572 Int_t nTimeBins = cal->GetNumberOfTimeBins();
573
e4f2f73d 574 Double_t chamberQA[kNPlanes];
575 for(int iplane=0; iplane<kNPlanes; iplane++){
fbb2ea06 576 chamberQA[iplane] = CookPlaneQA(&layers[iplane*nTimeBins]);
e4f2f73d 577 //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
578 }
579
580 Double_t tconfig[kNConfigs];
581 Int_t planes[4];
582 for(int iconf=0; iconf<kNConfigs; iconf++){
583 GetSeedingConfig(iconf, planes);
d76231c8 584 tconfig[iconf] = fgTopologicQA[iconf];
e4f2f73d 585 for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]];
586 }
587
588 TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
589 return tconfig[configs[0]];
590}
591
592//____________________________________________________________________
593Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
594 , AliTRDseedV1 *sseed
595 , Int_t *ipar)
596{
597 //
598 // Make tracklet seeds in the TRD stack.
599 //
600 // Parameters :
601 // layers : Array of stack propagation layers containing clusters
602 // sseed : Array of empty tracklet seeds. On exit they are filled.
603 // ipar : Control parameters:
604 // ipar[0] -> seeding chambers configuration
605 // ipar[1] -> stack index
606 // ipar[2] -> number of track candidates found so far
607 //
608 // Output :
609 // Number of tracks candidates found.
610 //
611 // Detailed description
612 //
613 // The following steps are performed:
614 // 1. Select seeding layers from seeding chambers
615 // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
616 // The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
617 // this order. The parameters controling the range of accepted clusters in
618 // layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
619 // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
620 // 4. Initialize seeding tracklets in the seeding chambers.
621 // 5. Filter 0.
622 // Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
623 // Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
624 // 6. Attach clusters to seeding tracklets and find linear approximation of
625 // the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
626 // clusters used by current seeds should not exceed ... (25).
627 // 7. Filter 1.
628 // All 4 seeding tracklets should be correctly constructed (see
629 // AliTRDseedV1::AttachClustersIter())
630 // 8. Helix fit of the seeding tracklets
631 // 9. Filter 2.
632 // Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
633 // 10. Extrapolation of the helix fit to the other 2 chambers:
634 // a) Initialization of extrapolation tracklet with fit parameters
635 // b) Helix fit of tracklets
636 // c) Attach clusters and linear interpolation to extrapolated tracklets
637 // d) Helix fit of tracklets
638 // 11. Improve seeding tracklets quality by reassigning clusters.
639 // See AliTRDtrackerV1::ImproveSeedQuality() for details.
640 // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
641 // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
642 // 14. Cooking labels for tracklets. Should be done only for MC
643 // 15. Register seeds.
644 //
645
fbb2ea06 646 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
647 Int_t nTimeBins = cal->GetNumberOfTimeBins();
e4f2f73d 648 AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
649 AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
650 Int_t ncl, mcl; // working variable for looping over clusters
651 Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
652 // chi2 storage
653 // chi2[0] = tracklet chi2 on the Z direction
654 // chi2[1] = tracklet chi2 on the R direction
655 Double_t chi2[4];
656
657
658 // this should be data member of AliTRDtrack
659 Double_t seedQuality[kMaxTracksStack];
660
661 // unpack control parameters
662 Int_t config = ipar[0];
663 Int_t istack = ipar[1];
664 Int_t ntracks = ipar[2];
665 Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);
666#ifdef DEBUG
667 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
668#endif
669
670 // Init chambers geometry
671 Double_t hL[kNPlanes]; // Tilting angle
672 Float_t padlength[kNPlanes]; // pad lenghts
673 AliTRDpadPlane *pp;
674 for(int il=0; il<kNPlanes; il++){
675 pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
676 hL[il] = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
fbb2ea06 677 padlength[il] = 10.; //pp->GetLengthIPad();
e4f2f73d 678 }
679
680 Double_t cond0[4], cond1[4], cond2[4];
681 // make seeding layers (to be moved in Clusters2TracksStack)
682 AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
fbb2ea06 683 for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * nTimeBins], planes[isl]);
e4f2f73d 684
685
686 // Start finding seeds
687 Int_t icl = 0;
688 while((c[3] = (*layer[3])[icl++])){
689 if(!c[3]) continue;
690 layer[0]->BuildCond(c[3], cond0, 0);
691 layer[0]->GetClusters(cond0, index, ncl);
692 Int_t jcl = 0;
693 while(jcl<ncl) {
694 c[0] = (*layer[0])[index[jcl++]];
695 if(!c[0]) continue;
696 Double_t dx = c[3]->GetX() - c[0]->GetX();
697 Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
698 Double_t phi = (c[3]->GetY() - c[0]->GetY())/dx;
699 layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
700 layer[1]->GetClusters(cond1, jndex, mcl);
701
702 Int_t kcl = 0;
703 while(kcl<mcl) {
704 c[1] = (*layer[1])[jndex[kcl++]];
705 if(!c[1]) continue;
706 layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
707 c[2] = layer[2]->GetNearestCluster(cond2);
708 if(!c[2]) continue;
709
710 //AliInfo("Seeding clusters found. Building seeds ...");
711 //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
712 for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
713
714 fFitter->Reset();
715
716 fFitter->FitRieman(c, kNSeedPlanes);
717
718 chi2[0] = 0.; chi2[1] = 0.;
719 AliTRDseedV1 *tseed = 0x0;
720 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
fbb2ea06 721 tseed = &cseed[planes[iLayer]];
e4f2f73d 722 tseed->SetRecoParam(fRecoParam);
fbb2ea06 723 tseed->SetLayer(planes[iLayer]);
724 tseed->SetTilt(hL[planes[iLayer]]);
725 tseed->SetPadLength(TMath::Sqrt(c[iLayer]->GetSigmaZ2()*12));
726 tseed->SetX0(layer[iLayer]->GetX());
727
728 tseed->Update(fFitter->GetRiemanFitter());
e4f2f73d 729 chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
730 chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
731 }
732
733 Bool_t isFake = kFALSE;
734 if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
735 if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
736 if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
737#ifdef DEBUG
738 if(AliTRDReconstructor::StreamLevel() >= 2){
739 Float_t yref[4], ycluster[4];
740 for(int il=0; il<4; il++){
741 tseed = &cseed[planes[il]];
742 yref[il] = tseed->GetYref(0);
743 ycluster[il] = c[il]->GetY();
744 }
745 Float_t threshold = .5;//1./(3. - sLayer);
746 Int_t ll = c[3]->GetLabel(0);
fbb2ea06 747 TTreeSRedirector &cs0 = *fDebugStreamerV1;
e4f2f73d 748 cs0 << "MakeSeeds0"
749 <<"isFake=" << isFake
750 <<"label=" << ll
751 <<"threshold=" << threshold
752 <<"chi2=" << chi2[1]
753 <<"yref0="<<yref[0]
754 <<"yref1="<<yref[1]
755 <<"yref2="<<yref[2]
756 <<"yref3="<<yref[3]
757 <<"ycluster0="<<ycluster[0]
758 <<"ycluster1="<<ycluster[1]
759 <<"ycluster2="<<ycluster[2]
760 <<"ycluster3="<<ycluster[3]
761 <<"\n";
762 }
763#endif
764
765 if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
766 //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
767 continue;
768 }
769 if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
770 //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
771 continue;
772 }
773 //AliInfo("Passed chi2 filter.");
774
775#ifdef DEBUG
776 if(AliTRDReconstructor::StreamLevel() >= 2){
777 Float_t minmax[2] = { -100.0, 100.0 };
778 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
779 Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
780 if (max < minmax[1]) minmax[1] = max;
781 Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
782 if (min > minmax[0]) minmax[0] = min;
783 }
784 Double_t xpos[4];
785 for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
fbb2ea06 786 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
e4f2f73d 787 cstreamer << "MakeSeeds1"
788 << "isFake=" << isFake
789 << "config=" << config
790 << "Cl0.=" << c[0]
791 << "Cl1.=" << c[1]
792 << "Cl2.=" << c[2]
793 << "Cl3.=" << c[3]
794 << "X0=" << xpos[0] //layer[sLayer]->GetX()
795 << "X1=" << xpos[1] //layer[sLayer + 1]->GetX()
796 << "X2=" << xpos[2] //layer[sLayer + 2]->GetX()
797 << "X3=" << xpos[3] //layer[sLayer + 3]->GetX()
798 << "Y2exp=" << cond2[0]
799 << "Z2exp=" << cond2[1]
800 << "Chi2R=" << chi2[0]
801 << "Chi2Z=" << chi2[1]
802 << "Seed0.=" << &cseed[planes[0]]
803 << "Seed1.=" << &cseed[planes[1]]
804 << "Seed2.=" << &cseed[planes[2]]
805 << "Seed3.=" << &cseed[planes[3]]
806 << "Zmin=" << minmax[0]
807 << "Zmax=" << minmax[1]
808 << "\n" ;
809 }
810#endif
811 // try attaching clusters to tracklets
812 Int_t nUsedCl = 0;
813 Int_t nlayers = 0;
814 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
fbb2ea06 815 AliTRDseedV1 tseed = cseed[planes[iLayer]];
816 if(!tseed.AttachClustersIter(&layers[planes[iLayer]*nTimeBins], 5., kFALSE, c[iLayer])) continue;
817 cseed[planes[iLayer]] = tseed;
818 nUsedCl += cseed[planes[iLayer]].GetNUsed();
e4f2f73d 819 if(nUsedCl > 25) break;
820 nlayers++;
821 }
822 if(nlayers < kNSeedPlanes){
823 //AliInfo("Failed updating all seeds.");
824 continue;
825 }
826 // fit tracklets and cook likelihood
827 chi2[0] = 0.; chi2[1] = 0.;
828 fFitter->FitRieman(&cseed[0], &planes[0]);
829 AliRieman *rim = fFitter->GetRiemanFitter();
830 for(int iLayer=0; iLayer<4; iLayer++){
fbb2ea06 831 cseed[planes[iLayer]].Update(rim);
e4f2f73d 832 chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
833 chi2[1] += cseed[planes[iLayer]].GetChi2Y();
834 }
835 Double_t chi2r = chi2[1], chi2z = chi2[0];
836 Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
837 if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
838 //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
839 continue;
840 }
841 //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
842
843
844 // book preliminary results
845 seedQuality[ntracks] = like;
846 fSeedLayer[ntracks] = config;/*sLayer;*/
847
848 // attach clusters to the extrapolation seeds
849 Int_t lextrap[2];
850 GetExtrapolationConfig(config, lextrap);
851 Int_t nusedf = 0; // debug value
852 for(int iLayer=0; iLayer<2; iLayer++){
853 Int_t jLayer = lextrap[iLayer];
854
855 // prepare extrapolated seed
856 cseed[jLayer].Reset();
857 cseed[jLayer].SetRecoParam(fRecoParam);
fbb2ea06 858 cseed[jLayer].SetLayer(jLayer);
e4f2f73d 859 cseed[jLayer].SetTilt(hL[jLayer]);
fbb2ea06 860 cseed[jLayer].SetX0(layers[jLayer * nTimeBins + (nTimeBins/2)].GetX()); // ????????
861 //cseed[jLayer].SetPadLength(??????????);
862 cseed[jLayer].Update(rim);
863
864 AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*nTimeBins], &cseed[jLayer]);
865 if(cd == 0x0) continue;
866// if(!cd) continue;
867 cseed[jLayer].SetPadLength(TMath::Sqrt(cd->GetSigmaZ2() * 12.));
868 cseed[jLayer].SetX0(cd->GetX()); // reference defined by a seedingLayer which is defined by the x-coordinate of the layers inside
e4f2f73d 869
870 // fit extrapolated seed
871 AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
872 if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
873 if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
874 AliTRDseedV1 tseed = cseed[jLayer];
fbb2ea06 875 if(!tseed.AttachClustersIter(&layers[jLayer*nTimeBins], 1000.)) continue;
e4f2f73d 876 cseed[jLayer] = tseed;
877 nusedf += cseed[jLayer].GetNUsed(); // debug value
878 AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
879 }
880 //AliInfo("Extrapolation done.");
881
882 ImproveSeedQuality(layers, cseed);
883 //AliInfo("Improve seed quality done.");
884
885 nlayers = 0;
886 Int_t nclusters = 0;
887 Int_t findable = 0;
888 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
889 if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
890 if (!cseed[iLayer].IsOK()) continue;
891 nclusters += cseed[iLayer].GetN2();
892 nlayers++;
893 }
894 if (nlayers < 3){
895 //AliInfo("Failed quality check on seeds.");
896 continue;
897 }
898
899 // fit full track and cook likelihoods
900 fFitter->FitRieman(&cseed[0]);
901 Double_t chi2ZF = 0., chi2RF = 0.;
902 for(int ilayer=0; ilayer<6; ilayer++){
fbb2ea06 903 cseed[ilayer].Update(fFitter->GetRiemanFitter());
e4f2f73d 904 if (!cseed[ilayer].IsOK()) continue;
905 //tchi2 = cseed[ilayer].GetChi2Z();
906 //printf("layer %d chi2 %e\n", ilayer, tchi2);
907 chi2ZF += cseed[ilayer].GetChi2Z();
908 chi2RF += cseed[ilayer].GetChi2Y();
909 }
910 chi2ZF /= TMath::Max((nlayers - 3.), 1.);
911 chi2RF /= TMath::Max((nlayers - 3.), 1.);
912
913 // do the final track fitting
914 fFitter->SetLayers(nlayers);
915#ifdef DEBUG
fbb2ea06 916 fFitter->SetDebugStream(fDebugStreamerV1);
e4f2f73d 917#endif
918 fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
919 Double_t param[3];
920 Double_t chi2[2];
921 fFitter->GetHyperplaneFitResults(param);
922 fFitter->GetHyperplaneFitChi2(chi2);
923 //AliInfo("Hyperplane fit done\n");
924
925 // finalize tracklets
926 Int_t labels[12];
927 Int_t outlab[24];
928 Int_t nlab = 0;
929 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
930 if (!cseed[iLayer].IsOK()) continue;
931
932 if (cseed[iLayer].GetLabels(0) >= 0) {
933 labels[nlab] = cseed[iLayer].GetLabels(0);
934 nlab++;
935 }
936
937 if (cseed[iLayer].GetLabels(1) >= 0) {
938 labels[nlab] = cseed[iLayer].GetLabels(1);
939 nlab++;
940 }
941 }
942 Freq(nlab,labels,outlab,kFALSE);
943 Int_t label = outlab[0];
944 Int_t frequency = outlab[1];
945 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
946 cseed[iLayer].SetFreq(frequency);
947 cseed[iLayer].SetC(param[1]/*cR*/);
948 cseed[iLayer].SetCC(param[0]/*cC*/);
949 cseed[iLayer].SetChi2(chi2[0]);
950 cseed[iLayer].SetChi2Z(chi2ZF);
951 }
952
953#ifdef DEBUG
954 if(AliTRDReconstructor::StreamLevel() >= 2){
955 Double_t curv = (fFitter->GetRiemanFitter())->GetC();
fbb2ea06 956 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
e4f2f73d 957 cstreamer << "MakeSeeds2"
958 << "C=" << curv
959 << "Chi2R=" << chi2r
960 << "Chi2Z=" << chi2z
961 << "Chi2TR=" << chi2[0]
962 << "Chi2TC=" << chi2[1]
963 << "Chi2RF=" << chi2RF
964 << "Chi2ZF=" << chi2ZF
965 << "Ncl=" << nclusters
966 << "Nlayers=" << nlayers
967 << "NUsedS=" << nUsedCl
968 << "NUsed=" << nusedf
969 << "Findable" << findable
970 << "Like=" << like
971 << "S0.=" << &cseed[0]
972 << "S1.=" << &cseed[1]
973 << "S2.=" << &cseed[2]
974 << "S3.=" << &cseed[3]
975 << "S4.=" << &cseed[4]
976 << "S5.=" << &cseed[5]
977 << "Label=" << label
978 << "Freq=" << frequency
979 << "\n";
980 }
981#endif
982
983 ntracks++;
984 if(ntracks == kMaxTracksStack){
985 AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
986 return ntracks;
987 }
988 cseed += 6;
989 }
990 }
991 }
992 for(int isl=0; isl<4; isl++) delete layer[isl];
993
994 return ntracks;
995}
996
997//_____________________________________________________________________________
fbb2ea06 998AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params)
e4f2f73d 999{
1000 //
1001 // Build a TRD track out of tracklet candidates
1002 //
1003 // Parameters :
1004 // seeds : array of tracklets
1005 // params : track parameters (see MakeSeeds() function body for a detailed description)
1006 //
1007 // Output :
1008 // The TRD track.
1009 //
1010 // Detailed description
1011 //
1012 // To be discussed with Marian !!
1013 //
1014
fbb2ea06 1015 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1016 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1017
e4f2f73d 1018 Double_t alpha = AliTRDgeometry::GetAlpha();
1019 Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
1020 Double_t c[15];
1021
1022 c[ 0] = 0.2;
1023 c[ 1] = 0.0; c[ 2] = 2.0;
1024 c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
1025 c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
1026 c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
1027
fbb2ea06 1028 Int_t index = 0;
1029 AliTRDcluster *cl = 0;
1030
1031 for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
1032 if (seeds[ilayer].IsOK()) {
1033 for (Int_t itime = nTimeBins - 1; itime > 0; itime--) {
1034 if (seeds[ilayer].GetIndexes(itime) > 0) {
1035 index = seeds[ilayer].GetIndexes(itime);
1036 cl = seeds[ilayer].GetClusters(itime);
1037 break;
1038 }
1039 }
1040 }
1041 if (index > 0) {
1042 break;
1043 }
1044 }
1045 if (cl == 0) return 0;
1046 AliTRDtrack *track = new AliTRDtrack(cl
1047 ,index
1048 ,&params[1]
1049 ,c
1050 ,params[0]
1051 ,params[6]*alpha+shift);
1052 // SetCluster(cl, 0); // A. Bercuci
e4f2f73d 1053 track->PropagateTo(params[0]-5.0);
1054 track->ResetCovariance(1);
fbb2ea06 1055 Int_t rc = FollowBackProlongation(*track);
1056 if (rc < 30) {
e4f2f73d 1057 delete track;
fbb2ea06 1058 track = 0;
1059 }
1060 else {
1061 track->CookdEdx();
1062 track->CookdEdxTimBin(-1);
1063 CookLabel(track,0.9);
e4f2f73d 1064 }
1065
1066 return track;
1067}
1068
1069//____________________________________________________________________
1070void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
1071 , AliTRDseedV1 *cseed)
1072{
1073 //
1074 // Sort tracklets according to "quality" and try to "improve" the first 4 worst
1075 //
1076 // Parameters :
1077 // layers : Array of propagation layers for a stack/supermodule
1078 // cseed : Array of 6 seeding tracklets which has to be improved
1079 //
1080 // Output :
1081 // cssed : Improved seeds
1082 //
1083 // Detailed description
1084 //
1085 // Iterative procedure in which new clusters are searched for each
1086 // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
1087 // can be maximized. If some optimization is found the old seeds are replaced.
1088 //
1089
1090 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1091 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1092
1093 // make a local working copy
1094 AliTRDseedV1 bseed[6];
1095 for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
1096
1097
1098 Float_t lastquality = 10000.0;
1099 Float_t lastchi2 = 10000.0;
1100 Float_t chi2 = 1000.0;
1101
1102 for (Int_t iter = 0; iter < 4; iter++) {
1103 Float_t sumquality = 0.0;
1104 Float_t squality[6];
1105 Int_t sortindexes[6];
1106
1107 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1108 squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
1109 sumquality +=squality[jLayer];
1110 }
1111 if ((sumquality >= lastquality) || (chi2 > lastchi2)) break;
1112
1113
1114 lastquality = sumquality;
1115 lastchi2 = chi2;
1116 if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
1117
1118
1119 TMath::Sort(6, squality, sortindexes, kFALSE);
1120 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1121 Int_t bLayer = sortindexes[jLayer];
1122 bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
1123 }
1124
1125 chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
1126 } // Loop: iter
1127}
1128
1129//____________________________________________________________________
fbb2ea06 1130Double_t AliTRDtrackerV1::CookPlaneQA(AliTRDstackLayer *layers)
e4f2f73d 1131{
1132 //
1133 // Calculate plane quality for seeding.
1134 //
1135 //
1136 // Parameters :
1137 // layers : Array of propagation layers for this plane.
1138 //
1139 // Output :
1140 // plane quality factor for seeding
1141 //
1142 // Detailed description
1143 //
1144 // The quality of the plane for seeding is higher if:
1145 // 1. the average timebin population is closer to an integer number
1146 // 2. the distribution of clusters/timebin is closer to a uniform distribution.
1147 // - the slope of the first derivative of a parabolic fit is small or
1148 // - the slope of a linear fit is small
1149 //
1150
1151 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1152 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1153
1154// Double_t x;
1155// TLinearFitter fitter(1, "pol1");
1156// fitter.ClearPoints();
1157 Int_t ncl = 0;
1158 Int_t nused = 0;
1159 Int_t nClLayer;
1160 for(int itb=0; itb<nTimeBins; itb++){
1161 //x = layer[itb].GetX();
1162 //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
1163 //if(!layer[itb].GetNClusters()) continue;
1164 //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
1165 nClLayer = layers[itb].GetNClusters();
1166 ncl += nClLayer;
1167 for(Int_t incl = 0; incl < nClLayer; incl++)
1168 if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
1169 }
1170
1171 // calculate the deviation of the mean number of clusters from the
1172 // closest integer values
d76231c8 1173 Float_t nclMed = float(ncl-nused)/nTimeBins;
1174 Int_t ncli = Int_t(nclMed);
1175 Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
1176 nclDev -= (nclDev>.5) && ncli ? .5 : 0.;
1177 /*Double_t quality = */ return TMath::Exp(-2.*nclDev);
e4f2f73d 1178
1179// // get slope of the derivative
1180// if(!fitter.Eval()) return quality;
1181// fitter.PrintResults(3);
1182// Double_t a = fitter.GetParameter(1);
1183//
fbb2ea06 1184// printf("nclDev(%f) a(%f)\n", nclDev, a);
e4f2f73d 1185// return quality*TMath::Exp(-a);
1186}
1187
1188//____________________________________________________________________
1189Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
1190 , Int_t planes[4]
1191 , Double_t *chi2)
1192{
1193 //
1194 // Calculate the probability of this track candidate.
1195 //
1196 // Parameters :
1197 // cseeds : array of candidate tracklets
1198 // planes : array of seeding planes (see seeding configuration)
1199 // chi2 : chi2 values (on the Z and Y direction) from the rieman fit of the track.
1200 //
1201 // Output :
1202 // likelihood value
1203 //
1204 // Detailed description
1205 //
1206 // The track quality is estimated based on the following 4 criteria:
1207 // 1. precision of the rieman fit on the Y direction (likea)
1208 // 2. chi2 on the Y direction (likechi2y)
1209 // 3. chi2 on the Z direction (likechi2z)
1210 // 4. number of attached clusters compared to a reference value
1211 // (see AliTRDrecoParam::fkFindable) (likeN)
1212 //
1213 // The distributions for each type of probabilities are given below as of
1214 // (date). They have to be checked to assure consistency of estimation.
1215 //
1216
1217 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1218 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1219 // ratio of the total number of clusters/track which are expected to be found by the tracker.
1220 Float_t fgFindable = fRecoParam->GetFindableClusters();
1221
1222
1223 Int_t nclusters = 0;
1224 Double_t sumda = 0.;
1225 for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
1226 Int_t jlayer = planes[ilayer];
1227 nclusters += cseed[jlayer].GetN2();
1228 sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
1229 }
1230 Double_t likea = TMath::Exp(-sumda*10.6);
1231 Double_t likechi2y = 0.0000000001;
1232 if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
1233 Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
1234 Int_t enc = Int_t(fgFindable*4.*nTimeBins); // Expected Number Of Clusters, normally 72
1235 Double_t likeN = TMath::Exp(-(enc - nclusters) * 0.19);
1236
1237 Double_t like = likea * likechi2y * likechi2z * likeN;
1238
1239#ifdef DEBUG
1240 //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));
1241 if(AliTRDReconstructor::StreamLevel() >= 2){
fbb2ea06 1242 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
e4f2f73d 1243 cstreamer << "CookLikelihood"
1244 << "sumda=" << sumda
1245 << "chi0=" << chi2[0]
1246 << "chi1=" << chi2[1]
1247 << "likea=" << likea
1248 << "likechi2y=" << likechi2y
1249 << "likechi2z=" << likechi2z
1250 << "nclusters=" << nclusters
1251 << "likeN=" << likeN
1252 << "like=" << like
1253 << "\n";
1254 }
1255#endif
1256
1257 return like;
1258}
1259
1260//___________________________________________________________________
1261void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
1262 , Int_t *planes
1263 , Double_t *params)
1264{
1265 //
1266 // Determines the Mean number of clusters per layer.
1267 // Needed to determine good Seeding Layers
1268 //
1269 // Parameters:
1270 // - Array of AliTRDstackLayers
1271 // - Container for the params
1272 //
1273 // Detailed description
1274 //
1275 // Two Iterations:
1276 // In the first Iteration the mean is calculted using all layers.
1277 // After this, all layers outside the 1-sigma-region are rejected.
1278 // Then the mean value and the standard-deviation are calculted a second
1279 // time in order to select all layers in the 1-sigma-region as good-candidates.
1280 //
1281
1282 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1283 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1284
1285 Float_t mean = 0, stdev = 0;
1286 Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
1287 Int_t position = 0;
1288 memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
1289 memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
1290 Int_t nused = 0;
1291 for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
1292 for(Int_t ils = 0; ils < nTimeBins; ils++){
1293 position = planes[ipl]*nTimeBins + ils;
1294 ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
1295 nused = 0;
1296 for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
1297 if((layers[position].GetCluster(icl))->IsUsed()) nused++;
1298 ncl[ipl * nTimeBins + ils] -= nused;
1299 }
1300 }
1301 // Declaration of quartils:
1302 //Double_t qvals[3] = {0.0, 0.0, 0.0};
1303 //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
1304 // Iterations
1305 Int_t counter;
1306 Double_t *array;
1307 Int_t *limit;
1308 Int_t nLayers = nTimeBins * kNSeedPlanes;
1309 for(Int_t iter = 0; iter < 2; iter++){
1310 array = (iter == 0) ? &ncl[0] : &mcl[0];
1311 limit = (iter == 0) ? &nLayers : &counter;
1312 counter = 0;
1313 if(iter == 1){
1314 for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
1315 if((ncl[i] > mean + stdev) || (ncl[i] < mean - stdev)) continue; // Outside 1-sigma region
1316// if((ncl[i] > qvals[2]) || (ncl[i] < qvals[0])) continue; // Outside 1-sigma region
1317 if(ncl[i] == 0) continue; // 0-Layers also rejected
1318 mcl[counter] = ncl[i];
1319 counter++;
1320 }
1321 }
1322 if(*limit == 0) break;
1323 printf("Limit = %d\n", *limit);
1324 //using quartils instead of mean and RMS
1325// TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
1326 mean = TMath::Median(*limit, array, 0x0);
1327 stdev = TMath::RMS(*limit, array);
1328 }
1329// printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
1330// memcpy(params,qvals,sizeof(Double_t)*3);
1331 params[1] = (Double_t)TMath::Nint(mean);
1332 params[0] = (Double_t)TMath::Nint(mean - stdev);
1333 params[2] = (Double_t)TMath::Nint(mean + stdev);
1334
1335}
1336
1337//___________________________________________________________________
1338Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
1339 , Double_t *params)
1340{
1341 //
1342 // Algorithm to find optimal seeding layer
1343 // Layers inside one sigma region (given by Quantiles) are sorted
1344 // according to their difference.
1345 // All layers outside are sorted according t
1346 //
1347 // Parameters:
1348 // - Array of AliTRDstackLayers (in the current plane !!!)
1349 // - Container for the Indices of the seeding Layer candidates
1350 //
1351 // Output:
1352 // - Number of Layers inside the 1-sigma-region
1353 //
1354 // The optimal seeding layer should contain the mean number of
1355 // custers in the layers in one chamber.
1356 //
1357
1358 //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
1359 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1360 const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
1361 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1362 Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
1363 memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
1364 memset(indices, 0, sizeof(Int_t)*kNTimeBins);
1365 memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
1366 Int_t nused = 0;
1367 for(Int_t ils = 0; ils < nTimeBins; ils++){
1368 ncl[ils] = layers[ils].GetNClusters();
1369 nused = 0;
1370 for(Int_t icl = 0; icl < ncl[ils]; icl++)
1371 if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
1372 ncl[ils] -= nused;
1373 }
1374
1375 Float_t mean = params[1];
1376 for(Int_t ils = 0; ils < nTimeBins; ils++){
1377 memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
1378 indices[bins[ncl[ils]+1]] = ils;
1379 for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
1380 bins[i]++;
1381 }
1382
1383 //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
1384 Int_t sbin = -1;
1385 Int_t nElements;
1386 Int_t position = 0;
1387 TRandom *r = new TRandom();
1388 Int_t iter = 0;
1389 while(1){
1390 while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
1391 // Randomly selecting one bin
1392 sbin = (Int_t)r->Poisson(mean);
1393 }
1394 printf("Bin = %d\n",sbin);
1395 //Randomly selecting one Layer in the bin
1396 nElements = bins[sbin + 1] - bins[sbin];
1397 printf("nElements = %d\n", nElements);
1398 if(iter == 5){
1399 position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
1400 break;
1401 }
1402 else if(nElements==0){
1403 iter++;
1404 continue;
1405 }
1406 position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
1407 break;
1408 }
1409 delete r;
1410 return indices[position];
1411}
1412
1413//____________________________________________________________________
1414AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
1415 , AliTRDseedV1/*AliRieman*/ *reference)
1416{
1417 //
1418 // Finds a seeding Cluster for the extrapolation chamber.
1419 //
1420 // The seeding cluster should be as close as possible to the assumed
1421 // track which is represented by a Rieman fit.
1422 // Therefore the selecting criterion is the minimum distance between
1423 // the best fitting cluster and the Reference which is derived from
1424 // the AliTRDseed. Because all layers are assumed to be equally good
1425 // a linear search is performed.
1426 //
1427 // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
1428 // - sfit: the reference
1429 //
1430 // Output: - the best seeding cluster
1431 //
1432
1433 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1434 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1435
1436 // distances as squared distances
1437 Int_t index = 0;
d76231c8 1438 Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0;
e4f2f73d 1439 ypos = reference->GetYref(0);
1440 zpos = reference->GetZref(0);
1441 AliTRDcluster *currentBest = 0x0, *temp = 0x0;
1442 for(Int_t ils = 0; ils < nTimeBins; ils++){
1443 // Reference positions
1444// ypos = reference->GetYat(layers[ils].GetX());
1445// zpos = reference->GetZat(layers[ils].GetX());
1446 index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
1447 if(index == -1) continue;
1448 temp = layers[ils].GetCluster(index);
1449 if(!temp) continue;
1450 distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
d76231c8 1451 if(distance < nearestDistance){
1452 nearestDistance = distance;
e4f2f73d 1453 currentBest = temp;
1454 }
1455 }
1456 return currentBest;
1457}
1458
1459//____________________________________________________________________
1460AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
1461 , Int_t plane)
1462{
1463 //
1464 // Creates a seeding layer
1465 //
1466
1467 // constants
1468 const Int_t kMaxRows = 16;
1469 const Int_t kMaxCols = 144;
1470 const Int_t kMaxPads = 2304;
1471
1472 // Get the calculation
1473 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1474 Int_t nTimeBins = cal->GetNumberOfTimeBins();
1475
1476 // Get the geometrical data of the chamber
1477 AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
1478 Int_t nCols = pp->GetNcols();
1479 Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
1480 Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
1481 Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
1482 Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
1483 Int_t nRows = pp->GetNrows();
1484 Float_t binlength = (ymax - ymin)/nCols;
1485 //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
1486
1487 // Fill the histogram
1488 Int_t arrpos;
1489 Float_t ypos;
1490 Int_t irow, nClusters;
1491 Int_t *histogram[kMaxRows]; // 2D-Histogram
1492 Int_t hvals[kMaxPads]; memset(hvals, 0, sizeof(Int_t)*kMaxPads);
1493 Float_t *sigmas[kMaxRows];
1494 Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads);
1495 AliTRDcluster *c = 0x0;
1496 for(Int_t irs = 0; irs < kMaxRows; irs++){
1497 histogram[irs] = &hvals[irs*kMaxCols];
1498 sigmas[irs] = &svals[irs*kMaxCols];
1499 }
1500 for(Int_t iTime = 0; iTime < nTimeBins; iTime++){
1501 nClusters = layers[iTime].GetNClusters();
1502 for(Int_t incl = 0; incl < nClusters; incl++){
1503 c = layers[iTime].GetCluster(incl);
1504 ypos = c->GetY();
1505 if(ypos > ymax && ypos < ymin) continue;
1506 irow = pp->GetPadRowNumber(c->GetZ()); // Zbin
1507 if(irow < 0)continue;
1508 arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
1509 if(ypos == ymax) arrpos = nCols - 1;
1510 histogram[irow][arrpos]++;
1511 sigmas[irow][arrpos] += c->GetSigmaZ2();
1512 }
1513 }
1514
1515// Now I have everything in the histogram, do the selection
1516// printf("Starting the analysis\n");
1517 //Int_t nPads = nCols * nRows;
1518 // This is what we are interested in: The center of gravity of the best candidates
1519 Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
1520 Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
1521 Float_t *cogy[kMaxRows];
1522 Float_t *cogz[kMaxRows];
1523 // Lookup-Table storing coordinates according ti the bins
1524 Float_t yLengths[kMaxCols];
1525 Float_t zLengths[kMaxRows];
1526 for(Int_t icnt = 0; icnt < nCols; icnt++){
1527 yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
1528 }
1529 for(Int_t icnt = 0; icnt < nRows; icnt++){
1530 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
1531 }
1532
1533 // A bitfield is used to mask the pads as usable
1534 Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
1535 for(UChar_t icount = 0; icount < nRows; icount++){
1536 cogy[icount] = &cogyvals[icount*kMaxCols];
1537 cogz[icount] = &cogzvals[icount*kMaxCols];
1538 }
1539 // In this array the array position of the best candidates will be stored
1540 Int_t cand[kMaxTracksStack];
1541 Float_t sigcands[kMaxTracksStack];
1542
1543 // helper variables
1544 Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
1545 Int_t nCandidates = 0;
1546 Float_t norm, cogv;
1547 // histogram filled -> Select best bins
1548 TMath::Sort(kMaxPads, hvals, indices); // bins storing a 0 should not matter
1549 // Set Threshold
1550 Int_t maximum = hvals[indices[0]]; // best
1551 Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
1552 Int_t col, row, lower, lower1, upper, upper1;
1553 for(Int_t ib = 0; ib < kMaxPads; ib++){
1554 if(nCandidates >= kMaxTracksStack){
1555 AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack));
1556 break;
1557 }
1558 // Positions
1559 row = indices[ib]/nCols;
1560 col = indices[ib]%nCols;
1561 // here will be the threshold condition:
1562 if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
1563 if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
1564 break; // number of clusters below threshold: break;
1565 }
1566 // passing: Mark the neighbors
1567 lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
1568 lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
1569 for(Int_t ic = lower; ic < upper; ++ic)
1570 for(Int_t ir = lower1; ir < upper1; ++ir){
1571 if(ic == col && ir == row) continue;
1572 mask[ic] |= (1 << ir);
1573 }
1574 // Storing the position in an array
1575 // testing for neigboring
1576 cogv = 0;
1577 norm = 0;
1578 lower = TMath::Max(col - 1,0);
1579 upper = TMath::Min(col + 2, nCols);
1580 for(Int_t inb = lower; inb < upper; ++inb){
1581 cogv += yLengths[inb] * histogram[row][inb];
1582 norm += histogram[row][inb];
1583 }
1584 cogy[row][col] = cogv / norm;
1585 cogv = 0; norm = 0;
1586 lower = TMath::Max(row - 1, 0);
1587 upper = TMath::Min(row + 2, nRows);
1588 for(Int_t inb = lower; inb < upper; ++inb){
1589 cogv += zLengths[inb] * histogram[inb][col];
1590 norm += histogram[inb][col];
1591 }
1592 cogz[row][col] = cogv / norm;
1593 // passed the filter
1594 cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array
1595 sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
1596 // Analysis output
1597 nCandidates++;
1598 }
1599 AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
1600 fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
1601 fakeLayer->SetSector(layers[0].GetSector());
1602 AliTRDcluster **fakeClusters = 0x0;
1603 UInt_t *fakeIndices = 0x0;
1604 if(nCandidates){
1605 fakeClusters = new AliTRDcluster*[nCandidates];
1606 fakeIndices = new UInt_t[nCandidates];
1607 UInt_t fakeIndex = 0;
1608 for(Int_t ican = 0; ican < nCandidates; ican++){
1609 fakeClusters[ican] = new AliTRDcluster();
1610 fakeClusters[ican]->SetX(fakeLayer->GetX());
1611 fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
1612 fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
1613 fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
1614 fakeIndices[ican] = fakeIndex++;// fantasy number
1615 }
1616 }
1617 fakeLayer->SetRecoParam(fRecoParam);
1618 fakeLayer->SetClustersArray(fakeClusters, nCandidates);
1619 fakeLayer->SetIndexArray(fakeIndices);
1620 fakeLayer->SetNRows(nRows);
1621 fakeLayer->BuildIndices();
1622 //fakeLayer->PrintClusters();
1623
1624#ifdef DEBUG
1625 if(AliTRDReconstructor::StreamLevel() >= 3){
fbb2ea06 1626 TMatrixT<double> hist(nRows, nCols);
e4f2f73d 1627 for(Int_t i = 0; i < nRows; i++)
1628 for(Int_t j = 0; j < nCols; j++)
1629 hist(i,j) = histogram[i][j];
fbb2ea06 1630 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
e4f2f73d 1631 cstreamer << "MakeSeedingLayer"
1632 << "Iteration=" << fSieveSeeding
1633 << "plane=" << plane
1634 << "ymin=" << ymin
1635 << "ymax=" << ymax
1636 << "zmin=" << zmin
1637 << "zmax=" << zmax
1638 << "L.=" << fakeLayer
1639 << "Histogram.=" << &hist
1640 << "\n";
1641 }
1642#endif
1643 return fakeLayer;
1644}
1645
1646//____________________________________________________________________
fbb2ea06 1647void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig
1648 , Int_t planes[4]) const
e4f2f73d 1649{
1650 //
1651 // Map seeding configurations to detector planes.
1652 //
1653 // Parameters :
1654 // iconfig : configuration index
1655 // planes : member planes of this configuration. On input empty.
1656 //
1657 // Output :
1658 // planes : contains the planes which are defining the configuration
1659 //
1660 // Detailed description
1661 //
1662 // Here is the list of seeding planes configurations together with
1663 // their topological classification:
1664 //
1665 // 0 - 5432 TQ 0
1666 // 1 - 4321 TQ 0
1667 // 2 - 3210 TQ 0
1668 // 3 - 5321 TQ 1
1669 // 4 - 4210 TQ 1
1670 // 5 - 5431 TQ 1
1671 // 6 - 4320 TQ 1
1672 // 7 - 5430 TQ 2
1673 // 8 - 5210 TQ 2
1674 // 9 - 5421 TQ 3
1675 // 10 - 4310 TQ 3
1676 // 11 - 5410 TQ 4
1677 // 12 - 5420 TQ 5
1678 // 13 - 5320 TQ 5
1679 // 14 - 5310 TQ 5
1680 //
1681 // The topologic quality is modeled as follows:
1682 // 1. The general model is define by the equation:
1683 // p(conf) = exp(-conf/2)
1684 // 2. According to the topologic classification, configurations from the same
1685 // class are assigned the agerage value over the model values.
1686 // 3. Quality values are normalized.
1687 //
1688 // The topologic quality distribution as function of configuration is given below:
1689 //Begin_Html
1690 // <img src="gif/topologicQA.gif">
1691 //End_Html
1692 //
1693
1694 switch(iconfig){
1695 case 0: // 5432 TQ 0
1696 planes[0] = 2;
1697 planes[1] = 3;
1698 planes[2] = 4;
1699 planes[3] = 5;
1700 break;
1701 case 1: // 4321 TQ 0
1702 planes[0] = 1;
1703 planes[1] = 2;
1704 planes[2] = 3;
1705 planes[3] = 4;
1706 break;
1707 case 2: // 3210 TQ 0
1708 planes[0] = 0;
1709 planes[1] = 1;
1710 planes[2] = 2;
1711 planes[3] = 3;
1712 break;
1713 case 3: // 5321 TQ 1
1714 planes[0] = 1;
1715 planes[1] = 2;
1716 planes[2] = 3;
1717 planes[3] = 5;
1718 break;
1719 case 4: // 4210 TQ 1
1720 planes[0] = 0;
1721 planes[1] = 1;
1722 planes[2] = 2;
1723 planes[3] = 4;
1724 break;
1725 case 5: // 5431 TQ 1
1726 planes[0] = 1;
1727 planes[1] = 3;
1728 planes[2] = 4;
1729 planes[3] = 5;
1730 break;
1731 case 6: // 4320 TQ 1
1732 planes[0] = 0;
1733 planes[1] = 2;
1734 planes[2] = 3;
1735 planes[3] = 4;
1736 break;
1737 case 7: // 5430 TQ 2
1738 planes[0] = 0;
1739 planes[1] = 3;
1740 planes[2] = 4;
1741 planes[3] = 5;
1742 break;
1743 case 8: // 5210 TQ 2
1744 planes[0] = 0;
1745 planes[1] = 1;
1746 planes[2] = 2;
1747 planes[3] = 5;
1748 break;
1749 case 9: // 5421 TQ 3
1750 planes[0] = 1;
1751 planes[1] = 2;
1752 planes[2] = 4;
1753 planes[3] = 5;
1754 break;
1755 case 10: // 4310 TQ 3
1756 planes[0] = 0;
1757 planes[1] = 1;
1758 planes[2] = 3;
1759 planes[3] = 4;
1760 break;
1761 case 11: // 5410 TQ 4
1762 planes[0] = 0;
1763 planes[1] = 1;
1764 planes[2] = 4;
1765 planes[3] = 5;
1766 break;
1767 case 12: // 5420 TQ 5
1768 planes[0] = 0;
1769 planes[1] = 2;
1770 planes[2] = 4;
1771 planes[3] = 5;
1772 break;
1773 case 13: // 5320 TQ 5
1774 planes[0] = 0;
1775 planes[1] = 2;
1776 planes[2] = 3;
1777 planes[3] = 5;
1778 break;
1779 case 14: // 5310 TQ 5
1780 planes[0] = 0;
1781 planes[1] = 1;
1782 planes[2] = 3;
1783 planes[3] = 5;
1784 break;
1785 }
1786}
1787
1788//____________________________________________________________________
fbb2ea06 1789void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig
1790 , Int_t planes[2]) const
e4f2f73d 1791{
1792 //
1793 // Returns the extrapolation planes for a seeding configuration.
1794 //
1795 // Parameters :
1796 // iconfig : configuration index
1797 // planes : planes which are not in this configuration. On input empty.
1798 //
1799 // Output :
1800 // planes : contains the planes which are not in the configuration
1801 //
1802 // Detailed description
1803 //
1804
1805 switch(iconfig){
1806 case 0: // 5432 TQ 0
1807 planes[0] = 1;
1808 planes[1] = 0;
1809 break;
1810 case 1: // 4321 TQ 0
1811 planes[0] = 5;
1812 planes[1] = 0;
1813 break;
1814 case 2: // 3210 TQ 0
1815 planes[0] = 4;
1816 planes[1] = 5;
1817 break;
1818 case 3: // 5321 TQ 1
1819 planes[0] = 4;
1820 planes[1] = 0;
1821 break;
1822 case 4: // 4210 TQ 1
1823 planes[0] = 5;
1824 planes[1] = 3;
1825 break;
1826 case 5: // 5431 TQ 1
1827 planes[0] = 2;
1828 planes[1] = 0;
1829 break;
1830 case 6: // 4320 TQ 1
1831 planes[0] = 5;
1832 planes[1] = 1;
1833 break;
1834 case 7: // 5430 TQ 2
1835 planes[0] = 2;
1836 planes[1] = 1;
1837 break;
1838 case 8: // 5210 TQ 2
1839 planes[0] = 4;
1840 planes[1] = 3;
1841 break;
1842 case 9: // 5421 TQ 3
1843 planes[0] = 3;
1844 planes[1] = 0;
1845 break;
1846 case 10: // 4310 TQ 3
1847 planes[0] = 5;
1848 planes[1] = 2;
1849 break;
1850 case 11: // 5410 TQ 4
1851 planes[0] = 3;
1852 planes[1] = 2;
1853 break;
1854 case 12: // 5420 TQ 5
1855 planes[0] = 3;
1856 planes[1] = 1;
1857 break;
1858 case 13: // 5320 TQ 5
1859 planes[0] = 4;
1860 planes[1] = 1;
1861 break;
1862 case 14: // 5310 TQ 5
1863 planes[0] = 4;
1864 planes[1] = 2;
1865 break;
1866 }
1867}