]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New Calibration and Trigger classes included
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
CommitLineData
46d29e70 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
0fa7dfa7 16/* $Id$ */
bbf92647 17
029cd327 18///////////////////////////////////////////////////////////////////////////////
19// //
20// The standard TRD tracker //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
a2cb5b3d 24#include <Riostream.h>
46d29e70 25#include <TFile.h>
46d29e70 26#include <TBranch.h>
5443e65e 27#include <TTree.h>
c630aafd 28#include <TObjArray.h>
46d29e70 29
46d29e70 30#include "AliTRDgeometry.h"
5443e65e 31#include "AliTRDparameter.h"
a5cadd36 32#include "AliTRDpadPlane.h"
dde59437 33#include "AliTRDgeometryFull.h"
46d29e70 34#include "AliTRDcluster.h"
35#include "AliTRDtrack.h"
b7a0917f 36#include "AliESD.h"
46d29e70 37
3551db50 38#include "AliTRDcalibDB.h"
39#include "AliTRDCommonParam.h"
40
7ad19338 41#include "TTreeStream.h"
42#include "TGraph.h"
46d29e70 43#include "AliTRDtracker.h"
69b55c55 44#include "TLinearFitter.h"
45#include "AliRieman.h"
3551db50 46#include "AliTrackPointArray.h"
47#include "AliAlignObj.h"
8979685e 48
7ad19338 49//
46d29e70 50
51ClassImp(AliTRDtracker)
69b55c55 52ClassImp(AliTRDseed)
46d29e70 53
5443e65e 54
bbf92647 55
029cd327 56 const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5;
029cd327 57 const Float_t AliTRDtracker::fgkLabelFraction = 0.8;
029cd327 58 const Double_t AliTRDtracker::fgkMaxChi2 = 12.;
f6625211 59//
69b55c55 60 const Double_t AliTRDtracker::fgkOffset = -0.019;
61 const Double_t AliTRDtracker::fgkOffsetX = 0.26; // "time offset"
69b55c55 62 const Double_t AliTRDtracker::fgkCoef = 0.0106; // angular shift
7ad19338 63 const Double_t AliTRDtracker::fgkMean = 0.;
f6625211 64 const Double_t AliTRDtracker::fgkDriftCorrection = 1.000; // drift coefficient correction
7ad19338 65 const Double_t AliTRDtracker::fgkExB = 0.072; // ExB angle - for error parameterization
66
67
68// poscorrection = fgkCoef*(GetLocalTimeBin() - fgkMean)+fgkOffset;
69
9c9d2487 70
89f05372 71//____________________________________________________________________
72AliTRDtracker::AliTRDtracker():AliTracker(),
73 fGeom(0),
74 fPar(0),
75 fNclusters(0),
76 fClusters(0),
77 fNseeds(0),
78 fSeeds(0),
79 fNtracks(0),
80 fTracks(0),
89f05372 81 fTimeBinsPerPlane(0),
89f05372 82 fAddTRDseeds(kFALSE),
83 fNoTilt(kFALSE)
84{
b7a0917f 85 // Default constructor
86
89f05372 87 for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
88 for(Int_t j=0;j<5;j++)
89 for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
7ad19338 90 fDebugStreamer = 0;
89f05372 91}
46d29e70 92//____________________________________________________________________
c630aafd 93AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
46d29e70 94{
5443e65e 95 //
96 // Main constructor
97 //
b8dc2353 98
5443e65e 99 fAddTRDseeds = kFALSE;
5443e65e 100 fGeom = NULL;
b8dc2353 101 fNoTilt = kFALSE;
5443e65e 102
103 TDirectory *savedir=gDirectory;
104 TFile *in=(TFile*)geomfile;
105 if (!in->IsOpen()) {
106 printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
dde59437 107 printf(" FULL TRD geometry and DEFAULT TRD parameter will be used\n");
5443e65e 108 }
109 else {
110 in->cd();
c630aafd 111// in->ls();
5443e65e 112 fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
113 fPar = (AliTRDparameter*) in->Get("TRDparameter");
c630aafd 114// fGeom->Dump();
5443e65e 115 }
46d29e70 116
5443e65e 117 if(fGeom) {
118 // fTzero = geo->GetT0();
7c1698cb 119 // printf("Found geometry version %d on file \n", fGeom->IsVersion());
5443e65e 120 }
121 else {
c630aafd 122 printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
dde59437 123 fGeom = new AliTRDgeometryFull();
3c625a9b 124 fGeom->SetPHOShole();
125 fGeom->SetRICHhole();
c630aafd 126 }
127
128 if (!fPar) {
129 printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n");
130 printf("The DEFAULT TRD parameter will be used\n");
7ad19338 131 fPar = new AliTRDparameter("Pica","Vyjebana");
5443e65e 132 }
7ad19338 133 fPar = new AliTRDparameter("Pica","Vyjebana");
598156ef 134 fPar->Init();
5443e65e 135
136 savedir->cd();
46d29e70 137
5443e65e 138 // fGeom->SetT0(fTzero);
0a29d0f1 139
46d29e70 140 fNclusters = 0;
141 fClusters = new TObjArray(2000);
142 fNseeds = 0;
5443e65e 143 fSeeds = new TObjArray(2000);
46d29e70 144 fNtracks = 0;
5443e65e 145 fTracks = new TObjArray(1000);
a819a5f7 146
029cd327 147 for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
148 Int_t trS = CookSectorIndex(geomS);
149 fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS, fPar);
3c625a9b 150 for (Int_t icham=0;icham<AliTRDgeometry::kNcham; icham++){
151 fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
152 }
5443e65e 153 }
3551db50 154 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
7ad19338 155 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
156 // Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle());
029cd327 157 if(tiltAngle < 0.1) {
b8dc2353 158 fNoTilt = kTRUE;
159 }
160
b8dc2353 161
bbf92647 162
5443e65e 163 // calculate max gap on track
a819a5f7 164
5443e65e 165 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
166 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
a819a5f7 167
7ad19338 168 Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
3551db50 169 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
7ad19338 170
5443e65e 171 Int_t tbAmp = fPar->GetTimeBefore();
172 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
b3a5a838 173 if(kTRUE) maxAmp = 0; // intentional until we change the parameter class
5443e65e 174 Int_t tbDrift = fPar->GetTimeMax();
8979685e 175 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx)+4; // MI change - take also last time bins
a819a5f7 176
8979685e 177 tbDrift = TMath::Min(tbDrift,maxDrift);
178 tbAmp = TMath::Min(tbAmp,maxAmp);
46d29e70 179
5443e65e 180 fTimeBinsPerPlane = tbAmp + tbDrift;
46d29e70 181
7ad19338 182 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
0a29d0f1 183
9c9d2487 184 savedir->cd();
5443e65e 185}
46d29e70 186
5443e65e 187//___________________________________________________________________
188AliTRDtracker::~AliTRDtracker()
46d29e70 189{
029cd327 190 //
191 // Destructor of AliTRDtracker
192 //
193
89f05372 194 if (fClusters) {
195 fClusters->Delete();
196 delete fClusters;
197 }
198 if (fTracks) {
199 fTracks->Delete();
200 delete fTracks;
201 }
202 if (fSeeds) {
203 fSeeds->Delete();
204 delete fSeeds;
205 }
5443e65e 206 delete fGeom;
207 delete fPar;
0a29d0f1 208
029cd327 209 for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
210 delete fTrSec[geomS];
5443e65e 211 }
7ad19338 212 if (fDebugStreamer) {
213 //fDebugStreamer->Close();
214 delete fDebugStreamer;
215 }
5443e65e 216}
46d29e70 217
9c9d2487 218//_____________________________________________________________________
219
9c9d2487 220Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
221 //
222 // Rotates the track when necessary
223 //
224
225 Double_t alpha = AliTRDgeometry::GetAlpha();
226 Double_t y = track->GetY();
227 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
228
c630aafd 229 //Int_t ns = AliTRDgeometry::kNsect;
9c9d2487 230 //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
231
232 if (y > ymax) {
233 //s = (s+1) % ns;
234 if (!track->Rotate(alpha)) return kFALSE;
235 } else if (y <-ymax) {
236 //s = (s-1+ns) % ns;
237 if (!track->Rotate(-alpha)) return kFALSE;
238 }
239
240 return kTRUE;
241}
242
46e2d86c 243
7ad19338 244AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index){
46e2d86c 245 //
246 //try to find cluster in the backup list
247 //
248 AliTRDcluster * cl =0;
249 UInt_t *indexes = track->GetBackupIndexes();
250 for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
251 if (indexes[i]==0) break;
252 AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
253 if (!cli) break;
254 if (cli->GetLocalTimeBin()!=timebin) continue;
255 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
256 if (iplane==plane) {
257 cl = cli;
7ad19338 258 index = indexes[i];
46e2d86c 259 break;
260 }
261 }
262 return cl;
263}
264
3c625a9b 265
266Int_t AliTRDtracker::GetLastPlane(AliTRDtrack * track){
267 //
268 //return last updated plane
269 Int_t lastplane=0;
270 UInt_t *indexes = track->GetBackupIndexes();
271 for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
272 AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
273 if (!cli) break;
274 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
275 if (iplane>lastplane) {
276 lastplane = iplane;
277 }
278 }
279 return lastplane;
280}
c630aafd 281//___________________________________________________________________
282Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
283{
284 //
285 // Finds tracks within the TRD. The ESD event is expected to contain seeds
286 // at the outer part of the TRD. The seeds
287 // are found within the TRD if fAddTRDseeds is TRUE.
288 // The tracks are propagated to the innermost time bin
289 // of the TRD and the ESD event is updated
290 //
291
292 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 293 Float_t foundMin = fgkMinClustersInTrack * timeBins;
c630aafd 294 Int_t nseed = 0;
295 Int_t found = 0;
296 Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
297
298 Int_t n = event->GetNumberOfTracks();
299 for (Int_t i=0; i<n; i++) {
300 AliESDtrack* seed=event->GetTrack(i);
301 ULong_t status=seed->GetStatus();
302 if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
303 if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
304 nseed++;
7ad19338 305
c630aafd 306 AliTRDtrack* seed2 = new AliTRDtrack(*seed);
46e2d86c 307 //seed2->ResetCovariance();
c630aafd 308 AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
309 AliTRDtrack &t=*pt;
310 FollowProlongation(t, innerTB);
311 if (t.GetNumberOfClusters() >= foundMin) {
312 UseClusters(&t);
029cd327 313 CookLabel(pt, 1-fgkLabelFraction);
c630aafd 314 // t.CookdEdx();
315 }
316 found++;
317// cout<<found<<'\r';
318
319 if(PropagateToTPC(t)) {
320 seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
321 }
322 delete seed2;
323 delete pt;
324 }
325
326 cout<<"Number of loaded seeds: "<<nseed<<endl;
327 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
328
329 // after tracks from loaded seeds are found and the corresponding
330 // clusters are used, look for additional seeds from TRD
c630aafd 331
c630aafd 332
333 cout<<"Total number of found tracks: "<<found<<endl;
334
335 return 0;
336}
5443e65e 337
c630aafd 338
5443e65e 339
c630aafd 340//_____________________________________________________________________________
341Int_t AliTRDtracker::PropagateBack(AliESD* event) {
342 //
343 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
344 // backpropagated by the TPC tracker. Each seed is first propagated
345 // to the TRD, and then its prolongation is searched in the TRD.
346 // If sufficiently long continuation of the track is found in the TRD
347 // the track is updated, otherwise it's stored as originaly defined
348 // by the TPC tracker.
349 //
350
351 Int_t found=0;
c5a8e3df 352 Float_t foundMin = 20;
c630aafd 353 Int_t n = event->GetNumberOfTracks();
4f1c04d3 354 //
355 //Sort tracks
356 Float_t *quality =new Float_t[n];
357 Int_t *index =new Int_t[n];
c630aafd 358 for (Int_t i=0; i<n; i++) {
359 AliESDtrack* seed=event->GetTrack(i);
4f1c04d3 360 Double_t covariance[15];
361 seed->GetExternalCovariance(covariance);
362 quality[i] = covariance[0]+covariance[2];
363 }
364 TMath::Sort(n,quality,index,kFALSE);
365 //
366 for (Int_t i=0; i<n; i++) {
367 // AliESDtrack* seed=event->GetTrack(i);
368 AliESDtrack* seed=event->GetTrack(index[i]);
369
c630aafd 370 ULong_t status=seed->GetStatus();
371 if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
372 if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
373
374 Int_t lbl = seed->GetLabel();
375 AliTRDtrack *track = new AliTRDtrack(*seed);
376 track->SetSeedLabel(lbl);
f4e9508c 377 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); //make backup
c630aafd 378 fNseeds++;
4f1c04d3 379 Float_t p4 = track->GetC();
7bed16a7 380 //
f6625211 381 Int_t expectedClr = FollowBackProlongation(*track);
f4e9508c 382 /*
383 // only debug purpose
3c625a9b 384 if (track->GetNumberOfClusters()<expectedClr/3){
385 AliTRDtrack *track1 = new AliTRDtrack(*seed);
386 track1->SetSeedLabel(lbl);
387 FollowBackProlongation(*track1);
388 AliTRDtrack *track2= new AliTRDtrack(*seed);
389 track->SetSeedLabel(lbl);
390 FollowBackProlongation(*track2);
391 delete track1;
392 delete track2;
393 }
f4e9508c 394 */
4f1c04d3 395 if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)<0.2 || TMath::Abs(track->GetPt())>0.8 ) {
396 //
397 //make backup for back propagation
398 //
399 Int_t foundClr = track->GetNumberOfClusters();
400 if (foundClr >= foundMin) {
401 track->CookdEdx();
8979685e 402 CookdEdxTimBin(*track);
4f1c04d3 403 CookLabel(track, 1-fgkLabelFraction);
69b55c55 404 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
4f1c04d3 405 if(track->GetChi2()/track->GetNumberOfClusters()<4) { // sign only gold tracks
406 if (seed->GetKinkIndex(0)==0&&TMath::Abs(track->GetPt())<1.5 ) UseClusters(track);
407 }
408 Bool_t isGold = kFALSE;
409
410 if (track->GetChi2()/track->GetNumberOfClusters()<5) { //full gold track
7ad19338 411 // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
412 if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
4f1c04d3 413 isGold = kTRUE;
414 }
415 if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track
7ad19338 416 // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
417 if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
f4e9508c 418 isGold = kTRUE;
419 }
4f1c04d3 420 if (!isGold && track->GetBackupTrack()){
421 if (track->GetBackupTrack()->GetNumberOfClusters()>foundMin&&
422 (track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1))<7){
423 seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
424 isGold = kTRUE;
425 }
426 }
7ad19338 427 if (track->StatusForTOF()>0 &&track->fNCross==0 && Float_t(track->fN)/Float_t(track->fNExpected)>0.4){
428 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
429 }
16d9fbba 430 }
c630aafd 431 }
8979685e 432 // Debug part of tracking
433 TTreeSRedirector& cstream = *fDebugStreamer;
434 Int_t eventNr = event->GetEventNumber();
435 if (track->GetBackupTrack()){
436 cstream<<"Tracks"<<
437 "EventNr="<<eventNr<<
438 "ESD.="<<seed<<
439 "trd.="<<track<<
440 "trdback.="<<track->GetBackupTrack()<<
441 "\n";
442 }else{
443 cstream<<"Tracks"<<
444 "EventNr="<<eventNr<<
445 "ESD.="<<seed<<
446 "trd.="<<track<<
447 "trdback.="<<track<<
448 "\n";
449 }
450 //
451 //Propagation to the TOF (I.Belikov)
3c625a9b 452 if (track->GetStop()==kFALSE){
4f1c04d3 453
b94f0a96 454 Double_t xtof=371.;
3c625a9b 455 Double_t c2=track->GetC()*xtof - track->GetEta();
4f1c04d3 456 if (TMath::Abs(c2)>=0.99) {
c5a8e3df 457 delete track;
458 continue;
459 }
4f1c04d3 460 Double_t xTOF0 = 365. ;
16d9fbba 461 PropagateToOuterPlane(*track,xTOF0);
4f1c04d3 462 //
463 //energy losses taken to the account - check one more time
464 c2=track->GetC()*xtof - track->GetEta();
465 if (TMath::Abs(c2)>=0.99) {
466 delete track;
467 continue;
468 }
469
7bed16a7 470 //
3c625a9b 471 Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
472 Double_t y=track->GetYat(xtof);
473 if (y > ymax) {
7ac6fa52 474 if (!track->Rotate(AliTRDgeometry::GetAlpha())) {
475 delete track;
7bed16a7 476 continue;
7ac6fa52 477 }
3c625a9b 478 } else if (y <-ymax) {
7ac6fa52 479 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
480 delete track;
7bed16a7 481 continue;
7ac6fa52 482 }
3c625a9b 483 }
484
485 if (track->PropagateTo(xtof)) {
eab5961e 486 seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
487 for (Int_t i=0;i<kNPlane;i++) {
488 seed->SetTRDsignals(track->GetPIDsignals(i),i);
489 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
490 }
7ad19338 491 // seed->SetTRDtrack(new AliTRDtrack(*track));
3c625a9b 492 if (track->GetNumberOfClusters()>foundMin) found++;
493 }
494 }else{
495 if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
496 seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
16d9fbba 497 //seed->SetStatus(AliESDtrack::kTRDStop);
eab5961e 498 for (Int_t i=0;i<kNPlane;i++) {
499 seed->SetTRDsignals(track->GetPIDsignals(i),i);
500 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
501 }
7ad19338 502 //seed->SetTRDtrack(new AliTRDtrack(*track));
3c625a9b 503 found++;
504 }
1e9bb598 505 }
7ad19338 506 seed->SetTRDQuality(track->StatusForTOF());
8979685e 507 seed->SetTRDBudget(track->fBudget[0]);
508
d9b8978b 509 delete track;
7ad19338 510 //
1e9bb598 511 //End of propagation to the TOF
3c625a9b 512 //if (foundClr>foundMin)
513 // seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
514
c630aafd 515
516 }
517
518 cerr<<"Number of seeds: "<<fNseeds<<endl;
519 cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
69b55c55 520
521 // MakeSeedsMI(3,5,event); //new seeding
7ad19338 522
523
1e9bb598 524 fSeeds->Clear(); fNseeds=0;
4f1c04d3 525 delete [] index;
526 delete [] quality;
527
1e9bb598 528 return 0;
529
530}
531
532//_____________________________________________________________________________
533Int_t AliTRDtracker::RefitInward(AliESD* event)
534{
535 //
536 // Refits tracks within the TRD. The ESD event is expected to contain seeds
537 // at the outer part of the TRD.
538 // The tracks are propagated to the innermost time bin
539 // of the TRD and the ESD event is updated
540 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
541 //
542
543 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
544 Float_t foundMin = fgkMinClustersInTrack * timeBins;
545 Int_t nseed = 0;
546 Int_t found = 0;
547 Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
4f1c04d3 548 AliTRDtrack seed2;
1e9bb598 549
550 Int_t n = event->GetNumberOfTracks();
551 for (Int_t i=0; i<n; i++) {
552 AliESDtrack* seed=event->GetTrack(i);
4f1c04d3 553 new(&seed2) AliTRDtrack(*seed);
554 if (seed2.GetX()<270){
555 seed->UpdateTrackParams(&seed2, AliESDtrack::kTRDbackup); // backup TPC track - only update
f4e9508c 556 continue;
557 }
558
1e9bb598 559 ULong_t status=seed->GetStatus();
0dd7d129 560 if ( (status & AliESDtrack::kTRDout ) == 0 ) {
0dd7d129 561 continue;
562 }
563 if ( (status & AliESDtrack::kTRDin) != 0 ) {
0dd7d129 564 continue;
565 }
f4e9508c 566 nseed++;
7ad19338 567// if (1/seed2.Get1Pt()>1.5&& seed2.GetX()>260.) {
568// Double_t oldx = seed2.GetX();
569// seed2.PropagateTo(500.);
570// seed2.ResetCovariance(1.);
571// seed2.PropagateTo(oldx);
572// }
573// else{
574// seed2.ResetCovariance(5.);
575// }
4f1c04d3 576
577 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
578 UInt_t * indexes2 = seed2.GetIndexes();
7ad19338 579 for (Int_t i=0;i<kNPlane;i++) {
580 pt->SetPIDsignals(seed2.GetPIDsignals(i),i);
581 pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
582 }
eab5961e 583
46e2d86c 584 UInt_t * indexes3 = pt->GetBackupIndexes();
585 for (Int_t i=0;i<200;i++) {
586 if (indexes2[i]==0) break;
587 indexes3[i] = indexes2[i];
588 }
589 //AliTRDtrack *pt = seed2;
1e9bb598 590 AliTRDtrack &t=*pt;
f6625211 591 FollowProlongation(t, innerTB);
1e9bb598 592 if (t.GetNumberOfClusters() >= foundMin) {
46e2d86c 593 // UseClusters(&t);
594 //CookLabel(pt, 1-fgkLabelFraction);
7ad19338 595 t.CookdEdx();
596 CookdEdxTimBin(t);
1e9bb598 597 }
598 found++;
599// cout<<found<<'\r';
600
601 if(PropagateToTPC(t)) {
0fa7dfa7 602 seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
7ad19338 603 for (Int_t i=0;i<kNPlane;i++) {
604 seed->SetTRDsignals(pt->GetPIDsignals(i),i);
605 seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
606 }
7bed16a7 607 }else{
608 //if not prolongation to TPC - propagate without update
609 AliTRDtrack* seed2 = new AliTRDtrack(*seed);
610 seed2->ResetCovariance(5.);
611 AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
612 delete seed2;
eab5961e 613 if (PropagateToTPC(*pt2)) {
3551db50 614 //pt2->CookdEdx(0.,1.);
615 pt2->CookdEdx( ); // Modification by PS
eab5961e 616 CookdEdxTimBin(*pt2);
7bed16a7 617 seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit);
7ad19338 618 for (Int_t i=0;i<kNPlane;i++) {
619 seed->SetTRDsignals(pt2->GetPIDsignals(i),i);
620 seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
621 }
eab5961e 622 }
7bed16a7 623 delete pt2;
1e9bb598 624 }
1e9bb598 625 delete pt;
eab5961e 626 }
1e9bb598 627
628 cout<<"Number of loaded seeds: "<<nseed<<endl;
629 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
630
c630aafd 631 return 0;
632
633}
634
bbf92647 635
8979685e 636
637
638//---------------------------------------------------------------------------
f6625211 639Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
8979685e 640{
641 // Starting from current position on track=t this function tries
642 // to extrapolate the track up to timeBin=0 and to confirm prolongation
643 // if a close cluster is found. Returns the number of clusters
644 // expected to be found in sensitive layers
645 // GeoManager used to estimate mean density
646 Int_t sector;
647 Int_t lastplane = GetLastPlane(&t);
8979685e 648 Double_t alpha=t.GetAlpha();
649 alpha = TVector2::Phi_0_2pi(alpha);
3551db50 650 Double_t radLength = 0.0;
651 Double_t rho = 0.0;
652 Double_t x, dx;
8979685e 653 //, y, ymax, z;
654 Int_t expectedNumberOfClusters = 0;
655 Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
3551db50 656 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
8979685e 657 //
658 //
659 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
660 Double_t tanmax = TMath::Tan(0.5*alpha);
661
662 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
663 //
664 // propagate track in non active layers
665 //
666 if (!(fTrSec[0]->GetLayer(nr)->IsSensitive())){
667 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
668 t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
669 while (nr >rf && (!(fTrSec[0]->GetLayer(nr)->IsSensitive()))){
670 x = fTrSec[0]->GetLayer(nr)->GetX();
671 nr--;
672 if (!t.GetProlongation(x,y,z)) break;
673 if (TMath::Abs(y)>x*tanmax){
674 nr--;
675 break;
676 }
677 }
678 nr++;
679 x = fTrSec[0]->GetLayer(nr)->GetX();
680 if (!t.GetProlongation(x,y,z)) break;
681 //
682 // minimal mean and maximal budget scan
683 Float_t minbudget =10000;
684 Float_t meanbudget =0;
685 Float_t maxbudget =-1;
686// Float_t normbudget =0;
687// for (Int_t idy=-1;idy<=1;idy++)
688// for (Int_t idz=-1;idz<=1;idz++){
689 for (Int_t idy=0;idy<1;idy++)
690 for (Int_t idz=0;idz<1;idz++){
691 Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
692 Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
693 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
694 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
695 xyz1[2] = z2;
696 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
697 Float_t budget = param[0]*param[4];
698 meanbudget+=budget;
699 if (budget<minbudget) minbudget=budget;
700 if (budget>maxbudget) maxbudget=budget;
701 }
702 t.fBudget[0]+=minbudget;
703 t.fBudget[1]+=meanbudget/9.;
704 t.fBudget[2]+=minbudget;
705 //
706 //
707 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
708 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
709 xyz1[2] = z;
710 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
711
712 t.PropagateTo(x,param[1],param[0]);
713 t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //end global position
714 AdjustSector(&t);
715 continue;
716 }
717 //
718 //
719 // stop tracking for highly inclined tracks
720 if (!AdjustSector(&t)) break;
721 if (TMath::Abs(t.GetSnp())>0.95) break;
722 //
723 // propagate and update track in active layers
724 //
725 Int_t nr0 = nr; //first active layer
726 if (nr >rf && (fTrSec[0]->GetLayer(nr)->IsSensitive())){
727 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
728 t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
729 while (nr >rf && ((fTrSec[0]->GetLayer(nr)->IsSensitive()))){
730 x = fTrSec[0]->GetLayer(nr)->GetX();
731 nr--;
732 if (!t.GetProlongation(x,y,z)) break;
733 if (TMath::Abs(y)>x*tanmax){
734 nr--;
735 break;
736 }
737 }
738 // nr++;
739 x = fTrSec[0]->GetLayer(nr)->GetX();
740 if (!t.GetProlongation(x,y,z)) break;
741 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
742 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
743 xyz1[2] = z;
744 // end global position
745 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
746 rho = param[0];
747 radLength = param[1]; // get mean propagation parameters
748 }
749 //
750 // propagate and update
751 if (nr0-nr<5){
752 // short tracklet - do not update - edge effect
753 x = fTrSec[0]->GetLayer(nr)->GetX();
754 t.PropagateTo(x,radLength,rho);
755 AdjustSector(&t);
756 continue;
757 }
758 sector = t.GetSector();
759 //
760 //
761 for (Int_t ilayer=nr0;ilayer>=nr;ilayer--) {
762 expectedNumberOfClusters++;
763 t.fNExpected++;
764 if (t.fX>345) t.fNExpectedLast++;
765 AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
766 AliTRDcluster *cl=0;
767 UInt_t index=0;
768 Double_t maxChi2=fgkMaxChi2;
769 dx = (fTrSec[sector]->GetLayer(ilayer+1))->GetX()-timeBin.GetX();
770 x = timeBin.GetX();
771 t.PropagateTo(x,radLength,rho);
772 // Now go for the real cluster search
773 if (timeBin) {
774 AliTRDcluster * cl0 = timeBin[0];
775 if (!cl0) continue; // no clusters in given time bin
776 Int_t plane = fGeom->GetPlane(cl0->GetDetector());
777 if (plane>lastplane) continue;
778 Int_t timebin = cl0->GetLocalTimeBin();
779 AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index);
780 //
781 if (cl2) {
782 cl =cl2;
783 Double_t h01 = GetTiltFactor(cl);
784 maxChi2=t.GetPredictedChi2(cl,h01);
785 }
786
787 if (cl) {
788 // if (cl->GetNPads()<5)
789 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
790 Double_t h01 = GetTiltFactor(cl);
791 Int_t det = cl->GetDetector();
792 Int_t plane = fGeom->GetPlane(det);
793 if (t.fX>345){
794 t.fNLast++;
795 t.fChi2Last+=maxChi2;
796 }
797 if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
798 if(!t.Update(cl,maxChi2,index,h01)) {
8979685e 799 }
800 }
8979685e 801 }
802 }
803 }
804 }
805 return expectedNumberOfClusters;
806
807
808}
809
8979685e 810
811
69b55c55 812
8979685e 813//___________________________________________________________________
f6625211 814Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
8979685e 815{
7ad19338 816
8979685e 817 // Starting from current radial position of track <t> this function
818 // extrapolates the track up to outer timebin and in the sensitive
819 // layers confirms prolongation if a close cluster is found.
820 // Returns the number of clusters expected to be found in sensitive layers
821 // Use GEO manager for material Description
8979685e 822 Double_t alpha=t.GetAlpha();
823 TVector2::Phi_0_2pi(alpha);
824 Int_t sector;
825 Int_t clusters[1000];
826 for (Int_t i=0;i<1000;i++) clusters[i]=-1;
827 Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
3551db50 828 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
8979685e 829 Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
3551db50 830 Double_t radLength = 0.0;
831 Double_t rho = 0.0;
832 Double_t x, dx; //y, z;
8979685e 833
834 Int_t expectedNumberOfClusters = 0;
835 x = t.GetX();
836
837 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
838 Double_t tanmax = TMath::Tan(0.5*alpha);
839 Int_t nr;
840 Float_t ratio0=0;
841 AliTRDtracklet tracklet;
842 //
843 //
844 //
845 for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) {
846 //
847 // propagate track in non active layers
848 //
849 if (!(fTrSec[0]->GetLayer(nr)->IsSensitive())){
850 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
851 t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
852 while (nr <outerTB && (!(fTrSec[0]->GetLayer(nr)->IsSensitive()))){
853 x = fTrSec[0]->GetLayer(nr)->GetX();
854 nr++;
855 if (!t.GetProlongation(x,y,z)) break;
856 if (TMath::Abs(y)>x*tanmax){
857 nr++;
858 break;
859 }
860 }
861 nr--;
862 x = fTrSec[0]->GetLayer(nr)->GetX();
863 if (!t.GetProlongation(x,y,z)) break;
864 // minimal mean and maximal budget scan
865 Float_t minbudget =10000;
866 Float_t meanbudget =0;
867 Float_t maxbudget =-1;
868// Float_t normbudget =0;
869// for (Int_t idy=-1;idy<=1;idy++)
870// for (Int_t idz=-1;idz<=1;idz++){
871 for (Int_t idy=0;idy<1;idy++)
872 for (Int_t idz=0;idz<1;idz++){
873 Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
3a946bfb 874 Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaZ2()),1.);
8979685e 875
876 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
877 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
878 xyz1[2] = z2;
879 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
880 Float_t budget = param[0]*param[4];
881 meanbudget+=budget;
882 if (budget<minbudget) minbudget=budget;
883 if (budget>maxbudget) maxbudget=budget;
884 }
885 t.fBudget[0]+=minbudget;
886 t.fBudget[1]+=meanbudget/9.;
887 t.fBudget[2]+=minbudget;
888
889 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
890 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
891 xyz1[2] = z;
892 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
893 t.PropagateTo(x,param[1],param[0]);
894 t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //end global position
895 AdjustSector(&t);
896 continue;
897 }
898 //
899 //
900 // stop tracking for highly inclined tracks
901 if (!AdjustSector(&t)) break;
902 if (TMath::Abs(t.GetSnp())>0.95) break;
903 //
904 // propagate and update track in active layers
905 //
906 Int_t nr0 = nr; //first active layer
907 if (nr <outerTB && (fTrSec[0]->GetLayer(nr)->IsSensitive())){
908 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
909 t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
910 while (nr <outerTB && ((fTrSec[0]->GetLayer(nr)->IsSensitive()))){
911 x = fTrSec[0]->GetLayer(nr)->GetX();
912 nr++;
913 if (!t.GetProlongation(x,y,z)) break;
914 if (TMath::Abs(y)>(x*tanmax)){
915 nr++;
916 break;
917 }
918 }
919 x = fTrSec[0]->GetLayer(nr)->GetX();
920 if (!t.GetProlongation(x,y,z)) break;
921 // minimal mean and maximal budget scan
922 Float_t minbudget =10000;
923 Float_t meanbudget =0;
924 Float_t maxbudget =-1;
925 // Float_t normbudget =0;
926 // for (Int_t idy=-1;idy<=1;idy++)
927 // for (Int_t idz=-1;idz<=1;idz++){
928 for (Int_t idy=0;idy<1;idy++)
929 for (Int_t idz=0;idz<1;idz++){
930 Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
3a946bfb 931 Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaZ2()),1.);
8979685e 932
933 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
934 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
935 xyz1[2] = z2;
936 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
937 Float_t budget = param[0]*param[4];
938 meanbudget+=budget;
939 if (budget<minbudget) minbudget=budget;
940 if (budget>maxbudget) maxbudget=budget;
941 }
942 t.fBudget[0]+=minbudget;
943 t.fBudget[1]+=meanbudget/9.;
944 t.fBudget[2]+=minbudget;
945 //
946 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
947 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
948 xyz1[2] = z;
949 // end global position
950 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
951 rho = param[0];
952 radLength = param[1]; // get mean propagation parameters
953 }
954 //
955 //
956 if (nr-nr0<5){
957 // short tracklet - do not update - edge effect
958 x = fTrSec[0]->GetLayer(nr+1)->GetX();
959 t.PropagateTo(x,radLength,rho);
960 AdjustSector(&t);
961 continue;
962 }
963 //
964 //
965 sector = t.GetSector();
966 Float_t ncl = FindClusters(sector,nr0,nr,&t,clusters,tracklet);
967 if (tracklet.GetN()-2*tracklet.GetNCross()<10) continue;
968 ratio0 = ncl/Float_t(fTimeBinsPerPlane);
969 Float_t ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.);
970 if (tracklet.GetChi2()<18.&&ratio0>0.8 && ratio1>0.6 && ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85&&t.fN>20){
971 t.MakeBackupTrack(); // make backup of the track until is gold
972 }
973 //
974 //
975 for (Int_t ilayer=nr0;ilayer<=nr;ilayer++) {
976 expectedNumberOfClusters++;
977 t.fNExpected++;
978 if (t.fX>345) t.fNExpectedLast++;
979 AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
980 AliTRDcluster *cl=0;
981 UInt_t index=0;
982 Double_t maxChi2=fgkMaxChi2;
983 dx = (fTrSec[sector]->GetLayer(ilayer-1))->GetX()-timeBin.GetX();
984 x = timeBin.GetX();
985 t.PropagateTo(x,radLength,rho);
986 // Now go for the real cluster search
987 if (timeBin) {
988 if (clusters[ilayer]>0) {
989 index = clusters[ilayer];
990 cl = (AliTRDcluster*)GetCluster(index);
991 Double_t h01 = GetTiltFactor(cl);
992 maxChi2=t.GetPredictedChi2(cl,h01);
993 }
994
995 if (cl) {
996 // if (cl->GetNPads()<5)
997 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
998 Double_t h01 = GetTiltFactor(cl);
999 Int_t det = cl->GetDetector();
1000 Int_t plane = fGeom->GetPlane(det);
1001 if (t.fX>345){
1002 t.fNLast++;
1003 t.fChi2Last+=maxChi2;
1004 }
1005 if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1006 if(!t.Update(cl,maxChi2,index,h01)) {
8979685e 1007 }
1008 }
8979685e 1009 //
1010
1011 if (cl->GetLocalTimeBin()==1&&t.fN>20 && float(t.fChi2)/float(t.fN)<5){
1012 Float_t ratio1 = Float_t(t.fN)/Float_t(t.fNExpected);
1013 if (tracklet.GetChi2()<18&&ratio0>0.8&&ratio1>0.6 &&ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85){
1014 t.MakeBackupTrack(); // make backup of the track until is gold
1015 }
1016 }
1017 // reset material budget if 2 consecutive gold
1018 if (plane>0)
1019 if (t.fTracklets[plane].GetN()+t.fTracklets[plane-1].GetN()>20){
1020 t.fBudget[2] = 0;
1021 }
1022 }
1023 }
1024 }
1025 }
1026 //
1027 if (nr<outerTB)
1028 t.SetStop(kTRUE);
1029 else
1030 t.SetStop(kFALSE);
1031 return expectedNumberOfClusters;
5443e65e 1032}
1033
1e9bb598 1034//---------------------------------------------------------------------------
1035Int_t AliTRDtracker::Refit(AliTRDtrack& t, Int_t rf)
1036{
1037 // Starting from current position on track=t this function tries
1038 // to extrapolate the track up to timeBin=0 and to reuse already
1039 // assigned clusters. Returns the number of clusters
1040 // expected to be found in sensitive layers
1041 // get indices of assigned clusters for each layer
1042 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
8979685e 1043 Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
3551db50 1044 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
1e9bb598 1045 Int_t iCluster[90];
1046 for (Int_t i = 0; i < 90; i++) iCluster[i] = 0;
1047 for (Int_t i = 0; i < t.GetNumberOfClusters(); i++) {
1048 Int_t index = t.GetClusterIndex(i);
1049 AliTRDcluster *cl=(AliTRDcluster*) GetCluster(index);
1050 if (!cl) continue;
1051 Int_t detector=cl->GetDetector();
1052 Int_t localTimeBin=cl->GetLocalTimeBin();
1053 Int_t sector=fGeom->GetSector(detector);
1054 Int_t plane=fGeom->GetPlane(detector);
1055
1056 Int_t trackingSector = CookSectorIndex(sector);
1057
1058 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1059 if(gtb < 0) continue;
1060 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1061 iCluster[layer] = index;
1062 }
1063 t.ResetClusters();
1064
1065 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1066
1067 Double_t alpha=t.GetAlpha();
1068 alpha = TVector2::Phi_0_2pi(alpha);
1069
1070 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1071 Double_t radLength, rho, x, dx, y, ymax, z;
1072
1073 Int_t expectedNumberOfClusters = 0;
1074 Bool_t lookForCluster;
1075
1076 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1077
1078
1079 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
1080
1081 y = t.GetY(); z = t.GetZ();
1082
1083 // first propagate to the inner surface of the current time bin
1084 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1085 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1086 if(!t.PropagateTo(x,radLength,rho)) break;
1087 y = t.GetY();
1088 ymax = x*TMath::Tan(0.5*alpha);
1089 if (y > ymax) {
1090 s = (s+1) % ns;
1091 if (!t.Rotate(alpha)) break;
1092 if(!t.PropagateTo(x,radLength,rho)) break;
1093 } else if (y <-ymax) {
1094 s = (s-1+ns) % ns;
1095 if (!t.Rotate(-alpha)) break;
1096 if(!t.PropagateTo(x,radLength,rho)) break;
1097 }
1098
1099 y = t.GetY(); z = t.GetZ();
1100
1101 // now propagate to the middle plane of the next time bin
1102 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1103 x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1104 if(!t.PropagateTo(x,radLength,rho)) break;
1105 y = t.GetY();
1106 ymax = x*TMath::Tan(0.5*alpha);
1107 if (y > ymax) {
1108 s = (s+1) % ns;
1109 if (!t.Rotate(alpha)) break;
1110 if(!t.PropagateTo(x,radLength,rho)) break;
1111 } else if (y <-ymax) {
1112 s = (s-1+ns) % ns;
1113 if (!t.Rotate(-alpha)) break;
1114 if(!t.PropagateTo(x,radLength,rho)) break;
1115 }
1116
1117 if(lookForCluster) expectedNumberOfClusters++;
1118
1119 // use assigned cluster
1120 if (!iCluster[nr-1]) continue;
1121 AliTRDcluster *cl=(AliTRDcluster*)GetCluster(iCluster[nr-1]);
1122 Double_t h01 = GetTiltFactor(cl);
1123 Double_t chi2=t.GetPredictedChi2(cl, h01);
8979685e 1124 //if (cl->GetNPads()<5)
1125 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
4f1c04d3 1126
1127 //t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1e9bb598 1128 t.Update(cl,chi2,iCluster[nr-1],h01);
1129 }
1130
1131 return expectedNumberOfClusters;
1132}
1133
5443e65e 1134//___________________________________________________________________
1135
1136Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1137{
1138 // Starting from current radial position of track <t> this function
1139 // extrapolates the track up to radial position <xToGo>.
1140 // Returns 1 if track reaches the plane, and 0 otherwise
1141
1142 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1143
1144 Double_t alpha=t.GetAlpha();
1145
1146 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1147 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1148
1149 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1150
1151 Bool_t lookForCluster;
029cd327 1152 Double_t radLength, rho, x, dx, y, ymax, z;
5443e65e 1153
1154 x = t.GetX();
1155
1156 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1157
1158 Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1159
1160 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
1161
1162 y = t.GetY(); z = t.GetZ();
1163
1164 // first propagate to the outer surface of the current time bin
029cd327 1165 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 1166 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
029cd327 1167 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1168 y = t.GetY();
1169 ymax = x*TMath::Tan(0.5*alpha);
1170 if (y > ymax) {
1171 s = (s+1) % ns;
1172 if (!t.Rotate(alpha)) return 0;
1173 } else if (y <-ymax) {
1174 s = (s-1+ns) % ns;
1175 if (!t.Rotate(-alpha)) return 0;
1176 }
029cd327 1177 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1178
1179 y = t.GetY(); z = t.GetZ();
1180
1181 // now propagate to the middle plane of the next time bin
029cd327 1182 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
5443e65e 1183 x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
029cd327 1184 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1185 y = t.GetY();
1186 ymax = x*TMath::Tan(0.5*alpha);
1187 if (y > ymax) {
1188 s = (s+1) % ns;
1189 if (!t.Rotate(alpha)) return 0;
1190 } else if (y <-ymax) {
1191 s = (s-1+ns) % ns;
1192 if (!t.Rotate(-alpha)) return 0;
1193 }
029cd327 1194 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1195 }
1196 return 1;
1197}
1198
1199//___________________________________________________________________
1200
1201Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1202{
1203 // Starting from current radial position of track <t> this function
1204 // extrapolates the track up to radial position of the outermost
1205 // padrow of the TPC.
1206 // Returns 1 if track reaches the TPC, and 0 otherwise
1207
c630aafd 1208 //Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
5443e65e 1209
1210 Double_t alpha=t.GetAlpha();
c630aafd 1211 alpha = TVector2::Phi_0_2pi(alpha);
5443e65e 1212
1213 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1214
1215 Bool_t lookForCluster;
029cd327 1216 Double_t radLength, rho, x, dx, y, /*ymax,*/ z;
5443e65e 1217
1218 x = t.GetX();
1219
1220 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
5443e65e 1221 Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1222
1223 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1224
9c9d2487 1225 y = t.GetY();
1226 z = t.GetZ();
5443e65e 1227
1228 // first propagate to the outer surface of the current time bin
029cd327 1229 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
9c9d2487 1230 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2;
1231
029cd327 1232 if(!t.PropagateTo(x,radLength,rho)) return 0;
9c9d2487 1233 AdjustSector(&t);
029cd327 1234 if(!t.PropagateTo(x,radLength,rho)) return 0;
5443e65e 1235
9c9d2487 1236 y = t.GetY();
1237 z = t.GetZ();
5443e65e 1238
1239 // now propagate to the middle plane of the next time bin
029cd327 1240 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
9c9d2487 1241 x = fTrSec[s]->GetLayer(nr-1)->GetX();
1242
029cd327 1243 if(!t.PropagateTo(x,radLength,rho)) return 0;
9c9d2487 1244 AdjustSector(&t);
029cd327 1245 if(!t.PropagateTo(x,radLength,rho)) return 0;
c5a8e3df 1246 }
5443e65e 1247 return 1;
1248}
1249
c630aafd 1250//_____________________________________________________________________________
1251Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1252{
1253 // Fills clusters into TRD tracking_sectors
1254 // Note that the numbering scheme for the TRD tracking_sectors
1255 // differs from that of TRD sectors
4f1c04d3 1256 cout<<"\n Read Sectors clusters"<<endl;
c630aafd 1257 if (ReadClusters(fClusters,cTree)) {
1258 Error("LoadClusters","Problem with reading the clusters !");
1259 return 1;
1260 }
1261 Int_t ncl=fClusters->GetEntriesFast();
b7a0917f 1262 fNclusters=ncl;
c630aafd 1263 cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1264
1265 UInt_t index;
3c625a9b 1266 for (Int_t ichamber=0;ichamber<5;ichamber++)
1267 for (Int_t isector=0;isector<18;isector++){
1268 fHoles[ichamber][isector]=kTRUE;
1269 }
1270
1271
c630aafd 1272 while (ncl--) {
1273// printf("\r %d left ",ncl);
1274 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
029cd327 1275 Int_t detector=c->GetDetector();
1276 Int_t localTimeBin=c->GetLocalTimeBin();
c630aafd 1277 Int_t sector=fGeom->GetSector(detector);
1278 Int_t plane=fGeom->GetPlane(detector);
3c625a9b 1279
029cd327 1280 Int_t trackingSector = CookSectorIndex(sector);
3c625a9b 1281 if (c->GetLabel(0)>0){
1282 Int_t chamber = fGeom->GetChamber(detector);
1283 fHoles[chamber][trackingSector]=kFALSE;
1284 }
c630aafd 1285
029cd327 1286 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
c630aafd 1287 if(gtb < 0) continue;
029cd327 1288 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1289
1290 index=ncl;
69b55c55 1291 //
1292 // apply pos correction
1293 Float_t poscor = fgkCoef*(c->GetLocalTimeBin() - fgkMean)+fgkOffset;
f6625211 1294 //Float_t poscor=0;
69b55c55 1295 c->SetY(c->GetY()-poscor);
029cd327 1296 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
c630aafd 1297 }
7bed16a7 1298 // printf("\r\n");
3c625a9b 1299 //
1300 //
1301 /*
1302 for (Int_t isector=0;isector<18;isector++){
1303 for (Int_t ichamber=0;ichamber<5;ichamber++)
1304 if (fHoles[ichamber][isector]!=fGeom->IsHole(0,ichamber,17-isector))
1305 printf("Problem \t%d\t%d\t%d\t%d\n",isector,ichamber,fHoles[ichamber][isector],
1306 fGeom->IsHole(0,ichamber,17-isector));
1307 }
1308 */
c630aafd 1309 return 0;
1310}
1311
5443e65e 1312//_____________________________________________________________________________
b7a0917f 1313void AliTRDtracker::UnloadClusters()
5443e65e 1314{
1315 //
1316 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1317 //
1318
1319 Int_t i, nentr;
1320
1321 nentr = fClusters->GetEntriesFast();
1322 for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
b7a0917f 1323 fNclusters = 0;
5443e65e 1324
1325 nentr = fSeeds->GetEntriesFast();
1326 for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1327
1328 nentr = fTracks->GetEntriesFast();
1329 for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1330
1331 Int_t nsec = AliTRDgeometry::kNsect;
1332
1333 for (i = 0; i < nsec; i++) {
1334 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1335 fTrSec[i]->GetLayer(pl)->Clear();
1336 }
1337 }
1338
1339}
1340
7ad19338 1341//__________________________________________________________________________
69b55c55 1342void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD * esd)
7ad19338 1343{
1344 //
1345 // Creates seeds using clusters between position inner plane and outer plane
1346 //
69b55c55 1347 const Double_t maxtheta = 1;
1348 const Double_t maxphi = 2.0;
1349 //
1350 const Double_t kRoad0y = 6; // road for middle cluster
1351 const Double_t kRoad0z = 8.5; // road for middle cluster
1352 //
1353 const Double_t kRoad1y = 2; // road in y for seeded cluster
1354 const Double_t kRoad1z = 20; // road in z for seeded cluster
1355 //
1356 const Double_t kRoad2y = 3; // road in y for extrapolated cluster
1357 const Double_t kRoad2z = 20; // road in z for extrapolated cluster
1358 const Int_t maxseed = 3000;
1359 Int_t maxSec=AliTRDgeometry::kNsect;
7ad19338 1360
69b55c55 1361 //
1362 // linear fitters in planes
1363 TLinearFitter fitterTC(2,"hyp2"); // fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1364 TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
1365 fitterTC.StoreData(kTRUE);
1366 fitterT2.StoreData(kTRUE);
1367 AliRieman rieman(1000); // rieman fitter
1368 AliRieman rieman2(1000); // rieman fitter
7ad19338 1369 //
1370 // find the maximal and minimal layer for the planes
7ad19338 1371 //
1372 Int_t layers[6][2];
69b55c55 1373 AliTRDpropagationLayer* reflayers[6];
7ad19338 1374 for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;}
7ad19338 1375 for (Int_t ns=0;ns<maxSec;ns++){
1376 for (Int_t ilayer=0;ilayer<fTrSec[ns]->GetNumberOfLayers();ilayer++){
1377 AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer));
1378 if (layer==0) continue;
1379 Int_t det = layer[0]->GetDetector();
1380 Int_t plane = fGeom->GetPlane(det);
1381 if (ilayer<layers[plane][0]) layers[plane][0] = ilayer;
1382 if (ilayer>layers[plane][1]) layers[plane][1] = ilayer;
1383 }
1384 }
1385 //
3551db50 1386 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
69b55c55 1387 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
1388 Double_t hL[6]; // tilting angle
1389 Double_t xcl[6]; // x - position of reference cluster
1390 Double_t ycl[6]; // y - position of reference cluster
1391 Double_t zcl[6]; // z - position of reference cluster
1392 AliTRDcluster *cl[6]={0,0,0,0,0,0}; // seeding clusters
1393 Float_t padlength[6]={10,10,10,10,10,10}; //current pad-length
1394 Double_t chi2R =0, chi2Z=0;
1395 Double_t chi2RF =0, chi2ZF=0;
1396 //
1397 Int_t nclusters; // total number of clusters
1398 for (Int_t i=0;i<6;i++) {hL[i]=h01; if (i%2==1) hL[i]*=-1.;}
1399 //
1400 //
1401 // registered seed
1402 AliTRDseed *pseed = new AliTRDseed[maxseed*6];
1403 AliTRDseed *seed[maxseed];
1404 for (Int_t iseed=0;iseed<maxseed;iseed++) seed[iseed]= &pseed[iseed*6];
1405 AliTRDseed *cseed = seed[0];
1406 //
1407 Double_t seedquality[maxseed];
1408 Double_t seedquality2[maxseed];
1409 Double_t seedparams[maxseed][7];
1410 Int_t seedlayer[maxseed];
1411 Int_t registered =0;
1412 Int_t sort[maxseed];
1413 //
1414 // seeding part
1415 //
1416 for (Int_t ns = 0; ns<maxSec; ns++){ //loop over sectors
1417 //for (Int_t ns = 0; ns<5; ns++){ //loop over sectors
1418 registered = 0; // reset registerd seed counter
1419 cseed = seed[registered];
1420 Float_t iter=0;
1421 for (Int_t sLayer=2; sLayer>=0;sLayer--){
1422 //for (Int_t dseed=5;dseed<15; dseed+=3){ //loop over central seeding time bins
1423 iter+=1.;
1424 Int_t dseed = 5+Int_t(iter)*3;
1425 // Initialize seeding layers
1426 for (Int_t ilayer=0;ilayer<6;ilayer++){
1427 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1428 xcl[ilayer] = reflayers[ilayer]->GetX();
1429 }
7ad19338 1430 //
69b55c55 1431 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2])*0.5;
1432 AliTRDpropagationLayer& layer0=*reflayers[sLayer+0];
1433 AliTRDpropagationLayer& layer1=*reflayers[sLayer+1];
1434 AliTRDpropagationLayer& layer2=*reflayers[sLayer+2];
1435 AliTRDpropagationLayer& layer3=*reflayers[sLayer+3];
1436 //
1437 Int_t maxn3 = layer3;
1438 for (Int_t icl3=0;icl3<maxn3;icl3++){
1439 AliTRDcluster *cl3 = layer3[icl3];
1440 if (!cl3) continue;
1441 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2()*12.);
1442 ycl[sLayer+3] = cl3->GetY();
1443 zcl[sLayer+3] = cl3->GetZ();
1444 Float_t yymin0 = ycl[sLayer+3] - 1- maxphi *(xcl[sLayer+3]-xcl[sLayer+0]);
1445 Float_t yymax0 = ycl[sLayer+3] + 1+ maxphi *(xcl[sLayer+3]-xcl[sLayer+0]);
1446 Int_t maxn0 = layer0; //
1447 for (Int_t icl0=layer0.Find(yymin0);icl0<maxn0;icl0++){
1448 AliTRDcluster *cl0 = layer0[icl0];
1449 if (!cl0) continue;
1450 if (cl3->IsUsed()&&cl0->IsUsed()) continue;
1451 ycl[sLayer+0] = cl0->GetY();
1452 zcl[sLayer+0] = cl0->GetZ();
1453 if ( ycl[sLayer+0]>yymax0) break;
1454 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
1455 if (TMath::Abs(tanphi)>maxphi) continue;
1456 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
1457 if (TMath::Abs(tantheta)>maxtheta) continue;
1458 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2()*12.);
1459 //
1460 // expected position in 1 layer
1461 Double_t y1exp = ycl[sLayer+0]+(tanphi) *(xcl[sLayer+1]-xcl[sLayer+0]);
1462 Double_t z1exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+1]-xcl[sLayer+0]);
1463 Float_t yymin1 = y1exp - kRoad0y-tanphi;
1464 Float_t yymax1 = y1exp + kRoad0y+tanphi;
1465 Int_t maxn1 = layer1; //
1466 //
1467 for (Int_t icl1=layer1.Find(yymin1);icl1<maxn1;icl1++){
1468 AliTRDcluster *cl1 = layer1[icl1];
1469 if (!cl1) continue;
1470 Int_t nusedCl = 0;
1471 if (cl3->IsUsed()) nusedCl++;
1472 if (cl0->IsUsed()) nusedCl++;
1473 if (cl1->IsUsed()) nusedCl++;
1474 if (nusedCl>1) continue;
1475 ycl[sLayer+1] = cl1->GetY();
1476 zcl[sLayer+1] = cl1->GetZ();
1477 if ( ycl[sLayer+1]>yymax1) break;
1478 if (TMath::Abs(ycl[sLayer+1]-y1exp)>kRoad0y+tanphi) continue;
1479 if (TMath::Abs(zcl[sLayer+1]-z1exp)>kRoad0z) continue;
1480 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2()*12.);
1481 //
1482 Double_t y2exp = ycl[sLayer+0]+(tanphi) *(xcl[sLayer+2]-xcl[sLayer+0])+(ycl[sLayer+1]-y1exp);
1483 Double_t z2exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+2]-xcl[sLayer+0]);
1484 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y, kRoad1z);
1485 if (index2<=0) continue;
1486 AliTRDcluster *cl2 = (AliTRDcluster*)GetCluster(index2);
1487 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2()*12.);
1488 ycl[sLayer+2] = cl2->GetY();
1489 zcl[sLayer+2] = cl2->GetZ();
1490 if (TMath::Abs(cl2->GetZ()-z2exp)>kRoad0z) continue;
1491 //
1492 rieman.Reset();
1493 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1494 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1495 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1496 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1497 rieman.Update();
1498 //
1499 // reset fitter
1500 for (Int_t iLayer=0;iLayer<6;iLayer++){
1501 cseed[iLayer].Reset();
1502 }
1503 chi2Z =0.; chi2R=0.;
1504 for (Int_t iLayer=0;iLayer<4;iLayer++){
1505 cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
1506 chi2Z += (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer])*
1507 (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer]);
1508 cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1509 cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
1510 chi2R += (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer])*
1511 (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer]);
1512 cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1513 }
1514 if (TMath::Sqrt(chi2R)>1./iter) continue;
1515 if (TMath::Sqrt(chi2Z)>7./iter) continue;
1516 //
1517 //
1518 //
1519 Float_t minmax[2]={-100,100};
1520 for (Int_t iLayer=0;iLayer<4;iLayer++){
1521 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer]*0.5+1 -cseed[sLayer+iLayer].fZref[0];
1522 if (max<minmax[1]) minmax[1]=max;
1523 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer]*0.5-1 -cseed[sLayer+iLayer].fZref[0];
1524 if (min>minmax[0]) minmax[0]=min;
1525 }
1526 Bool_t isFake = kFALSE;
1527 if (cl0->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1528 if (cl1->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1529 if (cl2->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1530 if ((!isFake) || (icl3%10)==0 ){ //debugging print
1531 TTreeSRedirector& cstream = *fDebugStreamer;
1532 cstream<<"Seeds0"<<
1533 "isFake="<<isFake<<
1534 "Cl0.="<<cl0<<
1535 "Cl1.="<<cl1<<
1536 "Cl2.="<<cl2<<
1537 "Cl3.="<<cl3<<
1538 "Xref="<<xref<<
1539 "X0="<<xcl[sLayer+0]<<
1540 "X1="<<xcl[sLayer+1]<<
1541 "X2="<<xcl[sLayer+2]<<
1542 "X3="<<xcl[sLayer+3]<<
1543 "Y2exp="<<y2exp<<
1544 "Z2exp="<<z2exp<<
1545 "Chi2R="<<chi2R<<
1546 "Chi2Z="<<chi2Z<<
1547 "Seed0.="<<&cseed[sLayer+0]<<
1548 "Seed1.="<<&cseed[sLayer+1]<<
1549 "Seed2.="<<&cseed[sLayer+2]<<
1550 "Seed3.="<<&cseed[sLayer+3]<<
1551 "Zmin="<<minmax[0]<<
1552 "Zmax="<<minmax[1]<<
1553 "\n";
1554 }
1555
1556 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1557 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1558 //<<<<<<<<<<<<<<<<<< FIT SEEDING PART <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1559 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1560 cl[sLayer+0] = cl0;
1561 cl[sLayer+1] = cl1;
1562 cl[sLayer+2] = cl2;
1563 cl[sLayer+3] = cl3;
1564 Bool_t isOK=kTRUE;
1565 for (Int_t jLayer=0;jLayer<4;jLayer++){
1566 cseed[sLayer+jLayer].fTilt = hL[sLayer+jLayer];
1567 cseed[sLayer+jLayer].fPadLength = padlength[sLayer+jLayer];
1568 cseed[sLayer+jLayer].fX0 = xcl[sLayer+jLayer];
1569 for (Int_t iter=0; iter<2; iter++){
1570 //
1571 // in iteration 0 we try only one pad-row
1572 // if quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
1573 //
1574 AliTRDseed tseed = cseed[sLayer+jLayer];
1575 Float_t roadz = padlength[sLayer+jLayer]*0.5;
1576 if (iter>0) roadz = padlength[sLayer+jLayer];
1577 //
1578 Float_t quality =10000;
1579 for (Int_t iTime=2;iTime<20;iTime++){
1580 AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1581 Double_t dxlayer= layer.GetX()-xcl[sLayer+jLayer];
1582 Double_t zexp = cl[sLayer+jLayer]->GetZ() ;
1583 if (iter>0){
1584 // try 2 pad-rows in second iteration
1585 zexp = tseed.fZref[0]+ tseed.fZref[1]*dxlayer;
1586 if (zexp>cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()+padlength[sLayer+jLayer]*0.5;
1587 if (zexp<cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()-padlength[sLayer+jLayer]*0.5;
1588 }
1589 //
1590 Double_t yexp = tseed.fYref[0]+
1591 tseed.fYref[1]*dxlayer;
1592 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
1593 if (index<=0) continue;
1594 AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1595 //
1596 tseed.fIndexes[iTime] = index;
1597 tseed.fClusters[iTime] = cl; // register cluster
1598 tseed.fX[iTime] = dxlayer; // register cluster
1599 tseed.fY[iTime] = cl->GetY(); // register cluster
1600 tseed.fZ[iTime] = cl->GetZ(); // register cluster
1601 }
1602 tseed.Update();
1603 //count the number of clusters and distortions into quality
1604 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1605 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1606 TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
1607 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1608 if (iter==0 && tseed.isOK()) {
1609 cseed[sLayer+jLayer] = tseed;
1610 quality = tquality;
1611 if (tquality<5) break;
1612 }
1613 if (tseed.isOK() && tquality<quality)
1614 cseed[sLayer+jLayer] = tseed;
1615 }
1616 if (!cseed[sLayer+jLayer].isOK()){
1617 isOK = kFALSE;
1618 break;
1619 }
1620 cseed[sLayer+jLayer].CookLabels();
1621 cseed[sLayer+jLayer].UpdateUsed();
1622 nusedCl+= cseed[sLayer+jLayer].fNUsed;
1623 if (nusedCl>25){
1624 isOK = kFALSE;
1625 break;
1626 }
1627 }
1628 //
1629 if (!isOK) continue;
1630 nclusters=0;
1631 for (Int_t iLayer=0;iLayer<4;iLayer++){
1632 if (cseed[sLayer+iLayer].isOK()){
1633 nclusters+=cseed[sLayer+iLayer].fN2;
1634 }
1635 }
1636 //
1637 // iteration 0
1638 rieman.Reset();
1639 for (Int_t iLayer=0;iLayer<4;iLayer++){
1640 rieman.AddPoint(xcl[sLayer+iLayer],cseed[sLayer+iLayer].fYfitR[0],
1641 cseed[sLayer+iLayer].fZProb,1,10);
1642 }
1643 rieman.Update();
1644 //
1645 //
1646 chi2R =0; chi2Z=0;
1647 for (Int_t iLayer=0;iLayer<4;iLayer++){
1648 cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
1649 chi2R += (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0])*
1650 (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0]);
1651 cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1652 cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
1653 chi2Z += (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz)*
1654 (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz);
1655 cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1656 }
1657 Double_t curv = rieman.GetC();
1658 //
1659 // likelihoods
1660 //
1661 Double_t sumda =
1662 TMath::Abs(cseed[sLayer+0].fYfitR[1]- cseed[sLayer+0].fYref[1])+
1663 TMath::Abs(cseed[sLayer+1].fYfitR[1]- cseed[sLayer+1].fYref[1])+
1664 TMath::Abs(cseed[sLayer+2].fYfitR[1]- cseed[sLayer+2].fYref[1])+
1665 TMath::Abs(cseed[sLayer+3].fYfitR[1]- cseed[sLayer+3].fYref[1]);
1666 Double_t likea = TMath::Exp(-sumda*10.6);
1667 Double_t likechi2 = 0.0000000001;
1668 if (chi2R<0.5) likechi2+=TMath::Exp(-TMath::Sqrt(chi2R)*7.73);
1669 Double_t likechi2z = TMath::Exp(-chi2Z*0.088)/TMath::Exp(-chi2Z*0.019);
1670 Double_t likeN = TMath::Exp(-(72-nclusters)*0.19);
1671 Double_t like = likea*likechi2*likechi2z*likeN;
1672 //
1673 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fYref[1]-130*curv)*1.9);
1674 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fZref[1]-
1675 cseed[sLayer+0].fZref[0]/xcl[sLayer+0])*5.9);
1676 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
1677
1678 seedquality[registered] = like;
1679 seedlayer[registered] = sLayer;
1680 if (TMath::Log(0.000000000000001+like)<-15) continue;
1681 AliTRDseed seedb[6];
1682 for (Int_t iLayer=0;iLayer<6;iLayer++){
1683 seedb[iLayer] = cseed[iLayer];
1684 }
1685 //
1686 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1687 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1688 //<<<<<<<<<<<<<<< FULL TRACK FIT PART <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1689 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1690 //
1691 Int_t nlayers = 0;
1692 Int_t nusedf = 0;
1693 Int_t findable = 0;
1694 //
1695 // add new layers - avoid long extrapolation
1696 //
1697 Int_t tLayer[2]={0,0};
1698 if (sLayer==2) {tLayer[0]=1; tLayer[1]=0;}
1699 if (sLayer==1) {tLayer[0]=5; tLayer[1]=0;}
1700 if (sLayer==0) {tLayer[0]=4; tLayer[1]=5;}
1701 //
1702 for (Int_t iLayer=0;iLayer<2;iLayer++){
1703 Int_t jLayer = tLayer[iLayer]; // set tracking layer
1704 cseed[jLayer].Reset();
1705 cseed[jLayer].fTilt = hL[jLayer];
1706 cseed[jLayer].fPadLength = padlength[jLayer];
1707 cseed[jLayer].fX0 = xcl[jLayer];
1708 // get pad length and rough cluster
1709 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].fYref[0],
1710 cseed[jLayer].fZref[0],kRoad2y,kRoad2z);
1711 if (indexdummy<=0) continue;
1712 AliTRDcluster *cldummy = (AliTRDcluster*)GetCluster(indexdummy);
1713 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2()*12.);
1714 }
1715 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1716 //
1717 for (Int_t iLayer=0;iLayer<2;iLayer++){
1718 Int_t jLayer = tLayer[iLayer]; // set tracking layer
1719 if ( (jLayer==0) && !(cseed[1].isOK())) continue; // break not allowed
1720 if ( (jLayer==5) && !(cseed[4].isOK())) continue; // break not allowed
1721 Float_t zexp = cseed[jLayer].fZref[0];
1722 Double_t zroad = padlength[jLayer]*0.5+1.;
1723 //
1724 //
1725 for (Int_t iter=0;iter<2;iter++){
1726 AliTRDseed tseed = cseed[jLayer];
1727 Float_t quality = 10000;
1728 for (Int_t iTime=2;iTime<20;iTime++){
1729 AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1730 Double_t dxlayer = layer.GetX()-xcl[jLayer];
1731 Double_t yexp = tseed.fYref[0]+tseed.fYref[1]*dxlayer;
1732 Float_t yroad = kRoad1y;
1733 Int_t index = layer.FindNearestCluster(yexp,zexp, yroad, zroad);
1734 if (index<=0) continue;
1735 AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1736 //
1737 tseed.fIndexes[iTime] = index;
1738 tseed.fClusters[iTime] = cl; // register cluster
1739 tseed.fX[iTime] = dxlayer; // register cluster
1740 tseed.fY[iTime] = cl->GetY(); // register cluster
1741 tseed.fZ[iTime] = cl->GetZ(); // register cluster
1742 }
1743 tseed.Update();
1744 if (tseed.isOK()){
1745 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1746 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1747 TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
1748 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1749 //
1750 if (tquality<quality){
1751 cseed[jLayer]=tseed;
1752 quality = tquality;
1753 }
1754 }
1755 zroad*=2.;
1756 }
1757 if ( cseed[jLayer].isOK()){
1758 cseed[jLayer].CookLabels();
1759 cseed[jLayer].UpdateUsed();
1760 nusedf+= cseed[jLayer].fNUsed;
1761 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1762 }
1763 }
1764 //
1765 //
1766 // make copy
1767 AliTRDseed bseed[6];
1768 for (Int_t jLayer=0;jLayer<6;jLayer++){
1769 bseed[jLayer] = cseed[jLayer];
1770 }
1771 Float_t lastquality = 10000;
1772 Float_t lastchi2 = 10000;
1773 Float_t chi2 = 1000;
1774
1775 //
1776 for (Int_t iter =0; iter<4;iter++){
1777 //
1778 // sort tracklets according "quality", try to "improve" 4 worst
1779 //
1780 Float_t sumquality = 0;
1781 Float_t squality[6];
1782 Int_t sortindexes[6];
1783 for (Int_t jLayer=0;jLayer<6;jLayer++){
1784 if (bseed[jLayer].isOK()){
1785 AliTRDseed &tseed = bseed[jLayer];
1786 Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1787 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1788 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1789 TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
1790 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1791 squality[jLayer] = tquality;
1792 }
1793 else squality[jLayer]=-1;
1794 sumquality +=squality[jLayer];
1795 }
1796
1797 if (sumquality>=lastquality || chi2>lastchi2) break;
1798 lastquality = sumquality;
1799 lastchi2 = chi2;
1800 if (iter>0){
1801 for (Int_t jLayer=0;jLayer<6;jLayer++){
1802 cseed[jLayer] = bseed[jLayer];
1803 }
1804 }
1805 TMath::Sort(6,squality,sortindexes,kFALSE);
1806 //
1807 //
1808 for (Int_t jLayer=5;jLayer>1;jLayer--){
1809 Int_t bLayer = sortindexes[jLayer];
1810 AliTRDseed tseed = bseed[bLayer];
1811 for (Int_t iTime=2;iTime<20;iTime++){
1812 AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1813 Double_t dxlayer= layer.GetX()-xcl[bLayer];
1814 //
1815 Double_t zexp = tseed.fZref[0];
1816 Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1817 //
1818 Float_t roadz = padlength[bLayer]+1;
1819 if (TMath::Abs(tseed.fZProb-zexp)> padlength[bLayer]*0.5) {roadz = padlength[bLayer]*0.5;}
1820 if (tseed.fZfit[1]*tseed.fZref[1]<0) {roadz = padlength[bLayer]*0.5;}
1821 if (TMath::Abs(tseed.fZProb-zexp)<0.1*padlength[bLayer]) {
1822 zexp = tseed.fZProb;
1823 roadz = padlength[bLayer]*0.5;
1824 }
1825 //
1826 Double_t yexp = tseed.fYref[0]+
1827 tseed.fYref[1]*dxlayer-zcor;
1828 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
1829 if (index<=0) continue;
1830 AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1831 //
1832 tseed.fIndexes[iTime] = index;
1833 tseed.fClusters[iTime] = cl; // register cluster
1834 tseed.fX[iTime] = dxlayer; // register cluster
1835 tseed.fY[iTime] = cl->GetY(); // register cluster
1836 tseed.fZ[iTime] = cl->GetZ(); // register cluster
1837 }
1838 tseed.Update();
1839 if (tseed.isOK()) {
1840 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1841 Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1842 //
1843 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1844 TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
1845 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1846 //
1847 if (tquality<squality[bLayer])
1848 bseed[bLayer] = tseed;
1849 }
1850 }
1851 chi2 = AliTRDseed::FitRiemanTilt(bseed, kTRUE);
1852 }
1853 //
1854 //
1855 //
1856 nclusters = 0;
1857 nlayers = 0;
1858 findable = 0;
1859 for (Int_t iLayer=0;iLayer<6;iLayer++) {
1860 if (TMath::Abs(cseed[iLayer].fYref[0]/cseed[iLayer].fX0)<0.15)
1861 findable++;
1862 if (cseed[iLayer].isOK()){
1863 nclusters+=cseed[iLayer].fN2;
1864 nlayers++;
1865 }
1866 }
1867 if (nlayers<3) continue;
1868 rieman.Reset();
1869 for (Int_t iLayer=0;iLayer<6;iLayer++){
1870 if (cseed[iLayer].isOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
1871 cseed[iLayer].fZProb,1,10);
1872 }
1873 rieman.Update();
1874 //
1875 chi2RF =0;
1876 chi2ZF =0;
1877 for (Int_t iLayer=0;iLayer<6;iLayer++){
1878 if (cseed[iLayer].isOK()){
1879 cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
1880 chi2RF += (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0])*
1881 (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0]);
1882 cseed[iLayer].fYref[1] = rieman.GetDYat(xcl[iLayer]);
1883 cseed[iLayer].fZref[0] = rieman.GetZat(xcl[iLayer]);
1884 chi2ZF += (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz)*
1885 (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz);
1886 cseed[iLayer].fZref[1] = rieman.GetDZat(xcl[iLayer]);
1887 }
1888 }
1889 chi2RF/=TMath::Max((nlayers-3.),1.);
1890 chi2ZF/=TMath::Max((nlayers-3.),1.);
1891 curv = rieman.GetC();
1892
1893 //
1894
1895 Double_t xref2 = (xcl[2]+xcl[3])*0.5; // middle of the chamber
1896 Double_t dzmf = rieman.GetDZat(xref2);
1897 Double_t zmf = rieman.GetZat(xref2);
1898 //
1899 // fit hyperplane
1900 //
1901 Int_t npointsT =0;
1902 fitterTC.ClearPoints();
1903 fitterT2.ClearPoints();
1904 rieman2.Reset();
1905 for (Int_t iLayer=0; iLayer<6;iLayer++){
1906 if (!cseed[iLayer].isOK()) continue;
1907 for (Int_t itime=0;itime<25;itime++){
1908 if (!cseed[iLayer].fUsable[itime]) continue;
1909 Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2; // x relative to the midle chamber
1910 Double_t y = cseed[iLayer].fY[itime];
1911 Double_t z = cseed[iLayer].fZ[itime];
1912 // ExB correction to the correction
1913 // tilted rieman
1914 //
1915 Double_t uvt[6];
1916 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
1917 //
1918 Double_t t = 1./(x2*x2+y*y);
1919 uvt[1] = t; // t
1920 uvt[0] = 2.*x2*uvt[1]; // u
1921 //
1922 uvt[2] = 2.0*hL[iLayer]*uvt[1];
1923 uvt[3] = 2.0*hL[iLayer]*x*uvt[1];
1924 uvt[4] = 2.0*(y+hL[iLayer]*z)*uvt[1];
1925 //
1926 Double_t error = 2*0.2*uvt[1];
1927 fitterT2.AddPoint(uvt,uvt[4],error);
1928 //
1929 // constrained rieman
1930 //
1931 z =cseed[iLayer].fZ[itime];
1932 uvt[0] = 2.*x2*t; // u
1933 uvt[1] = 2*hL[iLayer]*x2*uvt[1];
1934 uvt[2] = 2*(y+hL[iLayer]*(z-GetZ()))*t;
1935 fitterTC.AddPoint(uvt,uvt[2],error);
1936 //
1937 rieman2.AddPoint(x2,y,z,1,10);
1938 npointsT++;
1939 }
1940 }
1941 rieman2.Update();
1942 fitterTC.Eval();
1943 fitterT2.Eval();
1944 Double_t rpolz0 = fitterT2.GetParameter(3);
1945 Double_t rpolz1 = fitterT2.GetParameter(4);
1946 //
1947 // linear fitter - not possible to make boundaries
1948 // non accept non possible z and dzdx combination
1949 //
1950 Bool_t acceptablez =kTRUE;
1951 for (Int_t iLayer=0; iLayer<6;iLayer++){
1952 if (cseed[iLayer].isOK()){
1953 Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
1954 if (TMath::Abs(cseed[iLayer].fZProb-zT2)>padlength[iLayer]*0.5+1)
1955 acceptablez = kFALSE;
1956 }
1957 }
1958 if (!acceptablez){
1959 fitterT2.FixParameter(3,zmf);
1960 fitterT2.FixParameter(4,dzmf);
1961 fitterT2.Eval();
1962 fitterT2.ReleaseParameter(3);
1963 fitterT2.ReleaseParameter(4);
1964 rpolz0 = fitterT2.GetParameter(3);
1965 rpolz1 = fitterT2.GetParameter(4);
1966 }
1967 //
1968 Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
1969 Double_t chi2TC = fitterTC.GetChisquare()/Float_t(npointsT);
1970 //
1971 Double_t polz1c = fitterTC.GetParameter(2);
1972 Double_t polz0c = polz1c*xref2;
1973 //
1974 Double_t aC = fitterTC.GetParameter(0);
1975 Double_t bC = fitterTC.GetParameter(1);
1976 Double_t CC = aC/TMath::Sqrt(bC*bC+1.); // curvature
1977 //
1978 Double_t aR = fitterT2.GetParameter(0);
1979 Double_t bR = fitterT2.GetParameter(1);
1980 Double_t dR = fitterT2.GetParameter(2);
1981 Double_t CR = 1+bR*bR-dR*aR;
1982 Double_t dca = 0.;
1983 if (CR>0){
1984 dca = -dR/(TMath::Sqrt(1+bR*bR-dR*aR)+TMath::Sqrt(1+bR*bR));
1985 CR = aR/TMath::Sqrt(CR);
1986 }
1987 //
1988 Double_t chi2ZT2=0, chi2ZTC=0;
1989 for (Int_t iLayer=0; iLayer<6;iLayer++){
1990 if (cseed[iLayer].isOK()){
1991 Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
1992 Double_t zTC = polz0c+polz1c*(xcl[iLayer] - xref2);
1993 chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz-zT2);
1994 chi2ZTC += TMath::Abs(cseed[iLayer].fMeanz-zTC);
1995 }
1996 }
1997 chi2ZT2/=TMath::Max((nlayers-3.),1.);
1998 chi2ZTC/=TMath::Max((nlayers-3.),1.);
1999 //
2000 //
2001 //
2002 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
2003 Float_t sumdaf = 0;
2004 for (Int_t iLayer=0;iLayer<6;iLayer++){
2005 if (cseed[iLayer].isOK())
2006 sumdaf += TMath::Abs((cseed[iLayer].fYfit[1]-cseed[iLayer].fYref[1])/cseed[iLayer].fSigmaY2);
2007 }
2008 sumdaf /= Float_t (nlayers-2.);
2009 //
2010 // likelihoods for full track
2011 //
2012 Double_t likezf = TMath::Exp(-chi2ZF*0.14);
2013 Double_t likechi2C = TMath::Exp(-chi2TC*0.677);
2014 Double_t likechi2TR = TMath::Exp(-chi2TR*0.78);
2015 Double_t likeaf = TMath::Exp(-sumdaf*3.23);
2016 seedquality2[registered] = likezf*likechi2TR*likeaf;
2017// Bool_t isGold = kFALSE;
2018//
2019// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
2020// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
2021// if (isGold &&nusedf<10){
2022// for (Int_t jLayer=0;jLayer<6;jLayer++){
2023// if ( seed[index][jLayer].isOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
2024// seed[index][jLayer].UseClusters(); //sign gold
2025// }
2026// }
2027 //
2028 //
2029 //
2030 Int_t index0=0;
2031 if (!cseed[0].isOK()){
2032 index0 = 1;
2033 if (!cseed[1].isOK()) index0 = 2;
2034 }
2035 seedparams[registered][0] = cseed[index0].fX0;
2036 seedparams[registered][1] = cseed[index0].fYref[0];
2037 seedparams[registered][2] = cseed[index0].fZref[0];
2038 seedparams[registered][5] = CR;
2039 seedparams[registered][3] = cseed[index0].fX0*CR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
2040 seedparams[registered][4] = cseed[index0].fZref[1]/
2041 TMath::Sqrt(1+cseed[index0].fYref[1]*cseed[index0].fYref[1]);
2042 seedparams[registered][6] = ns;
2043 //
2044 //
2045 Int_t labels[12], outlab[24];
2046 Int_t nlab=0;
2047 for (Int_t iLayer=0;iLayer<6;iLayer++){
2048 if (!cseed[iLayer].isOK()) continue;
2049 if (cseed[iLayer].fLabels[0]>=0) {
2050 labels[nlab] = cseed[iLayer].fLabels[0];
2051 nlab++;
2052 }
2053 if (cseed[iLayer].fLabels[1]>=0) {
2054 labels[nlab] = cseed[iLayer].fLabels[1];
2055 nlab++;
2056 }
2057 }
2058 Freq(nlab,labels,outlab,kFALSE);
2059 Int_t label = outlab[0];
2060 Int_t frequency = outlab[1];
2061 for (Int_t iLayer=0;iLayer<6;iLayer++){
2062 cseed[iLayer].fFreq = frequency;
2063 cseed[iLayer].fC = CR;
2064 cseed[iLayer].fCC = CC;
2065 cseed[iLayer].fChi2 = chi2TR;
2066 cseed[iLayer].fChi2Z = chi2ZF;
2067 }
2068 //
2069 if (1||(!isFake)){ //debugging print
2070 Float_t zvertex = GetZ();
2071 TTreeSRedirector& cstream = *fDebugStreamer;
2072 cstream<<"Seeds1"<<
2073 "isFake="<<isFake<<
2074 "Vertex="<<zvertex<<
2075 "Rieman2.="<<&rieman2<<
2076 "Rieman.="<<&rieman<<
2077 "Xref="<<xref<<
2078 "X0="<<xcl[0]<<
2079 "X1="<<xcl[1]<<
2080 "X2="<<xcl[2]<<
2081 "X3="<<xcl[3]<<
2082 "X4="<<xcl[4]<<
2083 "X5="<<xcl[5]<<
2084 "Chi2R="<<chi2R<<
2085 "Chi2Z="<<chi2Z<<
2086 "Chi2RF="<<chi2RF<< //chi2 of trackletes on full track
2087 "Chi2ZF="<<chi2ZF<< //chi2 z on tracklets on full track
2088 "Chi2ZT2="<<chi2ZT2<< //chi2 z on tracklets on full track - rieman tilt
2089 "Chi2ZTC="<<chi2ZTC<< //chi2 z on tracklets on full track - rieman tilt const
2090 //
2091 "Chi2TR="<<chi2TR<< //chi2 without vertex constrain
2092 "Chi2TC="<<chi2TC<< //chi2 with vertex constrain
2093 "C="<<curv<< // non constrained - no tilt correction
2094 "DR="<<dR<< // DR parameter - tilt correction
2095 "DCA="<<dca<< // DCA - tilt correction
2096 "CR="<<CR<< // non constrained curvature - tilt correction
2097 "CC="<<CC<< // constrained curvature
2098 "Polz0="<<polz0c<<
2099 "Polz1="<<polz1c<<
2100 "RPolz0="<<rpolz0<<
2101 "RPolz1="<<rpolz1<<
2102 "Ncl="<<nclusters<<
2103 "Nlayers="<<nlayers<<
2104 "NUsedS="<<nusedCl<<
2105 "NUsed="<<nusedf<<
2106 "Findable="<<findable<<
2107 "Like="<<like<<
2108 "LikePrim="<<likePrim<<
2109 "Likechi2C="<<likechi2C<<
2110 "Likechi2TR="<<likechi2TR<<
2111 "Likezf="<<likezf<<
2112 "LikeF="<<seedquality2[registered]<<
2113 "S0.="<<&cseed[0]<<
2114 "S1.="<<&cseed[1]<<
2115 "S2.="<<&cseed[2]<<
2116 "S3.="<<&cseed[3]<<
2117 "S4.="<<&cseed[4]<<
2118 "S5.="<<&cseed[5]<<
2119 "SB0.="<<&seedb[0]<<
2120 "SB1.="<<&seedb[1]<<
2121 "SB2.="<<&seedb[2]<<
2122 "SB3.="<<&seedb[3]<<
2123 "SB4.="<<&seedb[4]<<
2124 "SB5.="<<&seedb[5]<<
2125 "Label="<<label<<
2126 "Freq="<<frequency<<
2127 "sLayer="<<sLayer<<
2128 "\n";
2129 }
2130 if (registered<maxseed-1) {
2131 registered++;
2132 cseed = seed[registered];
2133 }
2134 }// end of loop over layer 1
2135 } // end of loop over layer 0
2136 } // end of loop over layer 3
2137 } // end of loop over seeding time bins
2138 //
2139 // choos best
2140 //
2141 TMath::Sort(registered,seedquality2,sort,kTRUE);
2142 Bool_t signedseed[maxseed];
2143 for (Int_t i=0;i<registered;i++){
2144 signedseed[i]= kFALSE;
2145 }
2146 for (Int_t iter=0; iter<5; iter++){
2147 for (Int_t iseed=0;iseed<registered;iseed++){
2148 Int_t index = sort[iseed];
2149 if (signedseed[index]) continue;
2150 Int_t labelsall[1000];
2151 Int_t nlabelsall=0;
2152 Int_t naccepted=0;;
2153 Int_t sLayer = seedlayer[index];
2154 Int_t ncl = 0;
2155 Int_t nused = 0;
2156 Int_t nlayers =0;
2157 Int_t findable = 0;
2158 for (Int_t jLayer=0;jLayer<6;jLayer++){
2159 if (TMath::Abs(seed[index][jLayer].fYref[0]/xcl[jLayer])<0.15)
2160 findable++;
2161 if (seed[index][jLayer].isOK()){
2162 seed[index][jLayer].UpdateUsed();
2163 ncl +=seed[index][jLayer].fN2;
2164 nused +=seed[index][jLayer].fNUsed;
2165 nlayers++;
2166 //cooking label
2167 for (Int_t itime=0;itime<25;itime++){
2168 if (seed[index][jLayer].fUsable[itime]){
2169 naccepted++;
2170 for (Int_t ilab=0;ilab<3;ilab++){
2171 Int_t tindex = seed[index][jLayer].fClusters[itime]->GetLabel(ilab);
2172 if (tindex>=0){
2173 labelsall[nlabelsall] = tindex;
2174 nlabelsall++;
2175 }
2176 }
2177 }
2178 }
2179 }
2180 }
7ad19338 2181 //
69b55c55 2182 if (nused>30) continue;
7ad19338 2183 //
69b55c55 2184 if (iter==0){
2185 if (nlayers<6) continue;
2186 if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue; // gold
2187 }
2188 //
2189 if (iter==1){
2190 if (nlayers<findable) continue;
2191 if (TMath::Log(0.000000001+seedquality2[index])<-4.) continue; //
7ad19338 2192 }
7ad19338 2193 //
7ad19338 2194 //
69b55c55 2195 if (iter==2){
2196 if (nlayers==findable || nlayers==6) continue;
2197 if (TMath::Log(0.000000001+seedquality2[index])<-6.) continue;
2198 }
7ad19338 2199 //
69b55c55 2200 if (iter==3){
2201 if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue;
2202 }
7ad19338 2203 //
69b55c55 2204 if (iter==4){
2205 if (TMath::Log(0.000000001+seedquality2[index])-nused/(nlayers-3.)<-15.) continue;
2206 }
7ad19338 2207 //
69b55c55 2208 signedseed[index] = kTRUE;
2209 //
2210 Int_t labels[1000], outlab[1000];
2211 Int_t nlab=0;
2212 for (Int_t iLayer=0;iLayer<6;iLayer++){
2213 if (seed[index][iLayer].isOK()){
2214 if (seed[index][iLayer].fLabels[0]>=0) {
2215 labels[nlab] = seed[index][iLayer].fLabels[0];
2216 nlab++;
2217 }
2218 if (seed[index][iLayer].fLabels[1]>=0) {
2219 labels[nlab] = seed[index][iLayer].fLabels[1];
2220 nlab++;
2221 }
2222 }
7ad19338 2223 }
69b55c55 2224 Freq(nlab,labels,outlab,kFALSE);
2225 Int_t label = outlab[0];
2226 Int_t frequency = outlab[1];
2227 Freq(nlabelsall,labelsall,outlab,kFALSE);
2228 Int_t label1 = outlab[0];
2229 Int_t label2 = outlab[2];
2230 Float_t fakeratio = (naccepted-outlab[1])/Float_t(naccepted);
2231 Float_t ratio = Float_t(nused)/Float_t(ncl);
2232 if (ratio<0.25){
2233 for (Int_t jLayer=0;jLayer<6;jLayer++){
2234 if ( seed[index][jLayer].isOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.2 )
2235 seed[index][jLayer].UseClusters(); //sign gold
2236 }
7ad19338 2237 }
2238 //
69b55c55 2239 Int_t eventNr = esd->GetEventNumber();
2240 TTreeSRedirector& cstream = *fDebugStreamer;
2241 //
2242 // register seed
2243 //
2244 AliTRDtrack * track = RegisterSeed(seed[index],seedparams[index]);
2245 AliTRDtrack dummy;
2246 if (!track) track=&dummy;
2247 else{
2248 AliESDtrack esdtrack;
2249 esdtrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
2250 esdtrack.SetLabel(label);
2251 esd->AddTrack(&esdtrack);
7ad19338 2252 TTreeSRedirector& cstream = *fDebugStreamer;
69b55c55 2253 cstream<<"Tracks"<<
2254 "EventNr="<<eventNr<<
2255 "ESD.="<<&esdtrack<<
2256 "trd.="<<track<<
2257 "trdback.="<<track<<
7ad19338 2258 "\n";
2259 }
69b55c55 2260
2261 cstream<<"Seeds2"<<
2262 "Iter="<<iter<<
2263 "Track.="<<track<<
2264 "Like="<<seedquality[index]<<
2265 "LikeF="<<seedquality2[index]<<
2266 "S0.="<<&seed[index][0]<<
2267 "S1.="<<&seed[index][1]<<
2268 "S2.="<<&seed[index][2]<<
2269 "S3.="<<&seed[index][3]<<
2270 "S4.="<<&seed[index][4]<<
2271 "S5.="<<&seed[index][5]<<
2272 "Label="<<label<<
2273 "Label1="<<label1<<
2274 "Label2="<<label2<<
2275 "FakeRatio="<<fakeratio<<
2276 "Freq="<<frequency<<
2277 "Ncl="<<ncl<<
2278 "Nlayers="<<nlayers<<
2279 "Findable="<<findable<<
2280 "NUsed="<<nused<<
2281 "sLayer="<<sLayer<<
2282 "EventNr="<<eventNr<<
2283 "\n";
7ad19338 2284 }
2285 }
69b55c55 2286 } // end of loop over sectors
2287 delete [] pseed;
2288}
2289
5443e65e 2290//_____________________________________________________________________________
b7a0917f 2291Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const
5443e65e 2292{
2293 //
a819a5f7 2294 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2295 // from the file. The names of the cluster tree and branches
2296 // should match the ones used in AliTRDclusterizer::WriteClusters()
2297 //
4f1c04d3 2298 Int_t nsize = Int_t(ClusterTree->GetTotBytes()/(sizeof(AliTRDcluster)));
2299 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2300
c630aafd 2301 TBranch *branch=ClusterTree->GetBranch("TRDcluster");
2302 if (!branch) {
2303 Error("ReadClusters","Can't get the branch !");
2304 return 1;
2305 }
029cd327 2306 branch->SetAddress(&clusterArray);
5443e65e 2307
2308 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
19dd5b2f 2309 // printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
a819a5f7 2310
a819a5f7 2311 // Loop through all entries in the tree
eb187bed 2312 Int_t nbytes = 0;
a819a5f7 2313 AliTRDcluster *c = 0;
7bed16a7 2314 // printf("\n");
a819a5f7 2315 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2316
2317 // Import the tree
5443e65e 2318 nbytes += ClusterTree->GetEvent(iEntry);
2319
a819a5f7 2320 // Get the number of points in the detector
029cd327 2321 Int_t nCluster = clusterArray->GetEntriesFast();
e24ea474 2322// printf("\r Read %d clusters from entry %d", nCluster, iEntry);
5443e65e 2323
a819a5f7 2324 // Loop through all TRD digits
2325 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
029cd327 2326 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
4f1c04d3 2327 AliTRDcluster *co = c;
a819a5f7 2328 array->AddLast(co);
4f1c04d3 2329 // delete clusterArray->RemoveAt(iCluster);
2330 clusterArray->RemoveAt(iCluster);
a819a5f7 2331 }
2332 }
7c1698cb 2333// cout<<"Allocated"<<nsize<<"\tLoaded"<<array->GetEntriesFast()<<"\n";
a819a5f7 2334
029cd327 2335 delete clusterArray;
5443e65e 2336
c630aafd 2337 return 0;
a819a5f7 2338}
2339
3551db50 2340//__________________________________________________________________
2341Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint& p) const
2342{
2343 //
2344 // Get track space point with index i
2345 // Origin: C.Cheshkov
2346 //
2347
2348 AliTRDcluster *cl = (AliTRDcluster*)fClusters->UncheckedAt(index);
2349 Int_t idet = cl->GetDetector();
2350 Int_t isector = fGeom->GetSector(idet);
2351 Int_t ichamber= fGeom->GetChamber(idet);
2352 Int_t iplan = fGeom->GetPlane(idet);
2353 Double_t local[3];
2354 local[0]=GetX(isector,iplan,cl->GetLocalTimeBin());
2355 local[1]=cl->GetY();
2356 local[2]=cl->GetZ();
2357 Double_t global[3];
2358 fGeom->RotateBack(idet,local,global);
2359 p.SetXYZ(global[0],global[1],global[2]);
2360 AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
2361 switch (iplan) {
2362 case 0:
2363 iLayer = AliAlignObj::kTRD1;
2364 break;
2365 case 1:
2366 iLayer = AliAlignObj::kTRD2;
2367 break;
2368 case 2:
2369 iLayer = AliAlignObj::kTRD3;
2370 break;
2371 case 3:
2372 iLayer = AliAlignObj::kTRD4;
2373 break;
2374 case 4:
2375 iLayer = AliAlignObj::kTRD5;
2376 break;
2377 case 5:
2378 iLayer = AliAlignObj::kTRD6;
2379 break;
2380 };
2381 Int_t modId = isector*fGeom->Ncham()+ichamber;
2382 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
2383 p.SetVolumeID(volid);
2384
2385 return kTRUE;
2386
2387}
2388
46d29e70 2389//__________________________________________________________________
029cd327 2390void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
2391{
2392 //
2393 // This cooks a label. Mmmmh, smells good...
2394 //
46d29e70 2395
2396 Int_t label=123456789, index, i, j;
5443e65e 2397 Int_t ncl=pt->GetNumberOfClusters();
029cd327 2398 const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
5443e65e 2399
029cd327 2400 Bool_t labelAdded;
46d29e70 2401
029cd327 2402 // Int_t s[kRange][2];
2403 Int_t **s = new Int_t* [kRange];
2404 for (i=0; i<kRange; i++) {
d1b06c24 2405 s[i] = new Int_t[2];
2406 }
029cd327 2407 for (i=0; i<kRange; i++) {
46d29e70 2408 s[i][0]=-1;
2409 s[i][1]=0;
2410 }
2411
2412 Int_t t0,t1,t2;
2413 for (i=0; i<ncl; i++) {
5443e65e 2414 index=pt->GetClusterIndex(i);
46d29e70 2415 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
5443e65e 2416 t0=c->GetLabel(0);
2417 t1=c->GetLabel(1);
2418 t2=c->GetLabel(2);
46d29e70 2419 }
2420
2421 for (i=0; i<ncl; i++) {
5443e65e 2422 index=pt->GetClusterIndex(i);
46d29e70 2423 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2424 for (Int_t k=0; k<3; k++) {
5443e65e 2425 label=c->GetLabel(k);
029cd327 2426 labelAdded=kFALSE; j=0;
46d29e70 2427 if (label >= 0) {
029cd327 2428 while ( (!labelAdded) && ( j < kRange ) ) {
a9814c09 2429 if (s[j][0]==label || s[j][1]==0) {
2430 s[j][0]=label;
2431 s[j][1]=s[j][1]+1;
029cd327 2432 labelAdded=kTRUE;
a9814c09 2433 }
2434 j++;
2435 }
46d29e70 2436 }
2437 }
2438 }
2439
46d29e70 2440 Int_t max=0;
2441 label = -123456789;
2442
029cd327 2443 for (i=0; i<kRange; i++) {
46d29e70 2444 if (s[i][1]>max) {
2445 max=s[i][1]; label=s[i][0];
2446 }
2447 }
5443e65e 2448
029cd327 2449 for (i=0; i<kRange; i++) {
5443e65e 2450 delete []s[i];
2451 }
2452
d1b06c24 2453 delete []s;
5443e65e 2454
2455 if ((1.- Float_t(max)/ncl) > wrong) label=-label;
2456
2457 pt->SetLabel(label);
2458
46d29e70 2459}
2460
c630aafd 2461
5443e65e 2462//__________________________________________________________________
029cd327 2463void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
2464{
2465 //
2466 // Use clusters, but don't abuse them!
2467 //
69b55c55 2468 const Float_t kmaxchi2 =18;
2469 const Float_t kmincl =10;
2470 AliTRDtrack * track = (AliTRDtrack*)t;
2471 //
5443e65e 2472 Int_t ncl=t->GetNumberOfClusters();
2473 for (Int_t i=from; i<ncl; i++) {
2474 Int_t index = t->GetClusterIndex(i);
2475 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
69b55c55 2476 //
2477 Int_t iplane = fGeom->GetPlane(c->GetDetector());
2478 if (track->fTracklets[iplane].GetChi2()>kmaxchi2) continue;
2479 if (track->fTracklets[iplane].GetN()<kmincl) continue;
2480 if (!(c->IsUsed())) c->Use();
5443e65e 2481 }
2482}
2483
2484
2485//_____________________________________________________________________
029cd327 2486Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2487{
2488 // Parametrised "expected" error of the cluster reconstruction in Y
2489
2490 Double_t s = 0.08 * 0.08;
2491 return s;
2492}
2493
2494//_____________________________________________________________________
029cd327 2495Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
0a29d0f1 2496{
5443e65e 2497 // Parametrised "expected" error of the cluster reconstruction in Z
2498
a9814c09 2499 Double_t s = 9 * 9 /12.;
5443e65e 2500 return s;
2501}
2502
5443e65e 2503//_____________________________________________________________________
029cd327 2504Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2505{
2506 //
029cd327 2507 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2508 // in tracking sector <sector> and plane <plane>
2509 //
2510
029cd327 2511 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
5443e65e 2512 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2513 return fTrSec[sector]->GetLayer(pl)->GetX();
2514
2515}
2516
c630aafd 2517
5443e65e 2518//_______________________________________________________
2519AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
029cd327 2520 Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex)
5443e65e 2521{
0a29d0f1 2522 //
5443e65e 2523 // AliTRDpropagationLayer constructor
0a29d0f1 2524 //
46d29e70 2525
029cd327 2526 fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = radLength;
2527 fClusters = NULL; fIndex = NULL; fTimeBinIndex = tbIndex;
5443e65e 2528
46d29e70 2529
029cd327 2530 for(Int_t i=0; i < (Int_t) kZones; i++) {
5443e65e 2531 fZc[i]=0; fZmax[i] = 0;
a819a5f7 2532 }
5443e65e 2533
2534 fYmax = 0;
2535
2536 if(fTimeBinIndex >= 0) {
029cd327 2537 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
2538 fIndex = new UInt_t[kMaxClusterPerTimeBin];
a819a5f7 2539 }
46d29e70 2540
3c625a9b 2541 for (Int_t i=0;i<5;i++) fIsHole[i] = kFALSE;
5443e65e 2542 fHole = kFALSE;
2543 fHoleZc = 0;
2544 fHoleZmax = 0;
2545 fHoleYc = 0;
2546 fHoleYmax = 0;
2547 fHoleRho = 0;
2548 fHoleX0 = 0;
2549
2550}
2551
2552//_______________________________________________________
2553void AliTRDtracker::AliTRDpropagationLayer::SetHole(
a9814c09 2554 Double_t Zmax, Double_t Ymax, Double_t rho,
029cd327 2555 Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2556{
2557 //
2558 // Sets hole in the layer
2559 //
5443e65e 2560 fHole = kTRUE;
2561 fHoleZc = Zc;
2562 fHoleZmax = Zmax;
2563 fHoleYc = Yc;
2564 fHoleYmax = Ymax;
2565 fHoleRho = rho;
029cd327 2566 fHoleX0 = radLength;
5443e65e 2567}
2568
46d29e70 2569
5443e65e 2570//_______________________________________________________
2571AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
2572{
2573 //
2574 // AliTRDtrackingSector Constructor
2575 //
a5cadd36 2576 AliTRDpadPlane *padPlane = 0;
2577
5443e65e 2578 fGeom = geo;
2579 fPar = par;
2580 fGeomSector = gs;
5443e65e 2581 fN = 0;
3c625a9b 2582 //
2583 // get holes description from geometry
2584 Bool_t holes[AliTRDgeometry::kNcham];
2585 //printf("sector\t%d\t",gs);
2586 for (Int_t icham=0; icham<AliTRDgeometry::kNcham;icham++){
2587 holes[icham] = fGeom->IsHole(0,icham,gs);
2588 //printf("%d",holes[icham]);
2589 }
2590 //printf("\n");
2591
029cd327 2592 for(UInt_t i=0; i < kMaxTimeBinIndex; i++) fTimeBinIndex[i] = -1;
5443e65e 2593
2594
2595 AliTRDpropagationLayer* ppl;
2596
029cd327 2597 Double_t x, xin, xout, dx, rho, radLength;
5443e65e 2598 Int_t steps;
46d29e70 2599
5443e65e 2600 // set time bins in the gas of the TPC
46d29e70 2601
5443e65e 2602 xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
029cd327 2603 rho = 0.9e-3; radLength = 28.94;
5443e65e 2604
2605 for(Int_t i=0; i<steps; i++) {
2606 x = xin + i*dx + dx/2;
029cd327 2607 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2608 InsertLayer(ppl);
46d29e70 2609 }
2610
5443e65e 2611 // set time bins in the outer field cage vessel
46d29e70 2612
029cd327 2613 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2614 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2615 InsertLayer(ppl);
46d29e70 2616
029cd327 2617 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2618 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2619 InsertLayer(ppl);
2620
029cd327 2621 dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
5443e65e 2622 steps = 5; dx = (xout - xin)/steps;
2623 for(Int_t i=0; i<steps; i++) {
2624 x = xin + i*dx + dx/2;
029cd327 2625 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2626 InsertLayer(ppl);
2627 }
2628
029cd327 2629 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2630 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2631 InsertLayer(ppl);
2632
029cd327 2633 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2634 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2635 InsertLayer(ppl);
0a29d0f1 2636
5443e65e 2637
2638 // set time bins in CO2
2639
2640 xin = xout; xout = 275.0;
2641 steps = 50; dx = (xout - xin)/steps;
029cd327 2642 rho = 1.977e-3; radLength = 36.2;
5443e65e 2643
2644 for(Int_t i=0; i<steps; i++) {
2645 x = xin + i*dx + dx/2;
029cd327 2646 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2647 InsertLayer(ppl);
2648 }
2649
2650 // set time bins in the outer containment vessel
2651
029cd327 2652 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2653 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2654 InsertLayer(ppl);
2655
029cd327 2656 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2657 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2658 InsertLayer(ppl);
2659
029cd327 2660 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2661 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2662 InsertLayer(ppl);
2663
029cd327 2664 dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
5443e65e 2665 steps = 10; dx = (xout - xin)/steps;
2666 for(Int_t i=0; i<steps; i++) {
2667 x = xin + i*dx + dx/2;
029cd327 2668 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2669 InsertLayer(ppl);
2670 }
2671
029cd327 2672 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2673 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2674 InsertLayer(ppl);
2675
029cd327 2676 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2677 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2678 InsertLayer(ppl);
2679
029cd327 2680 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2681 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
5443e65e 2682 InsertLayer(ppl);
2683
2684 Double_t xtrd = (Double_t) fGeom->Rmin();
2685
2686 // add layers between TPC and TRD (Air temporarily)
2687 xin = xout; xout = xtrd;
2688 steps = 50; dx = (xout - xin)/steps;
029cd327 2689 rho = 1.2e-3; radLength = 36.66;
5443e65e 2690
2691 for(Int_t i=0; i<steps; i++) {
2692 x = xin + i*dx + dx/2;
029cd327 2693 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2694 InsertLayer(ppl);
2695 }
2696
2697
3c625a9b 2698 // Double_t alpha=AliTRDgeometry::GetAlpha();
5443e65e 2699
2700 // add layers for each of the planes
2701
2702 Double_t dxRo = (Double_t) fGeom->CroHght(); // Rohacell
2703 Double_t dxSpace = (Double_t) fGeom->Cspace(); // Spacing between planes
2704 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2705 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2706 Double_t dxRad = (Double_t) fGeom->CraHght(); // Radiator
2707 Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
2708 Double_t dxPlane = dxTEC + dxSpace;
2709
029cd327 2710 Int_t tb, tbIndex;
2711 const Int_t kNchambers = AliTRDgeometry::Ncham();
3c625a9b 2712 Double_t ymax = 0;
2713 //, holeYmax = 0;
2714 Double_t ymaxsensitive=0;
029cd327 2715 Double_t *zc = new Double_t[kNchambers];
2716 Double_t *zmax = new Double_t[kNchambers];
3c625a9b 2717 Double_t *zmaxsensitive = new Double_t[kNchambers];
2718 // Double_t holeZmax = 1000.; // the whole sector is missing
5443e65e 2719
3551db50 2720 AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
2721 if (!commonParam)
2722 {
2723 printf("<AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector> ");
2724 printf("Could not get common params\n");
2725 return;
2726 }
2727
5443e65e 2728 for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
3c625a9b 2729 //
5443e65e 2730 // Radiator
2731 xin = xtrd + plane * dxPlane; xout = xin + dxRad;
029cd327 2732 steps = 12; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
5443e65e 2733 for(Int_t i=0; i<steps; i++) {
2734 x = xin + i*dx + dx/2;
3c625a9b 2735 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2736 InsertLayer(ppl);
2737 }
2738
3c625a9b 2739 ymax = fGeom->GetChamberWidth(plane)/2.;
a5cadd36 2740 // Modidified for new pad plane class, 22.04.05 (C.B.)
2741 // ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
3551db50 2742 padPlane = commonParam->GetPadPlane(plane,0);
a5cadd36 2743 ymaxsensitive = (padPlane->GetColSize(1)*padPlane->GetNcols()-4)/2.;
7ad19338 2744
2745 // ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
3c625a9b 2746
029cd327 2747 for(Int_t ch = 0; ch < kNchambers; ch++) {
2748 zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
a5cadd36 2749 //
2750 // Modidified for new pad plane class, 22.04.05 (C.B.)
2751 //Float_t pad = fPar->GetRowPadSize(plane,ch,0);
2752 Float_t pad = padPlane->GetRowSize(1);
7ad19338 2753 //Float_t pad = fPar->GetRowPadSize(plane,ch,0);
3551db50 2754 Float_t row0 = commonParam->GetRow0(plane,ch,0);
2755 Int_t nPads = commonParam->GetRowMax(plane,ch,0);
857b3eb0 2756 zmaxsensitive[ch] = Float_t(nPads)*pad/2.;
3c625a9b 2757 // zc[ch] = (pad * nPads)/2 + row0 - pad/2;
7ad19338 2758 // zc[ch] = (pad * nPads)/2 + row0;
2759 zc[ch] = -(pad * nPads)/2 + row0;
3c625a9b 2760 //zc[ch] = row0+zmax[ch]-AliTRDgeometry::RpadW();
2761
5443e65e 2762 }
2763
7ad19338 2764 dx = fgkDriftCorrection*fPar->GetDriftVelocity()
3551db50 2765 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
029cd327 2766 rho = 0.00295 * 0.85; radLength = 11.0;
5443e65e 2767
3551db50 2768 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
5443e65e 2769 Double_t xbottom = x0 - dxDrift;
2770 Double_t xtop = x0 + dxAmp;
3c625a9b 2771 //
5443e65e 2772 // Amplification region
5443e65e 2773 steps = (Int_t) (dxAmp/dx);
2774
2775 for(tb = 0; tb < steps; tb++) {
7ad19338 2776 x = x0 + tb * dx + dx/2+ fgkOffsetX;
029cd327 2777 tbIndex = CookTimeBinIndex(plane, -tb-1);
2778 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3c625a9b 2779 ppl->SetYmax(ymax,ymaxsensitive);
2780 ppl->SetZ(zc, zmax, zmaxsensitive);
2781 ppl->SetHoles(holes);
5443e65e 2782 InsertLayer(ppl);
2783 }
029cd327 2784 tbIndex = CookTimeBinIndex(plane, -steps);
5443e65e 2785 x = (x + dx/2 + xtop)/2;
2786 dx = 2*(xtop-x);
029cd327 2787 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3c625a9b 2788 ppl->SetYmax(ymax,ymaxsensitive);
2789 ppl->SetZ(zc, zmax,zmaxsensitive);
2790 ppl->SetHoles(holes);
5443e65e 2791 InsertLayer(ppl);
2792
2793 // Drift region
7ad19338 2794
2795 dx = fgkDriftCorrection*fPar->GetDriftVelocity()
3551db50 2796 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
8979685e 2797 steps = (Int_t) (dxDrift/dx)+3;
5443e65e 2798
2799 for(tb = 0; tb < steps; tb++) {
7ad19338 2800 x = x0 - tb * dx - dx/2 + fgkOffsetX; //temporary fix - fix it the parameters
029cd327 2801 tbIndex = CookTimeBinIndex(plane, tb);
5443e65e 2802
029cd327 2803 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3c625a9b 2804 ppl->SetYmax(ymax,ymaxsensitive);
2805 ppl->SetZ(zc, zmax, zmaxsensitive);
2806 ppl->SetHoles(holes);
5443e65e 2807 InsertLayer(ppl);
2808 }
029cd327 2809 tbIndex = CookTimeBinIndex(plane, steps);
5443e65e 2810 x = (x - dx/2 + xbottom)/2;
2811 dx = 2*(x-xbottom);
029cd327 2812 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3c625a9b 2813 ppl->SetYmax(ymax,ymaxsensitive);
2814 ppl->SetZ(zc, zmax, zmaxsensitive);
2815 ppl->SetHoles(holes);
5443e65e 2816 InsertLayer(ppl);
2817
2818 // Pad Plane
029cd327 2819 xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; radLength = 33.0;
2820 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3c625a9b 2821 ppl->SetYmax(ymax,ymaxsensitive);
2822 ppl->SetZ(zc, zmax,zmax);
2823 ppl->SetHoles(holes);
5443e65e 2824 InsertLayer(ppl);
2825
2826 // Rohacell
2827 xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
029cd327 2828 steps = 5; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
5443e65e 2829 for(Int_t i=0; i<steps; i++) {
2830 x = xin + i*dx + dx/2;
029cd327 2831 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3c625a9b 2832 ppl->SetYmax(ymax,ymaxsensitive);
2833 ppl->SetZ(zc, zmax,zmax);
2834 ppl->SetHoles(holes);
5443e65e 2835 InsertLayer(ppl);
2836 }
2837
2838 // Space between the chambers, air
2839 xin = xout; xout = xtrd + (plane + 1) * dxPlane;
029cd327 2840 steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
5443e65e 2841 for(Int_t i=0; i<steps; i++) {
2842 x = xin + i*dx + dx/2;
029cd327 2843 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2844 InsertLayer(ppl);
2845 }
2846 }
2847
2848 // Space between the TRD and RICH
2849 Double_t xRICH = 500.;
2850 xin = xout; xout = xRICH;
029cd327 2851 steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
5443e65e 2852 for(Int_t i=0; i<steps; i++) {
2853 x = xin + i*dx + dx/2;
029cd327 2854 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
5443e65e 2855 InsertLayer(ppl);
2856 }
2857
2858 MapTimeBinLayers();
029cd327 2859 delete [] zc;
2860 delete [] zmax;
4f1c04d3 2861 delete [] zmaxsensitive;
5443e65e 2862
2863}
2864
2865//______________________________________________________
2866
029cd327 2867Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 2868{
2869 //
2870 // depending on the digitization parameters calculates "global"
029cd327 2871 // time bin index for timebin <localTB> in plane <plane>
5443e65e 2872 //
2873
2874 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2875 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
7ad19338 2876
2877 Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
3551db50 2878 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
5443e65e 2879
2880 Int_t tbAmp = fPar->GetTimeBefore();
2881 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
b3a5a838 2882 if(kTRUE) maxAmp = 0; // intentional until we change parameter class
5443e65e 2883 Int_t tbDrift = fPar->GetTimeMax();
8979685e 2884 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx)+4; // MI change - take also last time bins
5443e65e 2885
029cd327 2886 Int_t tbPerPlane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
5443e65e 2887
029cd327 2888 Int_t gtb = (plane+1) * tbPerPlane - localTB - 1 - TMath::Min(tbAmp,maxAmp);
5443e65e 2889
029cd327 2890 if((localTB < 0) &&
2891 (TMath::Abs(localTB) > TMath::Min(tbAmp,maxAmp))) return -1;
2892 if(localTB >= TMath::Min(tbDrift,maxDrift)) return -1;
5443e65e 2893
2894 return gtb;
2895
2896
2897}
2898
2899//______________________________________________________
2900
2901void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2902{
2903 //
2904 // For all sensitive time bins sets corresponding layer index
2905 // in the array fTimeBins
2906 //
2907
2908 Int_t index;
2909
2910 for(Int_t i = 0; i < fN; i++) {
2911 index = fLayers[i]->GetTimeBinIndex();
2912
2913 // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2914
2915 if(index < 0) continue;
029cd327 2916 if(index >= (Int_t) kMaxTimeBinIndex) {
5443e65e 2917 printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2918 printf(" index %d exceeds allowed maximum of %d!\n",
029cd327 2919 index, kMaxTimeBinIndex-1);
5443e65e 2920 continue;
2921 }
2922 fTimeBinIndex[index] = i;
2923 }
2924
2925 Double_t x1, dx1, x2, dx2, gap;
2926
2927 for(Int_t i = 0; i < fN-1; i++) {
2928 x1 = fLayers[i]->GetX();
2929 dx1 = fLayers[i]->GetdX();
2930 x2 = fLayers[i+1]->GetX();
2931 dx2 = fLayers[i+1]->GetdX();
2932 gap = (x2 - dx2/2) - (x1 + dx1/2);
7c1698cb 2933// if(gap < -0.01) {
2934// printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2935// printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2936// }
2937// if(gap > 0.01) {
2938// printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2939// printf(" (%f - %f) - (%f + %f) = %f\n",
2940// x2, dx2/2, x1, dx1, gap);
2941// }
5443e65e 2942 }
2943}
2944
2945
2946//______________________________________________________
2947
2948
2949Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2950{
2951 //
2952 // Returns the number of time bin which in radial position is closest to <x>
2953 //
2954
2955 if(x >= fLayers[fN-1]->GetX()) return fN-1;
2956 if(x <= fLayers[0]->GetX()) return 0;
2957
2958 Int_t b=0, e=fN-1, m=(b+e)/2;
2959 for (; b<e; m=(b+e)/2) {
2960 if (x > fLayers[m]->GetX()) b=m+1;
2961 else e=m;
2962 }
2963 if(TMath::Abs(x - fLayers[m]->GetX()) >
2964 TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2965 else return m;
2966
2967}
2968
2969//______________________________________________________
2970
2971Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2972{
2973 //
2974 // Returns number of the innermost SENSITIVE propagation layer
2975 //
2976
2977 return GetLayerNumber(0);
2978}
2979
2980//______________________________________________________
2981
2982Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2983{
2984 //
2985 // Returns number of the outermost SENSITIVE time bin
2986 //
2987
2988 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 2989}
2990
5443e65e 2991//______________________________________________________
2992
2993Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2994{
2995 //
2996 // Returns number of SENSITIVE time bins
2997 //
2998
2999 Int_t tb, layer;
029cd327 3000 for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
5443e65e 3001 layer = GetLayerNumber(tb);
3002 if(layer>=0) break;
3003 }
3004 return tb+1;
3005}
3006
3007//______________________________________________________
3008
3009void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
3010{
3011 //
3012 // Insert layer <pl> in fLayers array.
3013 // Layers are sorted according to X coordinate.
3014
029cd327 3015 if ( fN == ((Int_t) kMaxLayersPerSector)) {
5443e65e 3016 printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
3017 return;
3018 }
3019 if (fN==0) {fLayers[fN++] = pl; return;}
3020 Int_t i=Find(pl->GetX());
3021
3022 memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3023 fLayers[i]=pl; fN++;
3024
3025}
3026
3027//______________________________________________________
3028
3029Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
3030{
3031 //
3032 // Returns index of the propagation layer nearest to X
3033 //
3034
3035 if (x <= fLayers[0]->GetX()) return 0;
3036 if (x > fLayers[fN-1]->GetX()) return fN;
3037 Int_t b=0, e=fN-1, m=(b+e)/2;
3038 for (; b<e; m=(b+e)/2) {
3039 if (x > fLayers[m]->GetX()) b=m+1;
3040 else e=m;
3041 }
3042 return m;
3043}
3044
7ad19338 3045
3046
3047
3048
3c625a9b 3049//______________________________________________________
3050void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
3051{
3052 //
3053 // set centers and the width of sectors
3054 for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
3055 fZc[icham] = center[icham];
3056 fZmax[icham] = w[icham];
3057 fZmaxSensitive[icham] = wsensitive[icham];
3058 // printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
3059 }
3060}
5443e65e 3061//______________________________________________________
3062
3c625a9b 3063void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3064{
3065 //
3066 // set centers and the width of sectors
3067 fHole = kFALSE;
3068 for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
3069 fIsHole[icham] = holes[icham];
3070 if (holes[icham]) fHole = kTRUE;
3071 }
3072}
3073
3074
3075
7ad19338 3076Bool_t AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
029cd327 3077 Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength,
a9814c09 3078 Bool_t &lookForCluster) const
5443e65e 3079{
3080 //
029cd327 3081 // Returns radial step <dx>, density <rho>, rad. length <radLength>,
5443e65e 3082 // and sensitivity <lookForCluster> in point <y,z>
3083 //
3084
4f1c04d3 3085 Double_t alpha = AliTRDgeometry::GetAlpha();
3086 Double_t ymax = fX*TMath::Tan(0.5*alpha);
3087
3088
5443e65e 3089 dx = fdX;
3090 rho = fRho;
029cd327 3091 radLength = fX0;
5443e65e 3092 lookForCluster = kFALSE;
4f1c04d3 3093 Bool_t cross =kFALSE;
3094 //
3095 //
3096 if ( (ymax-TMath::Abs(y))<3.){ //cross material
3097 rho*=40.;
3098 radLength*=40.;
3099 cross=kTRUE;
3100 }
3c625a9b 3101 //
3102 // check dead regions in sensitive volume
3c625a9b 3103 //
4f1c04d3 3104 Int_t zone=-1;
3105 for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
3106 if (TMath::Abs(z - fZc[ch]) > fZmax[ch]) continue; //not in given zone
3107 //
3108 if (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){
3109 if (fTimeBinIndex>=0) lookForCluster = !(fIsHole[zone]);
3110 if(TMath::Abs(y) > fYmaxSensitive){
3111 lookForCluster = kFALSE;
3112 }
3113 if (fIsHole[zone]) {
3114 //if hole
3115 rho = 1.29e-3;
3116 radLength = 36.66;
3117 }
3118 }else{
7ad19338 3119 cross = kTRUE; rho = 2.7; radLength = 24.01; //aluminium in between
4f1c04d3 3120 }
5443e65e 3121 }
3c625a9b 3122 //
7ad19338 3123 if (fTimeBinIndex>=0) return cross;
4f1c04d3 3124 //
3c625a9b 3125 //
5443e65e 3126 // check hole
7ad19338 3127 if (fHole==kFALSE) return cross;
3c625a9b 3128 //
3129 for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
3130 if (TMath::Abs(z - fZc[ch]) < fZmax[ch]){
3131 if (fIsHole[ch]) {
3132 //if hole
3133 rho = 1.29e-3;
3134 radLength = 36.66;
4f1c04d3 3135 }
3c625a9b 3136 }
3137 }
7ad19338 3138 return cross;
5443e65e 3139}
3140
3c625a9b 3141Int_t AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const
3142{
3143 //
3144 //
3145 if (fTimeBinIndex < 0) return -20; //unknown
3146 Int_t zone=-10; // dead zone
3147 for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
3148 if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
3149 zone = ch;
3150 }
3151 return zone;
3152}
3153
3154
5443e65e 3155//______________________________________________________
3156
3157void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
a9814c09 3158 UInt_t index) {
5443e65e 3159
3160// Insert cluster in cluster array.
3161// Clusters are sorted according to Y coordinate.
3162
3163 if(fTimeBinIndex < 0) {
3164 printf("*** attempt to insert cluster into non-sensitive time bin!\n");
3165 return;
3166 }
3167
029cd327 3168 if (fN== (Int_t) kMaxClusterPerTimeBin) {
5443e65e 3169 printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
3170 return;
3171 }
3172 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
3173 Int_t i=Find(c->GetY());
3174 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3175 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
3176 fIndex[i]=index; fClusters[i]=c; fN++;
3177}
3178
3179//______________________________________________________
3180
69b55c55 3181Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const {
5443e65e 3182
3183// Returns index of the cluster nearest in Y
3184
69b55c55 3185 if (fN<=0) return 0;
5443e65e 3186 if (y <= fClusters[0]->GetY()) return 0;
3187 if (y > fClusters[fN-1]->GetY()) return fN;
3188 Int_t b=0, e=fN-1, m=(b+e)/2;
3189 for (; b<e; m=(b+e)/2) {
3190 if (y > fClusters[m]->GetY()) b=m+1;
3191 else e=m;
3192 }
3193 return m;
3194}
3195
69b55c55 3196Int_t AliTRDtracker::AliTRDpropagationLayer::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const
7ad19338 3197{
3198 //
3199 // Returns index of the cluster nearest to the given y,z
3200 //
3201 Int_t index = -1;
3202 Int_t maxn = fN;
69b55c55 3203 Float_t mindist = maxroad;
7ad19338 3204 //
3205 for (Int_t i=Find(y-maxroad); i<maxn; i++) {
3206 AliTRDcluster* c=(AliTRDcluster*)(fClusters[i]);
69b55c55 3207 Float_t ycl = c->GetY();
7ad19338 3208 //
69b55c55 3209 if (ycl > y+maxroad) break;
3210 if (TMath::Abs(c->GetZ()-z) > maxroadz) continue;
3211 if (TMath::Abs(ycl-y)<mindist){
3212 mindist = TMath::Abs(ycl-y);
3213 index = fIndex[i];
7ad19338 3214 }
3215 }
3216 return index;
3217}
3218
3219
fd621f36 3220//---------------------------------------------------------
5443e65e 3221
fd621f36 3222Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
3223//
3224// Returns correction factor for tilted pads geometry
3225//
fd621f36 3226 Int_t det = c->GetDetector();
3227 Int_t plane = fGeom->GetPlane(det);
3551db50 3228 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
de4b10e5 3229 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 3230
3231 if(fNoTilt) h01 = 0;
fd621f36 3232 return h01;
3233}
5443e65e 3234
c630aafd 3235
eab5961e 3236void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
3237{
3238 // *** ADDED TO GET MORE INFORMATION FOR TRD PID ---- PS
3239 // This is setting fdEdxPlane and fTimBinPlane
3240 // Sums up the charge in each plane for track TRDtrack and also get the
3241 // Time bin for Max. Cluster
3242 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
3243
eab5961e 3244 Double_t clscharge[kNPlane], maxclscharge[kNPlane];
3245 Int_t nCluster[kNPlane], timebin[kNPlane];
3246
3247 //Initialization of cluster charge per plane.
3248 for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3249 clscharge[iPlane] = 0.0;
3250 nCluster[iPlane] = 0;
3251 timebin[iPlane] = -1;
3252 maxclscharge[iPlane] = 0.0;
3253 }
3254
3255 // Loop through all clusters associated to track TRDtrack
3256 Int_t nClus = TRDtrack.GetNumberOfClusters(); // from Kalmantrack
3257 for (Int_t iClus = 0; iClus < nClus; iClus++) {
3258 Double_t charge = TRDtrack.GetClusterdQdl(iClus);
3259 Int_t index = TRDtrack.GetClusterIndex(iClus);
3260 AliTRDcluster *TRDcluster = (AliTRDcluster *) GetCluster(index);
3261 if (!TRDcluster) continue;
3262 Int_t tb = TRDcluster->GetLocalTimeBin();
3263 if (!tb) continue;
3264 Int_t detector = TRDcluster->GetDetector();
3265 Int_t iPlane = fGeom->GetPlane(detector);
3266 clscharge[iPlane] = clscharge[iPlane]+charge;
3267 if(charge > maxclscharge[iPlane]) {
3268 maxclscharge[iPlane] = charge;
3269 timebin[iPlane] = tb;
3270 }
3271 nCluster[iPlane]++;
3272 } // end of loop over cluster
3273
3274 // Setting the fdEdxPlane and fTimBinPlane variabales
3275 Double_t Total_ch = 0;
3276 for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
bd50219c 3277 // Quality control of TRD track.
3278 if (nCluster[iPlane]<= 5) {
3279 clscharge[iPlane]=0.0;
3280 timebin[iPlane]=-1;
3281 }
eab5961e 3282 if (nCluster[iPlane]) clscharge[iPlane] /= nCluster[iPlane];
3283 TRDtrack.SetPIDsignals(clscharge[iPlane], iPlane);
3284 TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
3285 Total_ch= Total_ch+clscharge[iPlane];
3286 }
3287 // Int_t i;
3288 // Int_t nc=TRDtrack.GetNumberOfClusters();
3289 // Float_t dedx=0;
3290 // for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
3291 // dedx /= nc;
3292 // for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3293 // TRDtrack.SetPIDsignals(dedx, iPlane);
3294 // TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
3295 // }
3296
3297} // end of function
3298
c630aafd 3299
7ad19338 3300Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters,AliTRDtracklet&tracklet)
4f1c04d3 3301{
3302 //
3303 //
3304 // try to find nearest clusters to the track in timebins from t0 to t1
3305 //
3306 //
7ad19338 3307 //
3308 // correction coeficients - depends on TRD parameters - to be changed according it
3309 //
3310
3311 Double_t x[100],yt[100],zt[100];
3312 Double_t xmean=0; //reference x
3313 Double_t dz[10][100],dy[10][100];
3314 Float_t zmean[100], nmean[100];
3315 Int_t clfound=0;
3316 Int_t indexes[10][100]; // indexes of the clusters in the road
3317 AliTRDcluster *cl[10][100]; // pointers to the clusters in the road
3318 Int_t best[10][100]; // index of best matching cluster
3319 //
3320 //
69b55c55 3321
8979685e 3322 for (Int_t it=0;it<=t1-t0; it++){
4f1c04d3 3323 x[it]=0;
3324 yt[it]=0;
3325 zt[it]=0;
7ad19338 3326 clusters[it+t0]=-2;
3327 zmean[it]=0;
3328 nmean[it]=0;
3329 //
3330 for (Int_t ih=0;ih<10;ih++){
3331 indexes[ih][it]=-2; //reset indexes1
3332 cl[ih][it]=0;
3333 dz[ih][it]=-100;
3334 dy[ih][it]=-100;
3335 best[ih][it]=0;
3336 }
4f1c04d3 3337 }
3338 //
3339 Double_t x0 = track->GetX();
69b55c55 3340 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
4f1c04d3 3341 Int_t nall=0;
3342 Int_t nfound=0;
7ad19338 3343 Double_t h01 =0;
3344 Int_t plane =-1;
3345 Float_t padlength=0;
3346 AliTRDtrack track2(*track);
3347 Float_t snpy = track->GetSnp();
3348 Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy));
3349 if (snpy<0) tany*=-1;
3350 //
3351 Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
3352 Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
3353 Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2);
3354 if (road>6.) road=6.;
4f1c04d3 3355
7ad19338 3356 //
3357 for (Int_t it=0;it<t1-t0;it++){
3358 Double_t maxChi2[2]={fgkMaxChi2,fgkMaxChi2};
3359 AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it+t0));
4f1c04d3 3360 if (timeBin==0) continue; // no indexes1
3361 Int_t maxn = timeBin;
3362 x[it] = timeBin.GetX();
7ad19338 3363 track2.PropagateTo(x[it]);
3364 yt[it] = track2.GetY();
3365 zt[it] = track2.GetZ();
3366
3367 Double_t y=yt[it],z=zt[it];
4f1c04d3 3368 Double_t chi2 =1000000;
3369 nall++;
3370 //
7ad19338 3371 // find 2 nearest cluster at given time bin
3372 //
3373 //
4f1c04d3 3374 for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
3375 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
7ad19338 3376 h01 = GetTiltFactor(c);
3377 if (plane<0){
3378 Int_t det = c->GetDetector();
3379 plane = fGeom->GetPlane(det);
3380 padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
3381 }
3382 // if (c->GetLocalTimeBin()==0) continue;
4f1c04d3 3383 if (c->GetY() > y+road) break;
7ad19338 3384 if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;
3385
3386 Double_t dist = TMath::Abs(c->GetZ()-z);
3387 if (dist> (0.5*padlength+6.*sigmaz)) continue; // 6 sigma boundary cut
3388 Double_t cost = 0;
3389 //
3390 if (dist> (0.5*padlength-sigmaz)){ // sigma boundary cost function
3391 cost = (dist-0.5*padlength)/(2.*sigmaz);
3392 if (cost>-1) cost= (cost+1.)*(cost+1.);
3393 else cost=0;
3394 }
3395 // Int_t label = TMath::Abs(track->GetLabel());
3396 // if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3397 chi2=track2.GetPredictedChi2(c,h01)+cost;
3398 //
3399 clfound++;
3400 if (chi2 > maxChi2[1]) continue;
3401
3402 for (Int_t ih=2;ih<9; ih++){ //store the clusters in the road
3403 if (cl[ih][it]==0){
3404 cl[ih][it] = c;
3405 indexes[ih][it] =timeBin.GetIndex(i); // index - 9 - reserved for outliers
3406 break;
3407 }
4f1c04d3 3408 }
7ad19338 3409 //
3410 if (chi2 <maxChi2[0]){
3411 maxChi2[1] = maxChi2[0];
3412 maxChi2[0] = chi2;
3413 indexes[1][it] = indexes[0][it];
3414 cl[1][it] = cl[0][it];
3415 indexes[0][it] = timeBin.GetIndex(i);
3416 cl[0][it] = c;
3417 continue;
3418 }
3419 maxChi2[1]=chi2;
3420 cl[1][it] = c;
3421 indexes[1][it] =timeBin.GetIndex(i);
3422 }
3423 if (cl[0][it]){
3424 nfound++;
3425 xmean += x[it];
3426 }
4f1c04d3 3427 }
3428 //
7ad19338 3429 if (nfound<4) return 0;
3430 xmean /=Float_t(nfound); // middle x
3431 track2.PropagateTo(xmean); // propagate track to the center
4f1c04d3 3432 //
3433 // choose one of the variants
3434 //
7ad19338 3435 Int_t changes[10];
3436 Float_t sumz = 0;
3437 Float_t sum = 0;
3438 Double_t sumdy = 0;
3439 Double_t sumdy2 = 0;
3440 Double_t sumx = 0;
3441 Double_t sumxy = 0;
3442 Double_t sumx2 = 0;
3443 Double_t mpads = 0;
3444 //
3445 Int_t ngood[10];
3446 Int_t nbad[10];
3447 //
3448 Double_t meanz[10];
3449 Double_t moffset[10]; // mean offset
3450 Double_t mean[10]; // mean value
3451 Double_t angle[10]; // angle
3452 //
3453 Double_t smoffset[10]; // sigma of mean offset
3454 Double_t smean[10]; // sigma of mean value
3455 Double_t sangle[10]; // sigma of angle
3456 Double_t smeanangle[10]; // correlation
3457 //
3458 Double_t sigmas[10];
3459 Double_t tchi2s[10]; // chi2s for tracklet
3460 //
3461 // calculate zmean
3462 //
3463 for (Int_t it=0;it<t1-t0;it++){
3464 if (!cl[0][it]) continue;
3465 for (Int_t dt=-3;dt<=3;dt++){
3466 if (it+dt<0) continue;
8979685e 3467 if (it+dt>t1-t0) continue;
7ad19338 3468 if (!cl[0][it+dt]) continue;
3469 zmean[it]+=cl[0][it+dt]->GetZ();
3470 nmean[it]+=1.;
3471 }
3472 zmean[it]/=nmean[it];
3473 }
3474 //
3475 for (Int_t it=0; it<t1-t0;it++){
3476 best[0][it]=0;
3477 for (Int_t ih=0;ih<10;ih++){
3478 dz[ih][it]=-100;
3479 dy[ih][it]=-100;
4f1c04d3 3480 if (!cl[ih][it]) continue;
69b55c55 3481 //Float_t poscor = fgkCoef*(cl[ih][it]->GetLocalTimeBin() - fgkMean)+fgkOffset;
3482 Float_t poscor = 0; // applied during loading of clusters
3483 if (cl[ih][it]->IsUsed()) poscor=0; // correction already applied
7ad19338 3484 dz[ih][it] = cl[ih][it]->GetZ()- zt[it]; // calculate distance from track in z
3485 dy[ih][it] = cl[ih][it]->GetY()+ dz[ih][it]*h01 - poscor -yt[it]; // in y
3486 }
3487 // minimize changes
3488 if (!cl[0][it]) continue;
3489 if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it])
3490 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){
3491 best[0][it]=1;
4f1c04d3 3492 }
7ad19338 3493 }
3494 //
3495 // iterative choosing of "best path"
3496 //
3497 //
3498 Int_t label = TMath::Abs(track->GetLabel());
3499 Int_t bestiter=0;
3500 //
3501 for (Int_t iter=0;iter<9;iter++){
3502 //
3503 changes[iter]= 0;
3504 sumz = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0;
3505 // linear fit
3506 for (Int_t it=0;it<t1-t0;it++){
3507 if (!cl[best[iter][it]][it]) continue;
3508 //calculates pad-row changes
3509 Double_t zbefore= cl[best[iter][it]][it]->GetZ();
3510 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3511 for (Int_t itd = it-1; itd>=0;itd--) {
3512 if (cl[best[iter][itd]][itd]) {
3513 zbefore= cl[best[iter][itd]][itd]->GetZ();
3514 break;
3515 }
3516 }
3517 for (Int_t itd = it+1; itd<t1-t0;itd++) {
3518 if (cl[best[iter][itd]][itd]) {
3519 zafter= cl[best[iter][itd]][itd]->GetZ();
3520 break;
3521 }
3522 }
3523 if (TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore)>0.1&&TMath::Abs(cl[best[iter][it]][it]->GetZ()-zafter)>0.1) changes[iter]++;
3524 //
3525 Double_t dx = x[it]-xmean; // distance to reference x
3526 sumz += cl[best[iter][it]][it]->GetZ();
4f1c04d3 3527 sum++;
7ad19338 3528 sumdy += dy[best[iter][it]][it];
3529 sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it];
3530 sumx += dx;
3531 sumx2 += dx*dx;
3532 sumxy += dx*dy[best[iter][it]][it];
3533 mpads += cl[best[iter][it]][it]->GetNPads();
3534 if (cl[best[iter][it]][it]->GetLabel(0)==label || cl[best[iter][it]][it]->GetLabel(1)==label||cl[best[iter][it]][it]->GetLabel(2)==label){
3535 ngood[iter]++;
4f1c04d3 3536 }
3537 else{
7ad19338 3538 nbad[iter]++;
4f1c04d3 3539 }
3540 }
7ad19338 3541 //
3542 // calculates line parameters
3543 //
3544 Double_t det = sum*sumx2-sumx*sumx;
3545 angle[iter] = (sum*sumxy-sumx*sumdy)/det;
3546 mean[iter] = (sumx2*sumdy-sumx*sumxy)/det;
3547 meanz[iter] = sumz/sum;
3548 moffset[iter] = sumdy/sum;
3549 mpads /= sum; // mean number of pads
3550 //
3551 //
3552 Double_t sigma2 = 0; // normalized residuals - for line fit
3553 Double_t sigma1 = 0; // normalized residuals - constant fit
3554 //
3555 for (Int_t it=0;it<t1-t0;it++){
3556 if (!cl[best[iter][it]][it]) continue;
3557 Double_t dx = x[it]-xmean;
3558 Double_t ytr = mean[iter]+angle[iter]*dx;
3559 sigma2 += (dy[best[iter][it]][it]-ytr)*(dy[best[iter][it]][it]-ytr);
3560 sigma1 += (dy[best[iter][it]][it]-moffset[iter])*(dy[best[iter][it]][it]-moffset[iter]);
3561 sum++;
4f1c04d3 3562 }
7ad19338 3563 sigma2 /=(sum-2); // normalized residuals
3564 sigma1 /=(sum-1); // normalized residuals
3565 //
3566 smean[iter] = sigma2*(sumx2/det); // estimated error2 of mean
3567 sangle[iter] = sigma2*(sum/det); // estimated error2 of angle
3568 smeanangle[iter] = sigma2*(-sumx/det); // correlation
3569 //
3570 //
3571 sigmas[iter] = TMath::Sqrt(sigma1); //
3572 smoffset[iter]= (sigma1/sum)+0.01*0.01; // sigma of mean offset + unisochronity sigma
3573 //
3574 // iterative choosing of "better path"
3575 //
3576 for (Int_t it=0;it<t1-t0;it++){
3577 if (!cl[best[iter][it]][it]) continue;
3578 //
3579 Double_t sigmatr2 = smoffset[iter]+0.5*tany*tany; //add unisochronity + angular effect contribution
69b55c55 3580 Double_t sweight = 1./sigmatr2+1./track->GetSigmaY2();
7ad19338 3581 Double_t weighty = (moffset[iter]/sigmatr2)/sweight; // weighted mean
69b55c55 3582 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1+track->GetSigmaY2()); //
7ad19338 3583 Double_t mindist=100000;
3584 Int_t ihbest=0;
3585 for (Int_t ih=0;ih<10;ih++){
3586 if (!cl[ih][it]) break;
3587 Double_t dist2 = (dy[ih][it]-weighty)/sigmacl;
3588 dist2*=dist2; //chi2 distance
3589 if (dist2<mindist){
3590 mindist = dist2;
3591 ihbest =ih;
3592 }
3593 }
3594 best[iter+1][it]=ihbest;
4f1c04d3 3595 }
4f1c04d3 3596 //
7ad19338 3597 // update best hypothesy if better chi2 according tracklet position and angle
3598 //
69b55c55 3599 Double_t sy2 = smean[iter] + track->GetSigmaY2();
7ad19338 3600 Double_t sa2 = sangle[iter] + track->fCee;
3601 Double_t say = track->fCey;
3602 // Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
3603 // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
3604
3605 Double_t detchi = sy2*sa2-say*say;
3606 Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi}; //inverse value of covariance matrix
4f1c04d3 3607
7ad19338 3608 Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+
3609 2.*mean[bestiter]*angle[bestiter]*invers[2];
3610 Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+
3611 2*mean[iter]*angle[iter]*invers[2];
3612 tchi2s[iter] =chi21;
3613 //
3614 if (changes[iter]<=changes[bestiter] && chi21<chi20) {
3615 bestiter =iter;
3616 }
3617 }
3618 //
3619 //set clusters
3620 //
3621 Double_t sigma2 = sigmas[0]; // choose as sigma from 0 iteration
8979685e 3622 Short_t maxpos = -1;
3623 Float_t maxcharge = 0;
3624 Short_t maxpos4 = -1;
3625 Float_t maxcharge4 = 0;
3626 Short_t maxpos5 = -1;
3627 Float_t maxcharge5 = 0;
3628
7ad19338 3629 //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3630 //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
3631
3632 Double_t expectederr = sigma2*sigma2+0.01*0.01;
3633 if (mpads>3.5) expectederr += (mpads-3.5)*0.04;
3634 if (changes[bestiter]>1) expectederr+= changes[bestiter]*0.01;
3635 expectederr+=(0.03*(tany-fgkExB)*(tany-fgkExB))*15;
3636 // if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
3637 //expectederr+=10000;
3638 for (Int_t it=0;it<t1-t0;it++){
3639 if (!cl[best[bestiter][it]][it]) continue;
69b55c55 3640 // Float_t poscor = fgkCoef*(cl[best[bestiter][it]][it]->GetLocalTimeBin() - fgkMean)+fgkOffset;
3641 Float_t poscor = 0; //applied during loading of cluster
7ad19338 3642 cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // set cluster error
3643 if (!cl[best[bestiter][it]][it]->IsUsed()){
3644 cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY()-poscor); // ExB corrction correction
69b55c55 3645 // cl[best[bestiter][it]][it]->Use();
3646 }
3647 //
3648 // time bins with maximal charge
3649 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
3650 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3651 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3652 }
3653
3654 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
3655 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
3656 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3657 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3658 }
3659 }
3660 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
3661 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
3662 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3663 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3664 }
7ad19338 3665 }
8979685e 3666 //
3667 // time bins with maximal charge
3668 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
3669 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3670 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3671 }
3672
3673 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
3674 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
3675 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3676 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3677 }
3678 }
3679 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
3680 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
3681 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3682 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3683 }
3684 }
7ad19338 3685 clusters[it+t0] = indexes[best[bestiter][it]][it];
69b55c55 3686 //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 && cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0] = indexes[best[bestiter][it]][it]; //Test
7ad19338 3687 }
3688 //
3689 // set tracklet parameters
3690 //
3691 Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01;
3692 if (mpads>3.5) trackleterr2 += (mpads-3.5)*0.04;
3693 trackleterr2+= changes[bestiter]*0.01;
3694 trackleterr2*= TMath::Max(14.-nfound,1.);
3695 trackleterr2+= 0.2*(tany-fgkExB)*(tany-fgkExB);
3696 //
3697 tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2); //set tracklet parameters
3698 tracklet.SetTilt(h01);
3699 tracklet.SetP0(mean[bestiter]);
3700 tracklet.SetP1(angle[bestiter]);
3701 tracklet.SetN(nfound);
3702 tracklet.SetNCross(changes[bestiter]);
3703 tracklet.SetPlane(plane);
3704 tracklet.SetSigma2(expectederr);
3705 tracklet.SetChi2(tchi2s[bestiter]);
8979685e 3706 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
7ad19338 3707 track->fTracklets[plane] = tracklet;
3708 track->fNWrong+=nbad[0];
3709 //
3710 // Debuging part
3711 //
69b55c55 3712 TClonesArray array0("AliTRDcluster");
3713 TClonesArray array1("AliTRDcluster");
3714 array0.ExpandCreateFast(t1-t0+1);
3715 array1.ExpandCreateFast(t1-t0+1);
7ad19338 3716 TTreeSRedirector& cstream = *fDebugStreamer;
3717 AliTRDcluster dummy;
3718 Double_t dy0[100];
8979685e 3719 Double_t dyb[100];
3720
7ad19338 3721 for (Int_t it=0;it<t1-t0;it++){
3722 dy0[it] = dy[0][it];
3723 dyb[it] = dy[best[bestiter][it]][it];
3724 if(cl[0][it]) {
3725 new(array0[it]) AliTRDcluster(*cl[0][it]);
3726 }
3727 else{
3728 new(array0[it]) AliTRDcluster(dummy);
3729 }
3730 if(cl[best[bestiter][it]][it]) {
3731 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3732 }
3733 else{
3734 new(array1[it]) AliTRDcluster(dummy);
3735 }
4f1c04d3 3736 }
7ad19338 3737 TGraph graph0(t1-t0,x,dy0);
3738 TGraph graph1(t1-t0,x,dyb);
3739 TGraph graphy(t1-t0,x,yt);
3740 TGraph graphz(t1-t0,x,zt);
3741 //
3742 //
3743 cstream<<"tracklet"<<
3744 "track.="<<track<< // track parameters
3745 "tany="<<tany<< // tangent of the local track angle
3746 "xmean="<<xmean<< // xmean - reference x of tracklet
3747 "tilt="<<h01<< // tilt angle
3748 "nall="<<nall<< // number of foundable clusters
3749 "nfound="<<nfound<< // number of found clusters
3750 "clfound="<<clfound<< // total number of found clusters in road
3751 "mpads="<<mpads<< // mean number of pads per cluster
3752 "plane="<<plane<< // plane number
3753 "road="<<road<< // the width of the used road
3754 "graph0.="<<&graph0<< // x - y = dy for closest cluster
3755 "graph1.="<<&graph1<< // x - y = dy for second closest cluster
3756 "graphy.="<<&graphy<< // y position of the track
3757 "graphz.="<<&graphz<< // z position of the track
69b55c55 3758 // "fCl.="<<&array0<< // closest cluster
f6625211 3759 //"fCl2.="<<&array1<< // second closest cluster
8979685e 3760 "maxpos="<<maxpos<< // maximal charge postion
3761 "maxcharge="<<maxcharge<< // maximal charge
3762 "maxpos4="<<maxpos4<< // maximal charge postion - after bin 4
3763 "maxcharge4="<<maxcharge4<< // maximal charge - after bin 4
3764 "maxpos5="<<maxpos5<< // maximal charge postion - after bin 5
3765 "maxcharge5="<<maxcharge5<< // maximal charge - after bin 5
7ad19338 3766 //
3767 "bestiter="<<bestiter<< // best iteration number
3768 "tracklet.="<<&tracklet<< // corrspond to the best iteration
3769 "tchi20="<<tchi2s[0]<< // chi2 of cluster in the 0 iteration
3770 "tchi2b="<<tchi2s[bestiter]<< // chi2 of cluster in the best iteration
3771 "sigmas0="<<sigmas[0]<< // residuals sigma
3772 "sigmasb="<<sigmas[bestiter]<< // residulas sigma
3773 //
3774 "ngood0="<<ngood[0]<< // number of good clusters in 0 iteration
3775 "nbad0="<<nbad[0]<< // number of bad clusters in 0 iteration
3776 "ngoodb="<<ngood[bestiter]<< // in best iteration
3777 "nbadb="<<nbad[bestiter]<< // in best iteration
3778 //
3779 "changes0="<<changes[0]<< // changes of pardrows in iteration number 0
3780 "changesb="<<changes[bestiter]<< // changes of pardrows in best iteration
3781 //
3782 "moffset0="<<moffset[0]<< // offset fixing angle in iter=0
3783 "smoffset0="<<smoffset[0]<< // sigma of offset fixing angle in iter=0
3784 "moffsetb="<<moffset[bestiter]<< // offset fixing angle in iter=best
3785 "smoffsetb="<<smoffset[bestiter]<< // sigma of offset fixing angle in iter=best
3786 //
3787 "mean0="<<mean[0]<< // mean dy in iter=0;
3788 "smean0="<<smean[0]<< // sigma of mean dy in iter=0
3789 "meanb="<<mean[bestiter]<< // mean dy in iter=best
3790 "smeanb="<<smean[bestiter]<< // sigma of mean dy in iter=best
3791 //
3792 "angle0="<<angle[0]<< // angle deviation in the iteration number 0
3793 "sangle0="<<sangle[0]<< // sigma of angular deviation in iteration number 0
3794 "angleb="<<angle[bestiter]<< // angle deviation in the best iteration
3795 "sangleb="<<sangle[bestiter]<< // sigma of angle deviation in the best iteration
3796 //
3797 "expectederr="<<expectederr<< // expected error of cluster position
3798 "\n";
3799 //
3800 //
4f1c04d3 3801 return nfound;
3802}
3803
3804
69b55c55 3805Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down)
3806{
3807 //
3808 // Sort eleements according occurancy
3809 // The size of output array has is 2*n
3810 //
3811 Int_t * sindexS = new Int_t[n]; // temp array for sorting
3812 Int_t * sindexF = new Int_t[2*n];
3813 for (Int_t i=0;i<n;i++) sindexF[i]=0;
3814 //
3815 TMath::Sort(n,inlist, sindexS, down);
3816 Int_t last = inlist[sindexS[0]];
3817 Int_t val = last;
3818 sindexF[0] = 1;
3819 sindexF[0+n] = last;
3820 Int_t countPos = 0;
3821 //
3822 // find frequency
3823 for(Int_t i=1;i<n; i++){
3824 val = inlist[sindexS[i]];
3825 if (last == val) sindexF[countPos]++;
3826 else{
3827 countPos++;
3828 sindexF[countPos+n] = val;
3829 sindexF[countPos]++;
3830 last =val;
3831 }
3832 }
3833 if (last==val) countPos++;
3834 // sort according frequency
3835 TMath::Sort(countPos, sindexF, sindexS, kTRUE);
3836 for (Int_t i=0;i<countPos;i++){
3837 outlist[2*i ] = sindexF[sindexS[i]+n];
3838 outlist[2*i+1] = sindexF[sindexS[i]];
3839 }
3840 delete [] sindexS;
3841 delete [] sindexF;
3842
3843 return countPos;
3844}
3845
3846AliTRDtrack * AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
3847{
3848 //
3849 //
3850 //
3851 Double_t alpha=AliTRDgeometry::GetAlpha();
3852 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
3853 Double_t c[15];
3854 c[0] = 0.2;
3855 c[1] = 0 ; c[2] = 2;
3856 c[3] = 0 ; c[4] = 0; c[5] = 0.02;
3857 c[6] = 0 ; c[7] = 0; c[8] = 0; c[9] = 0.1;
3858 c[10] = 0 ; c[11] = 0; c[12] = 0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
3859 //
3860 Int_t index =0;
3861 AliTRDcluster *cl =0;
3862 for (Int_t ilayer=0;ilayer<6;ilayer++){
3863 if (seeds[ilayer].isOK()){
3864 for (Int_t itime=22;itime>0;itime--){
3865 if (seeds[ilayer].fIndexes[itime]>0){
3866 index = seeds[ilayer].fIndexes[itime];
3867 cl = seeds[ilayer].fClusters[itime];
3868 break;
3869 }
3870 }
3871 }
3872 if (index>0) break;
3873 }
3874 if (cl==0) return 0;
3875 AliTRDtrack * track = new AliTRDtrack(cl,index,&params[1],c, params[0],params[6]*alpha+shift);
3876 track->PropagateTo(params[0]-5.);
3877 track->ResetCovariance(1);
3878 //
f6625211 3879 Int_t rc=FollowBackProlongation(*track);
69b55c55 3880 if (rc<30) {
3881 delete track;
3882 track =0;
3883 }else{
3884 track->CookdEdx();
3885 CookdEdxTimBin(*track);
3886 CookLabel(track, 0.9);
3887 }
3888 return track;
3889}
3890
3891
3892
3893
3894
3895
3896AliTRDseed::AliTRDseed()
3897{
3898 //
3899 //
3900 fTilt =0; // tilting angle
3901 fPadLength = 0; // pad length
3902 fX0 = 0; // x0 position
3903 for (Int_t i=0;i<25;i++){
3904 fX[i]=0; // !x position
3905 fY[i]=0; // !y position
3906 fZ[i]=0; // !z position
3907 fIndexes[i]=0; // !indexes
3908 fClusters[i]=0; // !clusters
3909 }
3910 for (Int_t i=0;i<2;i++){
3911 fYref[i]=0; // reference y
3912 fZref[i]=0; // reference z
3913 fYfit[i]=0; // y fit position +derivation
3914 fYfitR[i]=0; // y fit position +derivation
3915 fZfit[i]=0; // z fit position
3916 fZfitR[i]=0; // z fit position
3917 fLabels[i]=0; // labels
3918 }
3919 fSigmaY = 0;
3920 fSigmaY2 = 0;
3921 fMeanz=0; // mean vaue of z
3922 fZProb=0; // max probbable z
3923 fMPads=0;
3924 //
3925 fN=0; // number of associated clusters
3926 fN2=0; // number of not crossed
3927 fNUsed=0; // number of used clusters
3928 fNChange=0; // change z counter
3929}
3930
3931void AliTRDseed::Reset(){
3932 //
3933 // reset seed
3934 //
3935 for (Int_t i=0;i<25;i++){
3936 fX[i]=0; // !x position
3937 fY[i]=0; // !y position
3938 fZ[i]=0; // !z position
3939 fIndexes[i]=0; // !indexes
3940 fClusters[i]=0; // !clusters
3941 fUsable[i] = kFALSE;
3942 }
3943 for (Int_t i=0;i<2;i++){
3944 fYref[i]=0; // reference y
3945 fZref[i]=0; // reference z
3946 fYfit[i]=0; // y fit position +derivation
3947 fYfitR[i]=0; // y fit position +derivation
3948 fZfit[i]=0; // z fit position
3949 fZfitR[i]=0; // z fit position
3950 fLabels[i]=-1; // labels
3951 }
3952 fSigmaY =0; //"robust" sigma in y
3953 fSigmaY2=0; //"robust" sigma in y
3954 fMeanz =0; // mean vaue of z
3955 fZProb =0; // max probbable z
3956 fMPads =0;
3957 //
3958 fN=0; // number of associated clusters
3959 fN2=0; // number of not crossed
3960 fNUsed=0; // number of used clusters
3961 fNChange=0; // change z counter
3962}
3963
3964void AliTRDseed::CookLabels(){
3965 //
3966 // cook 2 labels for seed
3967 //
3968 Int_t labels[200];
3969 Int_t out[200];
3970 Int_t nlab =0;
3971 for (Int_t i=0;i<25;i++){
3972 if (!fClusters[i]) continue;
3973 for (Int_t ilab=0;ilab<3;ilab++){
3974 if (fClusters[i]->GetLabel(ilab)>=0){
3975 labels[nlab] = fClusters[i]->GetLabel(ilab);
3976 nlab++;
3977 }
3978 }
3979 }
3980 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
3981 fLabels[0] = out[0];
3982 if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
3983}
3984
3985void AliTRDseed::UseClusters()
3986{
3987 //
3988 // use clusters
3989 //
3990 for (Int_t i=0;i<25;i++){
3991 if (!fClusters[i]) continue;
3992 if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
3993 }
3994}
3995
3996
3997void AliTRDseed::Update(){
3998 //
3999 //
4000 //
4001 const Float_t ratio = 0.8;
4002 const Int_t kClmin = 6;
4003 const Float_t kmaxtan = 2;
4004 if (TMath::Abs(fYref[1])>kmaxtan) return; // too much inclined track
4005 //
4006 Float_t sigmaexp = 0.05+TMath::Abs(fYref[1]*0.25); // expected r.m.s in y direction
4007 Float_t ycrosscor = fPadLength*fTilt*0.5; // y correction for crossing
4008 fNChange =0;
4009 //
4010 Double_t sumw, sumwx,sumwx2;
4011 Double_t sumwy, sumwxy, sumwz,sumwxz;
4012 Int_t zints[25]; // histograming of the z coordinate - get 1 and second max probable coodinates in z
4013 Int_t zouts[50]; //
4014 Float_t allowedz[25]; // allowed z for given time bin
4015 Float_t yres[25]; // residuals from reference
4016 Float_t anglecor = fTilt*fZref[1]; //correction to the angle
4017 //
4018 //
4019 fN=0; fN2 =0;
4020 for (Int_t i=0;i<25;i++){
4021 yres[i] =10000;
4022 if (!fClusters[i]) continue;
4023 yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
4024 zints[fN] = Int_t(fZ[i]);
4025 fN++;
4026 }
4027 if (fN<kClmin) return;
4028 Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
4029 fZProb = zouts[0];
4030 if (nz<=1) zouts[3]=0;
4031 if (zouts[1]+zouts[3]<kClmin) return;
4032 //
4033 if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0; // z distance bigger than pad - length
4034 //
4035 Int_t breaktime = -1;
4036 Bool_t mbefore = kFALSE;
4037 Int_t cumul[25][2];
4038 Int_t counts[2]={0,0};
4039 //
4040 if (zouts[3]>=3){
4041 //
4042 // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
4043 //
4044 fNChange=1;
4045 for (Int_t i=0;i<25;i++){
4046 cumul[i][0] = counts[0];
4047 cumul[i][1] = counts[1];
4048 if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
4049 if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
4050 }
4051 Int_t maxcount = 0;
4052 for (Int_t i=0;i<24;i++) {
4053 Int_t after = cumul[24][0]-cumul[i][0];
4054 Int_t before = cumul[i][1];
4055 if (after+before>maxcount) {
4056 maxcount=after+before;
4057 breaktime=i;
4058 mbefore=kFALSE;
4059 }
4060 after = cumul[24][1]-cumul[i][1];
4061 before = cumul[i][0];
4062 if (after+before>maxcount) {
4063 maxcount=after+before;
4064 breaktime=i;
4065 mbefore=kTRUE;
4066 }
4067 }
4068 breaktime-=1;
4069 }
4070 for (Int_t i=0;i<25;i++){
4071 if (i>breaktime) allowedz[i] = mbefore ? zouts[2]:zouts[0];
4072 if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
4073 }
4074 if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] && fZref[1]>0)){
4075 //
4076 // tracklet z-direction not in correspondance with track z direction
4077 //
4078 fNChange =0;
4079 for (Int_t i=0;i<25;i++){
4080 allowedz[i] = zouts[0]; //only longest taken
4081 }
4082 }
4083 //
4084 if (fNChange>0){
4085 //
4086 // cross pad -row tracklet - take the step change into account
4087 //
4088 for (Int_t i=0;i<25;i++){
4089 if (!fClusters[i]) continue;
4090 if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
4091 yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
4092 if (TMath::Abs(fZ[i]-fZProb)>2){
4093 if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
4094 if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
4095 }
4096 }
4097 }
4098 //
4099 Double_t yres2[25];
4100 Double_t mean,sigma;
4101 for (Int_t i=0;i<25;i++){
4102 if (!fClusters[i]) continue;
4103 if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
4104 yres2[fN2] = yres[i];
4105 fN2++;
4106 }
4107 if (fN2<kClmin){
4108 fN2 = 0;
4109 return;
4110 }
4111 EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*ratio-2));
4112 if (sigma<sigmaexp*0.8) sigma=sigmaexp;
4113 fSigmaY = sigma;
4114 //
4115 //
4116 // reset sums
4117 sumw=0; sumwx=0; sumwx2=0;
4118 sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
4119 fN2 =0;
4120 fMeanz =0;
4121 fMPads =0;
4122 //
4123 for (Int_t i=0;i<25;i++){
4124 fUsable[i]=kFALSE;
4125 if (!fClusters[i]) continue;
4126 if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
4127 if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
4128 fUsable[i] = kTRUE;
4129 fN2++;
4130 fMPads+=fClusters[i]->GetNPads();
4131 Float_t weight =1;
4132 if (fClusters[i]->GetNPads()>4) weight=0.5;
4133 if (fClusters[i]->GetNPads()>5) weight=0.2;
4134 //
4135 Double_t x = fX[i];
4136 sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
4137 sumwy+=weight*yres[i]; sumwxy+=weight*(yres[i])*x;
4138 sumwz+=weight*fZ[i]; sumwxz+=weight*fZ[i]*x;
4139 }
4140 if (fN2<kClmin){
4141 fN2 = 0;
4142 return;
4143 }
4144 fMeanz = sumwz/sumw;
4145 Float_t correction =0;
4146 if (fNChange>0){
4147 // tracklet on boundary
4148 if (fMeanz<fZProb) correction = ycrosscor;
4149 if (fMeanz>fZProb) correction = -ycrosscor;
4150 }
4151 Double_t det = sumw*sumwx2-sumwx*sumwx;
4152 fYfitR[0] = (sumwx2*sumwy-sumwx*sumwxy)/det;
4153 fYfitR[1] = (sumw*sumwxy-sumwx*sumwy)/det;
4154 //
4155 fSigmaY2 =0;
4156 for (Int_t i=0;i<25;i++){
4157 if (!fUsable[i]) continue;
4158 Float_t delta = yres[i]-fYfitR[0]-fYfitR[1]*fX[i];
4159 fSigmaY2+=delta*delta;
4160 }
4161 fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
4162 //
4163 fZfitR[0] = (sumwx2*sumwz-sumwx*sumwxz)/det;
4164 fZfitR[1] = (sumw*sumwxz-sumwx*sumwz)/det;
4165 fZfit[0] = (sumwx2*sumwz-sumwx*sumwxz)/det;
4166 fZfit[1] = (sumw*sumwxz-sumwx*sumwz)/det;
4167 fYfitR[0] += fYref[0]+correction;
4168 fYfitR[1] += fYref[1];
4169 fYfit[0] = fYfitR[0];
4170 fYfit[1] = fYfitR[1];
4171 //
4172 //
4173 UpdateUsed();
4174}
4175
4176
4177
4178
4179
4180
4181void AliTRDseed::UpdateUsed(){
4182 //
4183 fNUsed =0;
4184 for (Int_t i=0;i<25;i++){
4185 if (!fClusters[i]) continue;
4186 if ((fClusters[i]->IsUsed())) fNUsed++;
4187 }
4188}
4189
4190
4191void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
4192{
4193 //
4194 // robust estimator in 1D case MI version
4195 //
4196 //for the univariate case
4197 //estimates of location and scatter are returned in mean and sigma parameters
4198 //the algorithm works on the same principle as in multivariate case -
4199 //it finds a subset of size hh with smallest sigma, and then returns mean and
4200 //sigma of this subset
4201
4202 if (hh==0)
4203 hh=(nvectors+2)/2;
4204 Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
4205 Int_t *index=new Int_t[nvectors];
4206 TMath::Sort(nvectors, data, index, kFALSE);
4207 //
4208 Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
4209 Double_t factor = faclts[nquant-1];
4210 //
4211 //
4212 Double_t sumx =0;
4213 Double_t sumx2 =0;
4214 Int_t bestindex = -1;
4215 Double_t bestmean = 0;
4216 Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
4217 for (Int_t i=0; i<hh; i++){
4218 sumx += data[index[i]];
4219 sumx2 += data[index[i]]*data[index[i]];
4220 }
4221 //
4222 Double_t norm = 1./Double_t(hh);
4223 Double_t norm2 = 1./Double_t(hh-1);
4224 for (Int_t i=hh; i<nvectors; i++){
4225 Double_t cmean = sumx*norm;
4226 Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
4227 if (csigma<bestsigma){
4228 bestmean = cmean;
4229 bestsigma = csigma;
4230 bestindex = i-hh;
4231 }
4232 //
4233 //
4234 sumx += data[index[i]]-data[index[i-hh]];
4235 sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
4236 }
4237
4238 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
4239 mean = bestmean;
4240 sigma = bstd;
4241 delete [] index;
4242}
4243
4244
4245Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror){
4246 //
4247 //
4248 //
4249 TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
4250 fitterT2.StoreData(kTRUE);
4251 Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
4252 //
4253 Int_t npointsT =0;
4254 fitterT2.ClearPoints();
4255 for (Int_t iLayer=0; iLayer<6;iLayer++){
4256 if (!cseed[iLayer].isOK()) continue;
4257 Double_t tilt = cseed[iLayer].fTilt;
4258
4259 for (Int_t itime=0;itime<25;itime++){
4260 if (!cseed[iLayer].fUsable[itime]) continue;
4261 Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2; // x relative to the midle chamber
4262 Double_t y = cseed[iLayer].fY[itime];
4263 Double_t z = cseed[iLayer].fZ[itime];
4264 // tilted rieman
4265 //
4266 Double_t uvt[6];
4267 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
4268 Double_t t = 1./(x2*x2+y*y);
4269 uvt[1] = t; // t
4270 uvt[0] = 2.*x2*uvt[1]; // u
4271 uvt[2] = 2.0*tilt*uvt[1];
4272 uvt[3] = 2.0*tilt*x*uvt[1];
4273 uvt[4] = 2.0*(y+tilt*z)*uvt[1];
4274 //
4275 Double_t error = 2*uvt[1];
4276 if (terror) error*=cseed[iLayer].fSigmaY;
4277 else {error *=0.2;} //default error
4278 fitterT2.AddPoint(uvt,uvt[4],error);
4279 npointsT++;
4280 }
4281 }
4282 fitterT2.Eval();
4283 Double_t rpolz0 = fitterT2.GetParameter(3);
4284 Double_t rpolz1 = fitterT2.GetParameter(4);
4285 //
4286 // linear fitter - not possible to make boundaries
4287 // non accept non possible z and dzdx combination
4288 //
4289 Bool_t acceptablez =kTRUE;
4290 for (Int_t iLayer=0; iLayer<6;iLayer++){
4291 if (cseed[iLayer].isOK()){
4292 Double_t zT2 = rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
4293 if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
4294 acceptablez = kFALSE;
4295 }
4296 }
4297 if (!acceptablez){
4298 Double_t zmf = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
4299 Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
4300 fitterT2.FixParameter(3,zmf);
4301 fitterT2.FixParameter(4,dzmf);
4302 fitterT2.Eval();
4303 fitterT2.ReleaseParameter(3);
4304 fitterT2.ReleaseParameter(4);
4305 rpolz0 = fitterT2.GetParameter(3);
4306 rpolz1 = fitterT2.GetParameter(4);
4307 }
4308 //
4309 Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
4310 Double_t params[3];
4311 params[0] = fitterT2.GetParameter(0);
4312 params[1] = fitterT2.GetParameter(1);
4313 params[2] = fitterT2.GetParameter(2);
4314 Double_t CR = 1+params[1]*params[1]-params[2]*params[0];
4315 for (Int_t iLayer = 0; iLayer<6;iLayer++){
4316 Double_t x = cseed[iLayer].fX0;
4317 Double_t y=0,dy=0, z=0, dz=0;
4318 // y
4319 Double_t res2 = (x*params[0]+params[1]);
4320 res2*=res2;
4321 res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
4322 if (res2>=0){
4323 res2 = TMath::Sqrt(res2);
4324 y = (1-res2)/params[0];
4325 }
4326 //dy
4327 Double_t x0 = -params[1]/params[0];
4328 if (-params[2]*params[0]+params[1]*params[1]+1>0){
4329 Double_t Rm1 = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
4330 if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)>0){
4331 Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
4332 if (params[0]<0) res*=-1.;
4333 dy = res;
4334 }
4335 }
4336 z = rpolz0+rpolz1*(x-xref2);
4337 dz = rpolz1;
4338 cseed[iLayer].fYref[0] = y;
4339 cseed[iLayer].fYref[1] = dy;
4340 cseed[iLayer].fZref[0] = z;
4341 cseed[iLayer].fZref[1] = dz;
4342 cseed[iLayer].fC = CR;
4343 //
4344 }
4345 return chi2TR;
4346}