]> 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
Possibility to use loadlibs.C on Solaris x86
[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// //
c6f438c0 20// The standard TRD tracker //
21// Based on Kalman filltering approach //
029cd327 22// //
23///////////////////////////////////////////////////////////////////////////////
24
a2cb5b3d 25#include <Riostream.h>
46d29e70 26#include <TFile.h>
46d29e70 27#include <TBranch.h>
5443e65e 28#include <TTree.h>
c630aafd 29#include <TObjArray.h>
46d29e70 30
46d29e70 31#include "AliTRDgeometry.h"
a5cadd36 32#include "AliTRDpadPlane.h"
0b2318ec 33#include "AliTRDgeometry.h"
46d29e70 34#include "AliTRDcluster.h"
35#include "AliTRDtrack.h"
75bd7f81 36#include "AliTRDseed.h"
6c94f330 37#include "AliESD.h"
38
3551db50 39#include "AliTRDcalibDB.h"
40#include "AliTRDCommonParam.h"
6c94f330 41
42#include "TTreeStream.h"
43#include "TGraph.h"
46d29e70 44#include "AliTRDtracker.h"
6c94f330 45#include "TLinearFitter.h"
46#include "AliRieman.h"
47#include "AliTrackPointArray.h"
48#include "AliAlignObj.h"
7b580082 49#include "AliTRDReconstructor.h"
46d29e70 50
51ClassImp(AliTRDtracker)
bbf92647 52
029cd327 53 const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5;
029cd327 54 const Float_t AliTRDtracker::fgkLabelFraction = 0.8;
029cd327 55 const Double_t AliTRDtracker::fgkMaxChi2 = 12.;
75bd7f81 56 const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // correspond to tan = 3
57 const Double_t AliTRDtracker::fgkMaxStep = 2.; // maximal step size in propagation
7ad19338 58
75bd7f81 59//_____________________________________________________________________________
6c94f330 60AliTRDtracker::AliTRDtracker():
61 AliTracker(),
62 fGeom(0),
63 fNclusters(0),
64 fClusters(0),
65 fNseeds(0),
66 fSeeds(0),
67 fNtracks(0),
68 fTracks(0),
69 fTimeBinsPerPlane(0),
70 fAddTRDseeds(kFALSE),
71 fNoTilt(kFALSE),
72 fDebugStreamer(0)
89f05372 73{
75bd7f81 74 //
b7a0917f 75 // Default constructor
75bd7f81 76 //
b7a0917f 77
6c94f330 78 for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
79 for(Int_t j=0;j<5;j++)
80 for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
75bd7f81 81
89f05372 82}
75bd7f81 83
84//_____________________________________________________________________________
6c94f330 85AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
86 :AliTracker(t)
ad670fba 87 ,fGeom(0)
88 ,fNclusters(0)
89 ,fClusters(0)
90 ,fNseeds(0)
91 ,fSeeds(0)
92 ,fNtracks(0)
93 ,fTracks(0)
94 ,fTimeBinsPerPlane(0)
95 ,fAddTRDseeds(kFALSE)
96 ,fNoTilt(kFALSE)
97 ,fDebugStreamer(0)
6c94f330 98{
ad670fba 99 //
100 // Copy constructor
101 //
ad670fba 102}
103
104//_____________________________________________________________________________
105AliTRDtracker::AliTRDtracker(const TFile *geomfile)
106 :AliTracker()
107 ,fGeom(0)
108 ,fNclusters(0)
109 ,fClusters(new TObjArray(2000))
110 ,fNseeds(0)
111 ,fSeeds(new TObjArray(2000))
112 ,fNtracks(0)
113 ,fTracks(new TObjArray(1000))
114 ,fTimeBinsPerPlane(0)
115 ,fAddTRDseeds(kFALSE)
116 ,fNoTilt(kFALSE)
117 ,fDebugStreamer(0)
46d29e70 118{
5443e65e 119 //
120 // Main constructor
121 //
b8dc2353 122
6c94f330 123 TDirectory *savedir=gDirectory;
124 TFile *in=(TFile*)geomfile;
5443e65e 125 if (!in->IsOpen()) {
6c94f330 126 printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
127 printf(" FULL TRD geometry and DEFAULT TRD parameter will be used\n");
5443e65e 128 }
129 else {
130 in->cd();
6c94f330 131 fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
5443e65e 132 }
46d29e70 133
6c94f330 134 if(fGeom) {
135 // printf("Found geometry version %d on file \n", fGeom->IsVersion());
136 }
137 else {
138 printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
0b2318ec 139 fGeom = new AliTRDgeometry();
c630aafd 140 }
c6f438c0 141 fGeom->ReadGeoMatrices();
c630aafd 142
5443e65e 143 savedir->cd();
46d29e70 144
6c94f330 145 for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
029cd327 146 Int_t trS = CookSectorIndex(geomS);
6c94f330 147 fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS);
148 for (Int_t icham=0;icham<AliTRDgeometry::kNcham; icham++){
149 fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
3c625a9b 150 }
5443e65e 151 }
3551db50 152 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
7ad19338 153 Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
6c94f330 154 if(tiltAngle < 0.1) {
b8dc2353 155 fNoTilt = kTRUE;
156 }
157
59393e34 158 fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
46d29e70 159
6c94f330 160 fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
0a29d0f1 161
9c9d2487 162 savedir->cd();
75bd7f81 163
5443e65e 164}
46d29e70 165
75bd7f81 166//_____________________________________________________________________________
5443e65e 167AliTRDtracker::~AliTRDtracker()
46d29e70 168{
029cd327 169 //
170 // Destructor of AliTRDtracker
171 //
172
89f05372 173 if (fClusters) {
174 fClusters->Delete();
175 delete fClusters;
176 }
177 if (fTracks) {
178 fTracks->Delete();
179 delete fTracks;
180 }
181 if (fSeeds) {
182 fSeeds->Delete();
183 delete fSeeds;
184 }
5443e65e 185 delete fGeom;
0a29d0f1 186
029cd327 187 for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
188 delete fTrSec[geomS];
5443e65e 189 }
7ad19338 190 if (fDebugStreamer) {
6c94f330 191 //fDebugStreamer->Close();
7ad19338 192 delete fDebugStreamer;
193 }
9c9d2487 194
75bd7f81 195}
59393e34 196
75bd7f81 197//_____________________________________________________________________________
198Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
199{
59393e34 200 //
75bd7f81 201 // Transform internal TRD ID to global detector ID
59393e34 202 //
75bd7f81 203
6c94f330 204 Int_t isector = fGeom->GetSector(lid);
205 Int_t ichamber= fGeom->GetChamber(lid);
206 Int_t iplan = fGeom->GetPlane(lid);
59393e34 207 //
208 AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
209 switch (iplan) {
210 case 0:
211 iLayer = AliAlignObj::kTRD1;
212 break;
213 case 1:
214 iLayer = AliAlignObj::kTRD2;
215 break;
216 case 2:
217 iLayer = AliAlignObj::kTRD3;
218 break;
219 case 3:
220 iLayer = AliAlignObj::kTRD4;
221 break;
222 case 4:
223 iLayer = AliAlignObj::kTRD5;
224 break;
225 case 5:
226 iLayer = AliAlignObj::kTRD6;
227 break;
228 };
6c94f330 229 Int_t modId = isector*fGeom->Ncham()+ichamber;
59393e34 230 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
75bd7f81 231
59393e34 232 return volid;
75bd7f81 233
59393e34 234}
235
75bd7f81 236//_____________________________________________________________________________
237Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
238{
59393e34 239 //
75bd7f81 240 // Transform global detector ID to local detector ID
59393e34 241 //
75bd7f81 242
6c94f330 243 Int_t modId=0;
244 AliAlignObj::ELayerID layerId = AliAlignObj::VolUIDToLayer(gid, modId);
245 Int_t isector = modId/fGeom->Ncham();
246 Int_t ichamber = modId%fGeom->Ncham();
247 Int_t iLayer = -1;
59393e34 248 switch (layerId) {
249 case AliAlignObj::kTRD1:
6c94f330 250 iLayer = 0;
59393e34 251 break;
252 case AliAlignObj::kTRD2:
6c94f330 253 iLayer = 1;
59393e34 254 break;
255 case AliAlignObj::kTRD3:
6c94f330 256 iLayer = 2;
59393e34 257 break;
258 case AliAlignObj::kTRD4:
6c94f330 259 iLayer = 3;
59393e34 260 break;
261 case AliAlignObj::kTRD5:
6c94f330 262 iLayer = 4;
59393e34 263 break;
264 case AliAlignObj::kTRD6:
6c94f330 265 iLayer = 5;
59393e34 266 break;
267 default:
6c94f330 268 iLayer =-1;
59393e34 269 }
6c94f330 270 if (iLayer<0) return -1;
59393e34 271 Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
75bd7f81 272
59393e34 273 return lid;
59393e34 274
75bd7f81 275}
59393e34 276
75bd7f81 277//_____________________________________________________________________________
6c94f330 278Bool_t AliTRDtracker::Transform(AliTRDcluster * cluster)
75bd7f81 279{
59393e34 280 //
75bd7f81 281 // Transform something ... whatever ...
c6f438c0 282 //
75bd7f81 283
6c94f330 284 const Double_t kX0shift = 2.52; // magic constants for geo manager transformation
285 const Double_t kX0shift5 = 3.05; //
59393e34 286 //
59393e34 287 //
6c94f330 288 // apply alignment and calibration to transform cluster
289 //
290 //
291 Int_t detector = cluster->GetDetector();
292 Int_t plane = fGeom->GetPlane(cluster->GetDetector());
293 Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
294 Int_t sector = fGeom->GetSector(cluster->GetDetector());
c6f438c0 295
6c94f330 296 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
297 Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance
59393e34 298 //
59393e34 299 // ExB correction
300 //
301 Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
6c94f330 302 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
303 //
304 AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
305 AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
c6f438c0 306 Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
6c94f330 307 Double_t localPos[3], localPosTracker[3];
c6f438c0 308 localPos[0] = -cluster->GetX();
309 localPos[1] = cluster->GetY() - driftX*exB;
6c94f330 310 localPos[2] = cluster->GetZ() -zshiftIdeal;
311 //
c6f438c0 312 cluster->SetY(cluster->GetY() - driftX*exB);
313 Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
314 cluster->SetX(xplane- cluster->GetX());
6c94f330 315 //
316 TGeoHMatrix * matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector());
c6f438c0 317 if (!matrix){
6c94f330 318 // no matrix found - if somebody used geometry with holes
c6f438c0 319 AliError("Invalid Geometry - Default Geometry used\n");
320 return kTRUE;
321 }
322 matrix->LocalToMaster(localPos, localPosTracker);
6c94f330 323 //
324 //
325 //
326 if (AliTRDReconstructor::StreamLevel()>1){
327 (*fDebugStreamer)<<"Transform"<<
328 "Cl.="<<cluster<<
329 "matrix.="<<matrix<<
330 "Detector="<<detector<<
331 "Sector="<<sector<<
332 "Plane="<<plane<<
333 "Chamber="<<chamber<<
334 "lx0="<<localPosTracker[0]<<
335 "ly0="<<localPosTracker[1]<<
336 "lz0="<<localPosTracker[2]<<
337 "\n";
c6f438c0 338 }
6c94f330 339 //
340 if (plane==5)
c6f438c0 341 cluster->SetX(localPosTracker[0]+kX0shift5);
6c94f330 342 else
c6f438c0 343 cluster->SetX(localPosTracker[0]+kX0shift);
6c94f330 344
c6f438c0 345 cluster->SetY(localPosTracker[1]);
346 cluster->SetZ(localPosTracker[2]);
75bd7f81 347
59393e34 348 return kTRUE;
75bd7f81 349
59393e34 350}
351
75bd7f81 352//_____________________________________________________________________________
353// Bool_t AliTRDtracker::Transform(AliTRDcluster * cluster)
354//{
c6f438c0 355// //
356// //
357// const Double_t kDriftCorrection = 1.01; // drift coeficient correction
358// const Double_t kTime0Cor = 0.32; // time0 correction
359// //
360// const Double_t kX0shift = 2.52;
361// const Double_t kX0shift5 = 3.05;
362
363// //
364// // apply alignment and calibration to transform cluster
365// //
366// //
367// Int_t detector = cluster->GetDetector();
368// Int_t plane = fGeom->GetPlane(cluster->GetDetector());
369// Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
370// Int_t sector = fGeom->GetSector(cluster->GetDetector());
371
372// Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
373// Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance
374// //
375// // ExB correction
376// //
377// Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
378// Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
379// //
380
381// AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
382// AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
383// Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
384// Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3];
385// localPos[2] = -cluster->GetX();
386// localPos[0] = cluster->GetY() - driftX*exB;
387// localPos[1] = cluster->GetZ() -zshiftIdeal;
388// TGeoHMatrix * matrix = fGeom->GetGeoMatrix(cluster->GetDetector());
389// matrix->LocalToMaster(localPos, globalPos);
390
391// Double_t sectorAngle = 20.*(sector%18)+10;
392// TGeoHMatrix rotSector;
393// rotSector.RotateZ(sectorAngle);
394// rotSector.LocalToMaster(globalPos, localPosTracker);
395// //
396// //
397// TGeoHMatrix matrix2(*matrix);
398// matrix2.MultiplyLeft(&rotSector);
399// matrix2.LocalToMaster(localPos,localPosTracker2);
400// //
401// //
402// //
403// cluster->SetY(cluster->GetY() - driftX*exB);
404// Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
405// cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
406// (*fDebugStreamer)<<"Transform"<<
407// "Cl.="<<cluster<<
408// "matrix.="<<matrix<<
409// "matrix2.="<<&matrix2<<
410// "Detector="<<detector<<
411// "Sector="<<sector<<
412// "Plane="<<plane<<
413// "Chamber="<<chamber<<
414// "lx0="<<localPosTracker[0]<<
415// "ly0="<<localPosTracker[1]<<
416// "lz0="<<localPosTracker[2]<<
417// "lx2="<<localPosTracker2[0]<<
418// "ly2="<<localPosTracker2[1]<<
419// "lz2="<<localPosTracker2[2]<<
420// "\n";
421// //
422// if (plane==5)
423// cluster->SetX(localPosTracker[0]+kX0shift5);
424// else
425// cluster->SetX(localPosTracker[0]+kX0shift);
426
427// cluster->SetY(localPosTracker[1]);
428// cluster->SetZ(localPosTracker[2]);
429// return kTRUE;
430// }
431
75bd7f81 432//_____________________________________________________________________________
433Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
434{
9c9d2487 435 //
436 // Rotates the track when necessary
437 //
438
439 Double_t alpha = AliTRDgeometry::GetAlpha();
6c94f330 440 Double_t y = track->GetY();
441 Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
9c9d2487 442
c630aafd 443 //Int_t ns = AliTRDgeometry::kNsect;
9c9d2487 444 //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
445
6c94f330 446 if (y > ymax) {
9c9d2487 447 //s = (s+1) % ns;
6c94f330 448 if (!track->Rotate(alpha)) return kFALSE;
449 } else if (y <-ymax) {
9c9d2487 450 //s = (s-1+ns) % ns;
6c94f330 451 if (!track->Rotate(-alpha)) return kFALSE;
9c9d2487 452 }
453
454 return kTRUE;
9c9d2487 455
75bd7f81 456}
46e2d86c 457
75bd7f81 458//_____________________________________________________________________________
459AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
460 , Int_t timebin, UInt_t &index)
461{
46e2d86c 462 //
75bd7f81 463 // Try to find cluster in the backup list
46e2d86c 464 //
75bd7f81 465
6c94f330 466 AliTRDcluster * cl =0;
467 Int_t *indexes = track->GetBackupIndexes();
468 for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
469 if (indexes[i]==0) break;
470 AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
471 if (!cli) break;
472 if (cli->GetLocalTimeBin()!=timebin) continue;
46e2d86c 473 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
6c94f330 474 if (iplane==plane) {
475 cl = cli;
7ad19338 476 index = indexes[i];
46e2d86c 477 break;
478 }
479 }
75bd7f81 480
46e2d86c 481 return cl;
46e2d86c 482
75bd7f81 483}
3c625a9b 484
75bd7f81 485//_____________________________________________________________________________
6c94f330 486Int_t AliTRDtracker::GetLastPlane(AliTRDtrack * track)
75bd7f81 487{
3c625a9b 488 //
75bd7f81 489 // Return last updated plane
490 //
491
6c94f330 492 Int_t lastplane=0;
493 Int_t *indexes = track->GetBackupIndexes();
494 for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
495 AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
496 if (!cli) break;
3c625a9b 497 Int_t iplane = fGeom->GetPlane(cli->GetDetector());
6c94f330 498 if (iplane>lastplane) {
3c625a9b 499 lastplane = iplane;
500 }
501 }
75bd7f81 502
3c625a9b 503 return lastplane;
75bd7f81 504
3c625a9b 505}
75bd7f81 506
507//_____________________________________________________________________________
c630aafd 508Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
509{
510 //
511 // Finds tracks within the TRD. The ESD event is expected to contain seeds
512 // at the outer part of the TRD. The seeds
513 // are found within the TRD if fAddTRDseeds is TRUE.
514 // The tracks are propagated to the innermost time bin
515 // of the TRD and the ESD event is updated
516 //
517
6c94f330 518 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
029cd327 519 Float_t foundMin = fgkMinClustersInTrack * timeBins;
6c94f330 520 Int_t nseed = 0;
521 Int_t found = 0;
522 // Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
c630aafd 523
524 Int_t n = event->GetNumberOfTracks();
6c94f330 525 for (Int_t i=0; i<n; i++) {
526 AliESDtrack* seed=event->GetTrack(i);
527 ULong_t status=seed->GetStatus();
528 if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
529 if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
c630aafd 530 nseed++;
7ad19338 531
6c94f330 532 AliTRDtrack* seed2 = new AliTRDtrack(*seed);
46e2d86c 533 //seed2->ResetCovariance();
c630aafd 534 AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
6c94f330 535 //AliTRDtrack *pt = new AliTRDtrack(*seed2); // !!!
c630aafd 536 AliTRDtrack &t=*pt;
7b580082 537 FollowProlongation(t);
c630aafd 538 if (t.GetNumberOfClusters() >= foundMin) {
539 UseClusters(&t);
6c94f330 540 CookLabel(pt, 1-fgkLabelFraction);
541 // t.CookdEdx();
c630aafd 542 }
543 found++;
6c94f330 544// cout<<found<<'\r';
c630aafd 545
59393e34 546 Double_t xTPC = 250;
547 if (PropagateToX(t,xTPC,fgkMaxStep)) {
c630aafd 548 seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
549 }
550 delete seed2;
551 delete pt;
552 }
553
6c94f330 554 cout<<"Number of loaded seeds: "<<nseed<<endl;
555 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
556
557 // after tracks from loaded seeds are found and the corresponding
558 // clusters are used, look for additional seeds from TRD
559
560
561 cout<<"Total number of found tracks: "<<found<<endl;
c630aafd 562
563 return 0;
75bd7f81 564
c630aafd 565}
5443e65e 566
c630aafd 567//_____________________________________________________________________________
6c94f330 568Int_t AliTRDtracker::PropagateBack(AliESD* event)
75bd7f81 569{
c630aafd 570 //
571 // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
572 // backpropagated by the TPC tracker. Each seed is first propagated
573 // to the TRD, and then its prolongation is searched in the TRD.
574 // If sufficiently long continuation of the track is found in the TRD
575 // the track is updated, otherwise it's stored as originaly defined
576 // by the TPC tracker.
577 //
578
6c94f330 579 Int_t found=0;
c5a8e3df 580 Float_t foundMin = 20;
6c94f330 581 Int_t n = event->GetNumberOfTracks();
582 //
583 //Sort tracks
584 Float_t *quality =new Float_t[n];
585 Int_t *index =new Int_t[n];
586 for (Int_t i=0; i<n; i++) {
587 AliESDtrack* seed=event->GetTrack(i);
4f1c04d3 588 Double_t covariance[15];
589 seed->GetExternalCovariance(covariance);
590 quality[i] = covariance[0]+covariance[2];
591 }
592 TMath::Sort(n,quality,index,kFALSE);
6c94f330 593 //
594 for (Int_t i=0; i<n; i++) {
595 // AliESDtrack* seed=event->GetTrack(i);
4f1c04d3 596 AliESDtrack* seed=event->GetTrack(index[i]);
597
6c94f330 598 ULong_t status=seed->GetStatus();
599 if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
600 if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
c630aafd 601
602 Int_t lbl = seed->GetLabel();
603 AliTRDtrack *track = new AliTRDtrack(*seed);
604 track->SetSeedLabel(lbl);
6c94f330 605 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); //make backup
c630aafd 606 fNseeds++;
6c94f330 607 Float_t p4 = track->GetC();
608 //
f6625211 609 Int_t expectedClr = FollowBackProlongation(*track);
6c94f330 610 if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)<0.2 || TMath::Abs(track->GetPt())>0.8 ) {
4f1c04d3 611 //
6c94f330 612 //make backup for back propagation
4f1c04d3 613 //
614 Int_t foundClr = track->GetNumberOfClusters();
615 if (foundClr >= foundMin) {
616 track->CookdEdx();
8979685e 617 CookdEdxTimBin(*track);
6c94f330 618 CookLabel(track, 1-fgkLabelFraction);
619 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
620 if(track->GetChi2()/track->GetNumberOfClusters()<4) { // sign only gold tracks
621 if (seed->GetKinkIndex(0)==0&&TMath::Abs(track->GetPt())<1.5 ) UseClusters(track);
4f1c04d3 622 }
623 Bool_t isGold = kFALSE;
624
6c94f330 625 if (track->GetChi2()/track->GetNumberOfClusters()<5) { //full gold track
626 // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
627 if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
4f1c04d3 628 isGold = kTRUE;
629 }
6c94f330 630 if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track
631 // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
632 if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
f4e9508c 633 isGold = kTRUE;
634 }
6c94f330 635 if (!isGold && track->GetBackupTrack()){
636 if (track->GetBackupTrack()->GetNumberOfClusters()>foundMin&&
637 (track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1))<7){
638 seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
4f1c04d3 639 isGold = kTRUE;
640 }
641 }
6c94f330 642 if (track->StatusForTOF()>0 &&track->fNCross==0 && Float_t(track->fN)/Float_t(track->fNExpected)>0.4){
643 //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
7ad19338 644 }
16d9fbba 645 }
c630aafd 646 }
8979685e 647 // Debug part of tracking
6c94f330 648 TTreeSRedirector& cstream = *fDebugStreamer;
8979685e 649 Int_t eventNr = event->GetEventNumber();
6c94f330 650 if (AliTRDReconstructor::StreamLevel()>0){
651 if (track->GetBackupTrack()){
652 cstream<<"Tracks"<<
653 "EventNr="<<eventNr<<
654 "ESD.="<<seed<<
655 "trd.="<<track<<
656 "trdback.="<<track->GetBackupTrack()<<
657 "\n";
658 }else{
659 cstream<<"Tracks"<<
660 "EventNr="<<eventNr<<
661 "ESD.="<<seed<<
662 "trd.="<<track<<
663 "trdback.="<<track<<
664 "\n";
d337ef8d 665 }
8979685e 666 }
6c94f330 667 //
668 //Propagation to the TOF (I.Belikov)
669 if (track->GetStop()==kFALSE){
4f1c04d3 670
6c94f330 671 Double_t xtof=371.;
672 Double_t c2=track->GetSnp() + track->GetC()*(xtof - track->GetX());
673 if (TMath::Abs(c2)>=0.99) {
c5a8e3df 674 delete track;
675 continue;
676 }
59393e34 677 Double_t xTOF0 = 370. ;
678 PropagateToX(*track,xTOF0,fgkMaxStep);
6c94f330 679 //
680 //energy losses taken to the account - check one more time
681 c2=track->GetSnp() + track->GetC()*(xtof - track->GetX());
682 if (TMath::Abs(c2)>=0.99) {
4f1c04d3 683 delete track;
684 continue;
685 }
686
6c94f330 687 //
688 Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
689 Double_t y; track->GetYAt(xtof,GetBz(),y);
690 if (y > ymax) {
7ac6fa52 691 if (!track->Rotate(AliTRDgeometry::GetAlpha())) {
692 delete track;
7bed16a7 693 continue;
7ac6fa52 694 }
6c94f330 695 } else if (y <-ymax) {
7ac6fa52 696 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
697 delete track;
7bed16a7 698 continue;
7ac6fa52 699 }
3c625a9b 700 }
701
702 if (track->PropagateTo(xtof)) {
6c94f330 703 seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
704 for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
705 for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
6d45eaef 706 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
707 }
708 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 709 }
6c94f330 710 // seed->SetTRDtrack(new AliTRDtrack(*track));
711 if (track->GetNumberOfClusters()>foundMin) found++;
3c625a9b 712 }
6c94f330 713 }else{
714 if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
715 seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
16d9fbba 716 //seed->SetStatus(AliESDtrack::kTRDStop);
6c94f330 717 for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
718 for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
6d45eaef 719 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
720 }
721 seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
eab5961e 722 }
7ad19338 723 //seed->SetTRDtrack(new AliTRDtrack(*track));
3c625a9b 724 found++;
725 }
1e9bb598 726 }
7ad19338 727 seed->SetTRDQuality(track->StatusForTOF());
8979685e 728 seed->SetTRDBudget(track->fBudget[0]);
729
d9b8978b 730 delete track;
6c94f330 731 //
732 //End of propagation to the TOF
733 //if (foundClr>foundMin)
734 // seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
735
c630aafd 736
737 }
ad670fba 738
6c94f330 739 cerr<<"Number of seeds: "<<fNseeds<<endl;
740 cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
741
742 if (AliTRDReconstructor::SeedingOn()) MakeSeedsMI(3,5,event); //new seeding
7ad19338 743
6c94f330 744 fSeeds->Clear(); fNseeds=0;
4f1c04d3 745 delete [] index;
746 delete [] quality;
747
1e9bb598 748 return 0;
749
750}
751
752//_____________________________________________________________________________
6c94f330 753Int_t AliTRDtracker::RefitInward(AliESD* event)
1e9bb598 754{
755 //
756 // Refits tracks within the TRD. The ESD event is expected to contain seeds
757 // at the outer part of the TRD.
758 // The tracks are propagated to the innermost time bin
759 // of the TRD and the ESD event is updated
760 // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
761 //
762
6c94f330 763 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
1e9bb598 764 Float_t foundMin = fgkMinClustersInTrack * timeBins;
6c94f330 765 Int_t nseed = 0;
766 Int_t found = 0;
767 // Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
4f1c04d3 768 AliTRDtrack seed2;
6c94f330 769
1e9bb598 770 Int_t n = event->GetNumberOfTracks();
6c94f330 771 for (Int_t i=0; i<n; i++) {
772 AliESDtrack* seed=event->GetTrack(i);
4f1c04d3 773 new(&seed2) AliTRDtrack(*seed);
6c94f330 774 if (seed2.GetX()<270){
775 seed->UpdateTrackParams(&seed2, AliESDtrack::kTRDbackup); // backup TPC track - only update
f4e9508c 776 continue;
777 }
778
6c94f330 779 ULong_t status=seed->GetStatus();
780 if ( (status & AliESDtrack::kTRDout ) == 0 ) {
0dd7d129 781 continue;
782 }
6c94f330 783 if ( (status & AliESDtrack::kTRDin) != 0 ) {
0dd7d129 784 continue;
785 }
6c94f330 786 nseed++;
787
788 seed2.ResetCovariance(50.);
789
790// if (1/seed2.Get1Pt()>1.5&& seed2.GetX()>260.) {
791// Double_t oldx = seed2.GetX();
792// seed2.PropagateTo(500.);
793// seed2.ResetCovariance(1.);
794// seed2.PropagateTo(oldx);
795// }
796// else{
797// seed2.ResetCovariance(5.);
798// }
4f1c04d3 799
800 AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
6c94f330 801 //AliTRDtrack *pt = new AliTRDtrack(seed2); // !!!
802 Int_t * indexes2 = seed2.GetIndexes();
803 for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
804 for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
6d45eaef 805 pt->SetPIDsignals(seed2.GetPIDsignals(i,j),i,j);
806 }
7ad19338 807 pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
808 }
eab5961e 809
6c94f330 810 Int_t * indexes3 = pt->GetBackupIndexes();
811 for (Int_t i=0;i<200;i++) {
812 if (indexes2[i]==0) break;
46e2d86c 813 indexes3[i] = indexes2[i];
814 }
815 //AliTRDtrack *pt = seed2;
6c94f330 816 AliTRDtrack &t=*pt;
7b580082 817 FollowProlongation(t);
1e9bb598 818 if (t.GetNumberOfClusters() >= foundMin) {
6c94f330 819 // UseClusters(&t);
820 //CookLabel(pt, 1-fgkLabelFraction);
7ad19338 821 t.CookdEdx();
822 CookdEdxTimBin(t);
1e9bb598 823 }
824 found++;
6c94f330 825// cout<<found<<'\r';
59393e34 826 Double_t xTPC = 250;
6c94f330 827 if(PropagateToX(t,xTPC,fgkMaxStep)) {
828 seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
829 for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
830 for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
6d45eaef 831 seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
832 }
7ad19338 833 seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
834 }
6c94f330 835 }else{
836 //if not prolongation to TPC - propagate without update
837 AliTRDtrack* seed2 = new AliTRDtrack(*seed);
7bed16a7 838 seed2->ResetCovariance(5.);
839 AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
6c94f330 840 //AliTRDtrack *pt2 = new AliTRDtrack(*seed2); // !!!
7bed16a7 841 delete seed2;
59393e34 842 if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
6c94f330 843 //pt2->CookdEdx(0.,1.);
844 pt2->CookdEdx( ); // Modification by PS
eab5961e 845 CookdEdxTimBin(*pt2);
6c94f330 846 seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit);
847 for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
848 for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
6d45eaef 849 seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
850 }
7ad19338 851 seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
852 }
eab5961e 853 }
7bed16a7 854 delete pt2;
1e9bb598 855 }
1e9bb598 856 delete pt;
eab5961e 857 }
1e9bb598 858
6c94f330 859 cout<<"Number of loaded seeds: "<<nseed<<endl;
860 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
1e9bb598 861
c630aafd 862 return 0;
863
864}
865
75bd7f81 866//_____________________________________________________________________________
6c94f330 867Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t)
8979685e 868{
75bd7f81 869 //
8979685e 870 // Starting from current position on track=t this function tries
871 // to extrapolate the track up to timeBin=0 and to confirm prolongation
872 // if a close cluster is found. Returns the number of clusters
873 // expected to be found in sensitive layers
874 // GeoManager used to estimate mean density
75bd7f81 875 //
876
6c94f330 877 Int_t sector;
878 Int_t lastplane = GetLastPlane(&t);
3551db50 879 Double_t radLength = 0.0;
6c94f330 880 Double_t rho = 0.0;
881 Int_t expectedNumberOfClusters = 0;
882 //
883 //
884 //
885 for (Int_t iplane = lastplane; iplane>=0; iplane--){
886 //
887 Int_t row0 = GetGlobalTimeBin(0, iplane,GetTimeBinsPerPlane()-1);
888 Int_t rowlast = GetGlobalTimeBin(0, iplane,0);
8979685e 889 //
6c94f330 890 // propagate track close to the plane if neccessary
7b580082 891 //
892 Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
6c94f330 893 if (currentx < -fgkMaxStep +t.GetX()){
894 //propagate closer to chamber - safety space fgkMaxStep
895 if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
8979685e 896 }
8979685e 897 if (!AdjustSector(&t)) break;
59393e34 898 //
6c94f330 899 // get material budget
8979685e 900 //
6c94f330 901 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
902 t.GetXYZ(xyz0); //starting global position
903 // end global position
7b580082 904 x = fTrSec[0]->GetLayer(row0)->GetX();
905 if (!t.GetProlongation(x,y,z)) break;
6c94f330 906 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
907 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
908 xyz1[2] = z;
7b580082 909 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
6c94f330 910 rho = param[0];
911 radLength = param[1]; // get mean propagation parameters
8979685e 912 //
6c94f330 913 // propagate nad update
8979685e 914 //
8979685e 915 sector = t.GetSector();
6c94f330 916 // for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
917 for (Int_t itime=0 ;itime<GetTimeBinsPerPlane();itime++) {
918 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 919 expectedNumberOfClusters++;
920 t.fNExpected++;
6c94f330 921 if (t.GetX()>345) t.fNExpectedLast++;
922 AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
923 AliTRDcluster *cl=0;
924 UInt_t index=0;
925 Double_t maxChi2=fgkMaxChi2;
8979685e 926 x = timeBin.GetX();
8979685e 927 if (timeBin) {
6c94f330 928 AliTRDcluster * cl0 = timeBin[0];
929 if (!cl0) continue; // no clusters in given time bin
930 Int_t plane = fGeom->GetPlane(cl0->GetDetector());
931 if (plane>lastplane) continue;
8979685e 932 Int_t timebin = cl0->GetLocalTimeBin();
6c94f330 933 AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index);
934 //
8979685e 935 if (cl2) {
6c94f330 936 cl =cl2;
937 //Double_t h01 = GetTiltFactor(cl);
938 //maxChi2=t.GetPredictedChi2(cl,h01);
939 //
940 //if (t.GetLabel()==277)
941 // printf("%e %e %x\n",t.GetX(),maxChi2,(UInt_t)cl);
942 //
7b580082 943 }
8979685e 944 if (cl) {
6c94f330 945 // if (cl->GetNPads()<5)
59393e34 946 Double_t dxsample = timeBin.GetdX();
947 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
6c94f330 948 Double_t h01 = GetTiltFactor(cl);
949 Int_t det = cl->GetDetector();
950 Int_t plane = fGeom->GetPlane(det);
951 if (t.GetX()>345){
8979685e 952 t.fNLast++;
953 t.fChi2Last+=maxChi2;
954 }
59393e34 955 Double_t xcluster = cl->GetX();
956 t.PropagateTo(xcluster,radLength,rho);
6c94f330 957 if (!AdjustSector(&t)) break;
958 maxChi2=t.GetPredictedChi2(cl,h01);
959
960 if (t.GetLabel()==277)
961 printf("%e %e %e %e\n",t.GetAlpha(),t.GetX(),t.GetY(),t.GetZ());
962
963 if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
7b580082 964 }
6c94f330 965 }
8979685e 966 }
6c94f330 967 }
8979685e 968 }
8979685e 969
75bd7f81 970 return expectedNumberOfClusters;
69b55c55 971
75bd7f81 972}
7b580082 973
75bd7f81 974//_____________________________________________________________________________
f6625211 975Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
8979685e 976{
75bd7f81 977 //
8979685e 978 // Starting from current radial position of track <t> this function
979 // extrapolates the track up to outer timebin and in the sensitive
980 // layers confirms prolongation if a close cluster is found.
981 // Returns the number of clusters expected to be found in sensitive layers
982 // Use GEO manager for material Description
75bd7f81 983 //
59393e34 984
6c94f330 985 Int_t sector;
986 Int_t clusters[1000];
987 for (Int_t i=0;i<1000;i++) clusters[i]=-1;
988 Double_t radLength = 0.0;
989 Double_t rho = 0.0;
990 Int_t expectedNumberOfClusters = 0;
991 Float_t ratio0=0;
8979685e 992 AliTRDtracklet tracklet;
6c94f330 993 //
994 //
995 for (Int_t iplane = 0; iplane<AliESDtrack::kNPlane; iplane++){
996 Int_t row0 = GetGlobalTimeBin(0, iplane,GetTimeBinsPerPlane()-1);
997 Int_t rowlast = GetGlobalTimeBin(0, iplane,0);
998 //
999 Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1000 if (currentx<t.GetX()) continue;
1001 //
1002 // propagate closer to chamber if neccessary
1003 //
1004 if (currentx > fgkMaxStep +t.GetX()){
1005 if (!PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
8979685e 1006 }
8979685e 1007 if (!AdjustSector(&t)) break;
59393e34 1008 if (TMath::Abs(t.GetSnp())>fgkMaxSnp) break;
8979685e 1009 //
6c94f330 1010 // get material budget inside of chamber
59393e34 1011 //
6c94f330 1012 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1013 t.GetXYZ(xyz0); //starting global position
1014 // end global position
7b580082 1015 x = fTrSec[0]->GetLayer(rowlast)->GetX();
1016 if (!t.GetProlongation(x,y,z)) break;
6c94f330 1017 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
1018 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
1019 xyz1[2] = z;
7b580082 1020 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
6c94f330 1021 rho = param[0];
1022 radLength = param[1]; // get mean propagation parameters
8979685e 1023 //
7b580082 1024 // Find clusters
8979685e 1025 //
6c94f330 1026 sector = t.GetSector();
1027 Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
1028 if (tracklet.GetN()<GetTimeBinsPerPlane()/3) continue;
8979685e 1029 //
7b580082 1030 // Propagate and update track
8979685e 1031 //
6c94f330 1032 for (Int_t itime= GetTimeBinsPerPlane()-1;itime>=0;itime--) {
1033 Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
8979685e 1034 expectedNumberOfClusters++;
1035 t.fNExpected++;
6c94f330 1036 if (t.GetX()>345) t.fNExpectedLast++;
1037 AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
1038 AliTRDcluster *cl=0;
1039 UInt_t index=0;
1040 Double_t maxChi2=fgkMaxChi2;
8979685e 1041 x = timeBin.GetX();
6c94f330 1042 //
8979685e 1043 if (timeBin) {
6c94f330 1044 if (clusters[ilayer]>0) {
8979685e 1045 index = clusters[ilayer];
1046 cl = (AliTRDcluster*)GetCluster(index);
6c94f330 1047 //Double_t h01 = GetTiltFactor(cl);
1048 //maxChi2=t.GetPredictedChi2(cl,h01);
1049 //
1050 //if (t.GetLabel()==277)
1051 // printf("%e %e %x\n",t.GetX(),maxChi2,(UInt_t)cl);
1052 //
8979685e 1053 }
1054
1055 if (cl) {
6c94f330 1056 // if (cl->GetNPads()<5)
59393e34 1057 Double_t dxsample = timeBin.GetdX();
1058 t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
6c94f330 1059 Double_t h01 = GetTiltFactor(cl);
1060 Int_t det = cl->GetDetector();
1061 Int_t plane = fGeom->GetPlane(det);
1062 if (t.GetX()>345){
8979685e 1063 t.fNLast++;
6c94f330 1064 t.fChi2Last+=maxChi2;
8979685e 1065 }
59393e34 1066 Double_t xcluster = cl->GetX();
1067 t.PropagateTo(xcluster,radLength,rho);
6c94f330 1068 maxChi2=t.GetPredictedChi2(cl,h01);
1069
1070 if (t.GetLabel()==277)
1071 printf("%e %e %e %e\n",t.GetAlpha(),t.GetX(),t.GetY(),t.GetZ());
1072
8979685e 1073 if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1074 if(!t.Update(cl,maxChi2,index,h01)) {
8979685e 1075 }
1076 }
6c94f330 1077 //
1078 // reset material budget if 2 consecutive gold
1079 if (plane>0)
1080 if (t.fTracklets[plane].GetN()+t.fTracklets[plane-1].GetN()>20){
8979685e 1081 t.fBudget[2] = 0;
1082 }
6c94f330 1083 }
8979685e 1084 }
59393e34 1085 }
6c94f330 1086 ratio0 = ncl/Float_t(fTimeBinsPerPlane);
1087 Float_t ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.);
1088 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){
1089 t.MakeBackupTrack(); // make backup of the track until is gold
59393e34 1090 }
7b580082 1091
8979685e 1092 }
5443e65e 1093
75bd7f81 1094 return expectedNumberOfClusters;
1e9bb598 1095
75bd7f81 1096}
1e9bb598 1097
75bd7f81 1098//_____________________________________________________________________________
6c94f330 1099Int_t AliTRDtracker::PropagateToX(AliTRDtrack& t, Double_t xToGo, Double_t maxStep)
5443e65e 1100{
75bd7f81 1101 //
5443e65e 1102 // Starting from current radial position of track <t> this function
1103 // extrapolates the track up to radial position <xToGo>.
1104 // Returns 1 if track reaches the plane, and 0 otherwise
75bd7f81 1105 //
1106
59393e34 1107 const Double_t kEpsilon = 0.00001;
6c94f330 1108 // Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1109 Double_t xpos = t.GetX();
1110 Double_t dir = (xpos<xToGo) ? 1.:-1.;
1111 //
1112 while ( (xToGo-xpos)*dir > kEpsilon){
1113 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1114 //
1115 Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1116 t.GetXYZ(xyz0); //starting global position
1117 x = xpos+step;
1118 //
1119 if (!t.GetProlongation(x,y,z)) return 0; // no prolongation
1120 //
1121 xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
1122 xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
1123 xyz1[2] = z;
1124 //
59393e34 1125 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1126 if (!t.PropagateTo(x,param[1],param[0])) return 0;
1127 AdjustSector(&t);
1128 xpos = t.GetX();
5443e65e 1129 }
75bd7f81 1130
5443e65e 1131 return 1;
5443e65e 1132
59393e34 1133}
5443e65e 1134
c630aafd 1135//_____________________________________________________________________________
1136Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1137{
75bd7f81 1138 //
c630aafd 1139 // Fills clusters into TRD tracking_sectors
1140 // Note that the numbering scheme for the TRD tracking_sectors
1141 // differs from that of TRD sectors
75bd7f81 1142 //
1143
6c94f330 1144 cout<<"\n Read Sectors clusters"<<endl;
c630aafd 1145 if (ReadClusters(fClusters,cTree)) {
6c94f330 1146 Error("LoadClusters","Problem with reading the clusters !");
c630aafd 1147 return 1;
1148 }
6c94f330 1149 Int_t ncl=fClusters->GetEntriesFast();
1150 fNclusters=ncl;
1151 cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
c630aafd 1152
1153 UInt_t index;
6c94f330 1154 for (Int_t ichamber=0;ichamber<5;ichamber++)
1155 for (Int_t isector=0;isector<18;isector++){
1156 fHoles[ichamber][isector]=kTRUE;
3c625a9b 1157 }
3c625a9b 1158
ad670fba 1159
6c94f330 1160 while (ncl--) {
1161// printf("\r %d left ",ncl);
1162 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1163 Int_t detector=c->GetDetector();
1164 Int_t localTimeBin=c->GetLocalTimeBin();
1165 Int_t sector=fGeom->GetSector(detector);
1166 Int_t plane=fGeom->GetPlane(detector);
1167
029cd327 1168 Int_t trackingSector = CookSectorIndex(sector);
6c94f330 1169 if (c->GetLabel(0)>0){
3c625a9b 1170 Int_t chamber = fGeom->GetChamber(detector);
6c94f330 1171 fHoles[chamber][trackingSector]=kFALSE;
3c625a9b 1172 }
c630aafd 1173
6c94f330 1174 Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
c630aafd 1175 if(gtb < 0) continue;
029cd327 1176 Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
c630aafd 1177
6c94f330 1178 index=ncl;
1179 //
1180 // apply pos correction
59393e34 1181 Transform(c);
029cd327 1182 fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
c630aafd 1183 }
75bd7f81 1184
c630aafd 1185 return 0;
75bd7f81 1186
c630aafd 1187}
1188
5443e65e 1189//_____________________________________________________________________________
b7a0917f 1190void AliTRDtracker::UnloadClusters()
5443e65e 1191{
1192 //
1193 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1194 //
1195
6c94f330 1196 Int_t i, nentr;
5443e65e 1197
1198 nentr = fClusters->GetEntriesFast();
6c94f330 1199 for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
b7a0917f 1200 fNclusters = 0;
5443e65e 1201
1202 nentr = fSeeds->GetEntriesFast();
6c94f330 1203 for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
5443e65e 1204
1205 nentr = fTracks->GetEntriesFast();
6c94f330 1206 for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
5443e65e 1207
1208 Int_t nsec = AliTRDgeometry::kNsect;
6c94f330 1209
5443e65e 1210 for (i = 0; i < nsec; i++) {
1211 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1212 fTrSec[i]->GetLayer(pl)->Clear();
1213 }
1214 }
1215
1216}
1217
75bd7f81 1218//_____________________________________________________________________________
69b55c55 1219void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD * esd)
7ad19338 1220{
1221 //
1222 // Creates seeds using clusters between position inner plane and outer plane
1223 //
75bd7f81 1224
6c94f330 1225 const Double_t kMaxTheta = 1;
1226 const Double_t kMaxPhi = 2.0;
1227 //
1228 const Double_t kRoad0y = 6; // road for middle cluster
1229 const Double_t kRoad0z = 8.5; // road for middle cluster
1230 //
1231 const Double_t kRoad1y = 2; // road in y for seeded cluster
1232 const Double_t kRoad1z = 20; // road in z for seeded cluster
1233 //
1234 const Double_t kRoad2y = 3; // road in y for extrapolated cluster
1235 const Double_t kRoad2z = 20; // road in z for extrapolated cluster
c6f438c0 1236 const Int_t kMaxSeed = 3000;
6c94f330 1237 Int_t maxSec=AliTRDgeometry::kNsect;
7ad19338 1238
6c94f330 1239 //
1240 // linear fitters in planes
1241 TLinearFitter fitterTC(2,"hyp2"); // fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1242 TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
69b55c55 1243 fitterTC.StoreData(kTRUE);
1244 fitterT2.StoreData(kTRUE);
6c94f330 1245 AliRieman rieman(1000); // rieman fitter
1246 AliRieman rieman2(1000); // rieman fitter
7ad19338 1247 //
6c94f330 1248 // find the maximal and minimal layer for the planes
7ad19338 1249 //
1250 Int_t layers[6][2];
6c94f330 1251 AliTRDpropagationLayer* reflayers[6];
1252 for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;}
1253 for (Int_t ns=0;ns<maxSec;ns++){
1254 for (Int_t ilayer=0;ilayer<fTrSec[ns]->GetNumberOfLayers();ilayer++){
1255 AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer));
1256 if (layer==0) continue;
7ad19338 1257 Int_t det = layer[0]->GetDetector();
1258 Int_t plane = fGeom->GetPlane(det);
6c94f330 1259 if (ilayer<layers[plane][0]) layers[plane][0] = ilayer;
1260 if (ilayer>layers[plane][1]) layers[plane][1] = ilayer;
7ad19338 1261 }
1262 }
6c94f330 1263 //
3551db50 1264 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
6c94f330 1265 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
1266 Double_t hL[6]; // tilting angle
1267 Double_t xcl[6]; // x - position of reference cluster
1268 Double_t ycl[6]; // y - position of reference cluster
1269 Double_t zcl[6]; // z - position of reference cluster
1270 AliTRDcluster *cl[6]={0,0,0,0,0,0}; // seeding clusters
1271 Float_t padlength[6]={10,10,10,10,10,10}; //current pad-length
1272 Double_t chi2R =0, chi2Z=0;
1273 Double_t chi2RF =0, chi2ZF=0;
1274 //
1275 Int_t nclusters; // total number of clusters
1276 for (Int_t i=0;i<6;i++) {hL[i]=h01; if (i%2==1) hL[i]*=-1.;}
1277 //
1278 //
1279 // registered seed
c6f438c0 1280 AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1281 AliTRDseed *seed[kMaxSeed];
6c94f330 1282 for (Int_t iseed=0;iseed<kMaxSeed;iseed++) seed[iseed]= &pseed[iseed*6];
69b55c55 1283 AliTRDseed *cseed = seed[0];
6c94f330 1284 //
1285 Double_t seedquality[kMaxSeed];
1286 Double_t seedquality2[kMaxSeed];
1287 Double_t seedparams[kMaxSeed][7];
1288 Int_t seedlayer[kMaxSeed];
1289 Int_t registered =0;
1290 Int_t sort[kMaxSeed];
1291 //
1292 // seeding part
1293 //
1294 for (Int_t ns = 0; ns<maxSec; ns++){ //loop over sectors
1295 //for (Int_t ns = 0; ns<5; ns++){ //loop over sectors
1296 registered = 0; // reset registerd seed counter
1297 cseed = seed[registered];
1298 Float_t iter=0;
1299 for (Int_t sLayer=2; sLayer>=0;sLayer--){
1300 //for (Int_t dseed=5;dseed<15; dseed+=3){ //loop over central seeding time bins
1301 iter+=1.;
1302 Int_t dseed = 5+Int_t(iter)*3;
69b55c55 1303 // Initialize seeding layers
6c94f330 1304 for (Int_t ilayer=0;ilayer<6;ilayer++){
69b55c55 1305 reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1306 xcl[ilayer] = reflayers[ilayer]->GetX();
1307 }
6c94f330 1308 //
1309 Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2])*0.5;
1310 AliTRDpropagationLayer& layer0=*reflayers[sLayer+0];
1311 AliTRDpropagationLayer& layer1=*reflayers[sLayer+1];
1312 AliTRDpropagationLayer& layer2=*reflayers[sLayer+2];
1313 AliTRDpropagationLayer& layer3=*reflayers[sLayer+3];
1314 //
1315 Int_t maxn3 = layer3;
1316 for (Int_t icl3=0;icl3<maxn3;icl3++){
69b55c55 1317 AliTRDcluster *cl3 = layer3[icl3];
1318 if (!cl3) continue;
1319 padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2()*12.);
1320 ycl[sLayer+3] = cl3->GetY();
1321 zcl[sLayer+3] = cl3->GetZ();
6c94f330 1322 Float_t yymin0 = ycl[sLayer+3] - 1- kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
1323 Float_t yymax0 = ycl[sLayer+3] + 1+ kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
1324 Int_t maxn0 = layer0; //
1325 for (Int_t icl0=layer0.Find(yymin0);icl0<maxn0;icl0++){
69b55c55 1326 AliTRDcluster *cl0 = layer0[icl0];
1327 if (!cl0) continue;
6c94f330 1328 if (cl3->IsUsed()&&cl0->IsUsed()) continue;
69b55c55 1329 ycl[sLayer+0] = cl0->GetY();
1330 zcl[sLayer+0] = cl0->GetZ();
6c94f330 1331 if ( ycl[sLayer+0]>yymax0) break;
1332 Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
1333 if (TMath::Abs(tanphi)>kMaxPhi) continue;
1334 Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
1335 if (TMath::Abs(tantheta)>kMaxTheta) continue;
1336 padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2()*12.);
1337 //
1338 // expected position in 1 layer
1339 Double_t y1exp = ycl[sLayer+0]+(tanphi) *(xcl[sLayer+1]-xcl[sLayer+0]);
1340 Double_t z1exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+1]-xcl[sLayer+0]);
1341 Float_t yymin1 = y1exp - kRoad0y-tanphi;
1342 Float_t yymax1 = y1exp + kRoad0y+tanphi;
1343 Int_t maxn1 = layer1; //
1344 //
1345 for (Int_t icl1=layer1.Find(yymin1);icl1<maxn1;icl1++){
69b55c55 1346 AliTRDcluster *cl1 = layer1[icl1];
1347 if (!cl1) continue;
1348 Int_t nusedCl = 0;
1349 if (cl3->IsUsed()) nusedCl++;
1350 if (cl0->IsUsed()) nusedCl++;
1351 if (cl1->IsUsed()) nusedCl++;
6c94f330 1352 if (nusedCl>1) continue;
69b55c55 1353 ycl[sLayer+1] = cl1->GetY();
1354 zcl[sLayer+1] = cl1->GetZ();
6c94f330 1355 if ( ycl[sLayer+1]>yymax1) break;
1356 if (TMath::Abs(ycl[sLayer+1]-y1exp)>kRoad0y+tanphi) continue;
1357 if (TMath::Abs(zcl[sLayer+1]-z1exp)>kRoad0z) continue;
69b55c55 1358 padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2()*12.);
6c94f330 1359 //
1360 Double_t y2exp = ycl[sLayer+0]+(tanphi) *(xcl[sLayer+2]-xcl[sLayer+0])+(ycl[sLayer+1]-y1exp);
1361 Double_t z2exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+2]-xcl[sLayer+0]);
1362 Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y, kRoad1z);
1363 if (index2<=0) continue;
1364 AliTRDcluster *cl2 = (AliTRDcluster*)GetCluster(index2);
1365 padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2()*12.);
1366 ycl[sLayer+2] = cl2->GetY();
1367 zcl[sLayer+2] = cl2->GetZ();
1368 if (TMath::Abs(cl2->GetZ()-z2exp)>kRoad0z) continue;
1369 //
69b55c55 1370 rieman.Reset();
1371 rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1372 rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1373 rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1374 rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1375 rieman.Update();
6c94f330 1376 //
1377 // reset fitter
1378 for (Int_t iLayer=0;iLayer<6;iLayer++){
69b55c55 1379 cseed[iLayer].Reset();
1380 }
6c94f330 1381 chi2Z =0.; chi2R=0.;
1382 for (Int_t iLayer=0;iLayer<4;iLayer++){
69b55c55 1383 cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
6c94f330 1384 chi2Z += (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer])*
1385 (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer]);
69b55c55 1386 cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1387 cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
6c94f330 1388 chi2R += (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer])*
1389 (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer]);
69b55c55 1390 cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1391 }
6c94f330 1392 if (TMath::Sqrt(chi2R)>1./iter) continue;
1393 if (TMath::Sqrt(chi2Z)>7./iter) continue;
1394 //
1395 //
1396 //
1397 Float_t minmax[2]={-100,100};
1398 for (Int_t iLayer=0;iLayer<4;iLayer++){
1399 Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer]*0.5+1 -cseed[sLayer+iLayer].fZref[0];
1400 if (max<minmax[1]) minmax[1]=max;
1401 Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer]*0.5-1 -cseed[sLayer+iLayer].fZref[0];
1402 if (min>minmax[0]) minmax[0]=min;
69b55c55 1403 }
1404 Bool_t isFake = kFALSE;
6c94f330 1405 if (cl0->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1406 if (cl1->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1407 if (cl2->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1408 if (AliTRDReconstructor::StreamLevel()>0){
1409 if ((!isFake) || (icl3%10)==0 ){ //debugging print
1410 TTreeSRedirector& cstream = *fDebugStreamer;
1411 cstream<<"Seeds0"<<
1412 "isFake="<<isFake<<
1413 "Cl0.="<<cl0<<
1414 "Cl1.="<<cl1<<
1415 "Cl2.="<<cl2<<
1416 "Cl3.="<<cl3<<
1417 "Xref="<<xref<<
1418 "X0="<<xcl[sLayer+0]<<
1419 "X1="<<xcl[sLayer+1]<<
1420 "X2="<<xcl[sLayer+2]<<
1421 "X3="<<xcl[sLayer+3]<<
1422 "Y2exp="<<y2exp<<
1423 "Z2exp="<<z2exp<<
1424 "Chi2R="<<chi2R<<
1425 "Chi2Z="<<chi2Z<<
1426 "Seed0.="<<&cseed[sLayer+0]<<
1427 "Seed1.="<<&cseed[sLayer+1]<<
1428 "Seed2.="<<&cseed[sLayer+2]<<
1429 "Seed3.="<<&cseed[sLayer+3]<<
1430 "Zmin="<<minmax[0]<<
1431 "Zmax="<<minmax[1]<<
1432 "\n";
d337ef8d 1433 }
69b55c55 1434 }
1435
6c94f330 1436 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1437 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1438 //<<<<<<<<<<<<<<<<<< FIT SEEDING PART <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1439 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
69b55c55 1440 cl[sLayer+0] = cl0;
1441 cl[sLayer+1] = cl1;
1442 cl[sLayer+2] = cl2;
1443 cl[sLayer+3] = cl3;
6c94f330 1444 Bool_t isOK=kTRUE;
1445 for (Int_t jLayer=0;jLayer<4;jLayer++){
1446 cseed[sLayer+jLayer].fTilt = hL[sLayer+jLayer];
69b55c55 1447 cseed[sLayer+jLayer].fPadLength = padlength[sLayer+jLayer];
6c94f330 1448 cseed[sLayer+jLayer].fX0 = xcl[sLayer+jLayer];
1449 for (Int_t iter=0; iter<2; iter++){
69b55c55 1450 //
6c94f330 1451 // in iteration 0 we try only one pad-row
1452 // if quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
69b55c55 1453 //
1454 AliTRDseed tseed = cseed[sLayer+jLayer];
6c94f330 1455 Float_t roadz = padlength[sLayer+jLayer]*0.5;
1456 if (iter>0) roadz = padlength[sLayer+jLayer];
1457 //
1458 Float_t quality =10000;
1459 for (Int_t iTime=2;iTime<20;iTime++){
1460 AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1461 Double_t dxlayer= layer.GetX()-xcl[sLayer+jLayer];
1462 Double_t zexp = cl[sLayer+jLayer]->GetZ() ;
1463 if (iter>0){
1464 // try 2 pad-rows in second iteration
1465 zexp = tseed.fZref[0]+ tseed.fZref[1]*dxlayer;
1466 if (zexp>cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()+padlength[sLayer+jLayer]*0.5;
1467 if (zexp<cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()-padlength[sLayer+jLayer]*0.5;
69b55c55 1468 }
6c94f330 1469 //
1470 Double_t yexp = tseed.fYref[0]+
1471 tseed.fYref[1]*dxlayer;
1472 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
1473 if (index<=0) continue;
1474 AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1475 //
69b55c55 1476 tseed.fIndexes[iTime] = index;
6c94f330 1477 tseed.fClusters[iTime] = cl; // register cluster
1478 tseed.fX[iTime] = dxlayer; // register cluster
1479 tseed.fY[iTime] = cl->GetY(); // register cluster
1480 tseed.fZ[iTime] = cl->GetZ(); // register cluster
69b55c55 1481 }
1482 tseed.Update();
6c94f330 1483 //count the number of clusters and distortions into quality
1484 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1485 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1486 TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
1487 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1488 if (iter==0 && tseed.IsOK()) {
69b55c55 1489 cseed[sLayer+jLayer] = tseed;
6c94f330 1490 quality = tquality;
1491 if (tquality<5) break;
69b55c55 1492 }
6c94f330 1493 if (tseed.IsOK() && tquality<quality)
1494 cseed[sLayer+jLayer] = tseed;
69b55c55 1495 }
6c94f330 1496 if (!cseed[sLayer+jLayer].IsOK()){
69b55c55 1497 isOK = kFALSE;
1498 break;
1499 }
1500 cseed[sLayer+jLayer].CookLabels();
1501 cseed[sLayer+jLayer].UpdateUsed();
6c94f330 1502 nusedCl+= cseed[sLayer+jLayer].fNUsed;
1503 if (nusedCl>25){
69b55c55 1504 isOK = kFALSE;
1505 break;
1506 }
1507 }
6c94f330 1508 //
69b55c55 1509 if (!isOK) continue;
6c94f330 1510 nclusters=0;
1511 for (Int_t iLayer=0;iLayer<4;iLayer++){
1512 if (cseed[sLayer+iLayer].IsOK()){
1513 nclusters+=cseed[sLayer+iLayer].fN2;
69b55c55 1514 }
1515 }
6c94f330 1516 //
1517 // iteration 0
69b55c55 1518 rieman.Reset();
6c94f330 1519 for (Int_t iLayer=0;iLayer<4;iLayer++){
1520 rieman.AddPoint(xcl[sLayer+iLayer],cseed[sLayer+iLayer].fYfitR[0],
1521 cseed[sLayer+iLayer].fZProb,1,10);
69b55c55 1522 }
1523 rieman.Update();
6c94f330 1524 //
1525 //
1526 chi2R =0; chi2Z=0;
1527 for (Int_t iLayer=0;iLayer<4;iLayer++){
69b55c55 1528 cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
6c94f330 1529 chi2R += (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0])*
1530 (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0]);
69b55c55 1531 cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1532 cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
6c94f330 1533 chi2Z += (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz)*
1534 (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz);
69b55c55 1535 cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1536 }
1537 Double_t curv = rieman.GetC();
1538 //
6c94f330 1539 // likelihoods
69b55c55 1540 //
6c94f330 1541 Double_t sumda =
1542 TMath::Abs(cseed[sLayer+0].fYfitR[1]- cseed[sLayer+0].fYref[1])+
1543 TMath::Abs(cseed[sLayer+1].fYfitR[1]- cseed[sLayer+1].fYref[1])+
1544 TMath::Abs(cseed[sLayer+2].fYfitR[1]- cseed[sLayer+2].fYref[1])+
1545 TMath::Abs(cseed[sLayer+3].fYfitR[1]- cseed[sLayer+3].fYref[1]);
1546 Double_t likea = TMath::Exp(-sumda*10.6);
1547 Double_t likechi2 = 0.0000000001;
1548 if (chi2R<0.5) likechi2+=TMath::Exp(-TMath::Sqrt(chi2R)*7.73);
1549 Double_t likechi2z = TMath::Exp(-chi2Z*0.088)/TMath::Exp(-chi2Z*0.019);
1550 Double_t likeN = TMath::Exp(-(72-nclusters)*0.19);
1551 Double_t like = likea*likechi2*likechi2z*likeN;
1552 //
1553 Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fYref[1]-130*curv)*1.9);
1554 Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fZref[1]-
1555 cseed[sLayer+0].fZref[0]/xcl[sLayer+0])*5.9);
1556 Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
69b55c55 1557
6c94f330 1558 seedquality[registered] = like;
1559 seedlayer[registered] = sLayer;
1560 if (TMath::Log(0.000000000000001+like)<-15) continue;
69b55c55 1561 AliTRDseed seedb[6];
6c94f330 1562 for (Int_t iLayer=0;iLayer<6;iLayer++){
69b55c55 1563 seedb[iLayer] = cseed[iLayer];
1564 }
ad670fba 1565 //
6c94f330 1566 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1567 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1568 //<<<<<<<<<<<<<<< FULL TRACK FIT PART <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1569 //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
69b55c55 1570 //
6c94f330 1571 Int_t nlayers = 0;
1572 Int_t nusedf = 0;
1573 Int_t findable = 0;
69b55c55 1574 //
6c94f330 1575 // add new layers - avoid long extrapolation
69b55c55 1576 //
6c94f330 1577 Int_t tLayer[2]={0,0};
1578 if (sLayer==2) {tLayer[0]=1; tLayer[1]=0;}
1579 if (sLayer==1) {tLayer[0]=5; tLayer[1]=0;}
1580 if (sLayer==0) {tLayer[0]=4; tLayer[1]=5;}
69b55c55 1581 //
6c94f330 1582 for (Int_t iLayer=0;iLayer<2;iLayer++){
1583 Int_t jLayer = tLayer[iLayer]; // set tracking layer
69b55c55 1584 cseed[jLayer].Reset();
6c94f330 1585 cseed[jLayer].fTilt = hL[jLayer];
69b55c55 1586 cseed[jLayer].fPadLength = padlength[jLayer];
6c94f330 1587 cseed[jLayer].fX0 = xcl[jLayer];
1588 // get pad length and rough cluster
1589 Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].fYref[0],
1590 cseed[jLayer].fZref[0],kRoad2y,kRoad2z);
1591 if (indexdummy<=0) continue;
1592 AliTRDcluster *cldummy = (AliTRDcluster*)GetCluster(indexdummy);
1593 padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2()*12.);
69b55c55 1594 }
6c94f330 1595 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1596 //
1597 for (Int_t iLayer=0;iLayer<2;iLayer++){
1598 Int_t jLayer = tLayer[iLayer]; // set tracking layer
1599 if ( (jLayer==0) && !(cseed[1].IsOK())) continue; // break not allowed
1600 if ( (jLayer==5) && !(cseed[4].IsOK())) continue; // break not allowed
69b55c55 1601 Float_t zexp = cseed[jLayer].fZref[0];
6c94f330 1602 Double_t zroad = padlength[jLayer]*0.5+1.;
1603 //
1604 //
1605 for (Int_t iter=0;iter<2;iter++){
69b55c55 1606 AliTRDseed tseed = cseed[jLayer];
6c94f330 1607 Float_t quality = 10000;
1608 for (Int_t iTime=2;iTime<20;iTime++){
1609 AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1610 Double_t dxlayer = layer.GetX()-xcl[jLayer];
1611 Double_t yexp = tseed.fYref[0]+tseed.fYref[1]*dxlayer;
1612 Float_t yroad = kRoad1y;
1613 Int_t index = layer.FindNearestCluster(yexp,zexp, yroad, zroad);
1614 if (index<=0) continue;
1615 AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1616 //
69b55c55 1617 tseed.fIndexes[iTime] = index;
6c94f330 1618 tseed.fClusters[iTime] = cl; // register cluster
1619 tseed.fX[iTime] = dxlayer; // register cluster
1620 tseed.fY[iTime] = cl->GetY(); // register cluster
1621 tseed.fZ[iTime] = cl->GetZ(); // register cluster
69b55c55 1622 }
1623 tseed.Update();
6c94f330 1624 if (tseed.IsOK()){
1625 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1626 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1627 TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
1628 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1629 //
1630 if (tquality<quality){
1631 cseed[jLayer]=tseed;
1632 quality = tquality;
69b55c55 1633 }
1634 }
6c94f330 1635 zroad*=2.;
69b55c55 1636 }
6c94f330 1637 if ( cseed[jLayer].IsOK()){
69b55c55 1638 cseed[jLayer].CookLabels();
1639 cseed[jLayer].UpdateUsed();
6c94f330 1640 nusedf+= cseed[jLayer].fNUsed;
1641 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
69b55c55 1642 }
1643 }
6c94f330 1644 //
1645 //
1646 // make copy
69b55c55 1647 AliTRDseed bseed[6];
6c94f330 1648 for (Int_t jLayer=0;jLayer<6;jLayer++){
69b55c55 1649 bseed[jLayer] = cseed[jLayer];
1650 }
6c94f330 1651 Float_t lastquality = 10000;
1652 Float_t lastchi2 = 10000;
1653 Float_t chi2 = 1000;
69b55c55 1654
6c94f330 1655 //
1656 for (Int_t iter =0; iter<4;iter++){
69b55c55 1657 //
6c94f330 1658 // sort tracklets according "quality", try to "improve" 4 worst
69b55c55 1659 //
6c94f330 1660 Float_t sumquality = 0;
69b55c55 1661 Float_t squality[6];
1662 Int_t sortindexes[6];
6c94f330 1663 for (Int_t jLayer=0;jLayer<6;jLayer++){
1664 if (bseed[jLayer].IsOK()){
69b55c55 1665 AliTRDseed &tseed = bseed[jLayer];
6c94f330 1666 Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1667 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1668 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1669 TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
1670 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1671 squality[jLayer] = tquality;
ad670fba 1672 }
6c94f330 1673 else squality[jLayer]=-1;
1674 sumquality +=squality[jLayer];
69b55c55 1675 }
1676
6c94f330 1677 if (sumquality>=lastquality || chi2>lastchi2) break;
69b55c55 1678 lastquality = sumquality;
1679 lastchi2 = chi2;
6c94f330 1680 if (iter>0){
1681 for (Int_t jLayer=0;jLayer<6;jLayer++){
69b55c55 1682 cseed[jLayer] = bseed[jLayer];
1683 }
1684 }
1685 TMath::Sort(6,squality,sortindexes,kFALSE);
6c94f330 1686 //
1687 //
1688 for (Int_t jLayer=5;jLayer>1;jLayer--){
1689 Int_t bLayer = sortindexes[jLayer];
1690 AliTRDseed tseed = bseed[bLayer];
69b55c55 1691 for (Int_t iTime=2;iTime<20;iTime++){
6c94f330 1692 AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1693 Double_t dxlayer= layer.GetX()-xcl[bLayer];
1694 //
1695 Double_t zexp = tseed.fZref[0];
1696 Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1697 //
1698 Float_t roadz = padlength[bLayer]+1;
1699 if (TMath::Abs(tseed.fZProb-zexp)> padlength[bLayer]*0.5) {roadz = padlength[bLayer]*0.5;}
1700 if (tseed.fZfit[1]*tseed.fZref[1]<0) {roadz = padlength[bLayer]*0.5;}
1701 if (TMath::Abs(tseed.fZProb-zexp)<0.1*padlength[bLayer]) {
1702 zexp = tseed.fZProb;
1703 roadz = padlength[bLayer]*0.5;
69b55c55 1704 }
6c94f330 1705 //
1706 Double_t yexp = tseed.fYref[0]+
1707 tseed.fYref[1]*dxlayer-zcor;
1708 Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
1709 if (index<=0) continue;
1710 AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1711 //
69b55c55 1712 tseed.fIndexes[iTime] = index;
6c94f330 1713 tseed.fClusters[iTime] = cl; // register cluster
1714 tseed.fX[iTime] = dxlayer; // register cluster
1715 tseed.fY[iTime] = cl->GetY(); // register cluster
1716 tseed.fZ[iTime] = cl->GetZ(); // register cluster
69b55c55 1717 }
1718 tseed.Update();
c6f438c0 1719 if (tseed.IsOK()) {
6c94f330 1720 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1721 Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1722 //
1723 Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1724 TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
1725 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1726 //
1727 if (tquality<squality[bLayer])
69b55c55 1728 bseed[bLayer] = tseed;
1729 }
1730 }
6c94f330 1731 chi2 = AliTRDseed::FitRiemanTilt(bseed, kTRUE);
69b55c55 1732 }
6c94f330 1733 //
1734 //
1735 //
1736 nclusters = 0;
1737 nlayers = 0;
1738 findable = 0;
1739 for (Int_t iLayer=0;iLayer<6;iLayer++) {
1740 if (TMath::Abs(cseed[iLayer].fYref[0]/cseed[iLayer].fX0)<0.15)
69b55c55 1741 findable++;
6c94f330 1742 if (cseed[iLayer].IsOK()){
1743 nclusters+=cseed[iLayer].fN2;
69b55c55 1744 nlayers++;
1745 }
1746 }
6c94f330 1747 if (nlayers<3) continue;
69b55c55 1748 rieman.Reset();
6c94f330 1749 for (Int_t iLayer=0;iLayer<6;iLayer++){
1750 if (cseed[iLayer].IsOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
1751 cseed[iLayer].fZProb,1,10);
69b55c55 1752 }
1753 rieman.Update();
6c94f330 1754 //
1755 chi2RF =0;
1756 chi2ZF =0;
1757 for (Int_t iLayer=0;iLayer<6;iLayer++){
1758 if (cseed[iLayer].IsOK()){
69b55c55 1759 cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
6c94f330 1760 chi2RF += (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0])*
1761 (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0]);
69b55c55 1762 cseed[iLayer].fYref[1] = rieman.GetDYat(xcl[iLayer]);
1763 cseed[iLayer].fZref[0] = rieman.GetZat(xcl[iLayer]);
6c94f330 1764 chi2ZF += (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz)*
1765 (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz);
69b55c55 1766 cseed[iLayer].fZref[1] = rieman.GetDZat(xcl[iLayer]);
1767 }
1768 }
6c94f330 1769 chi2RF/=TMath::Max((nlayers-3.),1.);
1770 chi2ZF/=TMath::Max((nlayers-3.),1.);
69b55c55 1771 curv = rieman.GetC();
6c94f330 1772
1773 //
69b55c55 1774
6c94f330 1775 Double_t xref2 = (xcl[2]+xcl[3])*0.5; // middle of the chamber
1776 Double_t dzmf = rieman.GetDZat(xref2);
1777 Double_t zmf = rieman.GetZat(xref2);
69b55c55 1778 //
6c94f330 1779 // fit hyperplane
69b55c55 1780 //
6c94f330 1781 Int_t npointsT =0;
69b55c55 1782 fitterTC.ClearPoints();
1783 fitterT2.ClearPoints();
1784 rieman2.Reset();
6c94f330 1785 for (Int_t iLayer=0; iLayer<6;iLayer++){
c6f438c0 1786 if (!cseed[iLayer].IsOK()) continue;
6c94f330 1787 for (Int_t itime=0;itime<25;itime++){
69b55c55 1788 if (!cseed[iLayer].fUsable[itime]) continue;
6c94f330 1789 Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2; // x relative to the midle chamber
1790 Double_t y = cseed[iLayer].fY[itime];
1791 Double_t z = cseed[iLayer].fZ[itime];
69b55c55 1792 // ExB correction to the correction
6c94f330 1793 // tilted rieman
1794 //
69b55c55 1795 Double_t uvt[6];
6c94f330 1796 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
1797 //
1798 Double_t t = 1./(x2*x2+y*y);
1799 uvt[1] = t; // t
1800 uvt[0] = 2.*x2*uvt[1]; // u
1801 //
1802 uvt[2] = 2.0*hL[iLayer]*uvt[1];
1803 uvt[3] = 2.0*hL[iLayer]*x*uvt[1];
1804 uvt[4] = 2.0*(y+hL[iLayer]*z)*uvt[1];
1805 //
1806 Double_t error = 2*0.2*uvt[1];
69b55c55 1807 fitterT2.AddPoint(uvt,uvt[4],error);
1808 //
6c94f330 1809 // constrained rieman
69b55c55 1810 //
6c94f330 1811 z =cseed[iLayer].fZ[itime];
1812 uvt[0] = 2.*x2*t; // u
1813 uvt[1] = 2*hL[iLayer]*x2*uvt[1];
1814 uvt[2] = 2*(y+hL[iLayer]*(z-GetZ()))*t;
69b55c55 1815 fitterTC.AddPoint(uvt,uvt[2],error);
6c94f330 1816 //
69b55c55 1817 rieman2.AddPoint(x2,y,z,1,10);
1818 npointsT++;
1819 }
1820 }
1821 rieman2.Update();
1822 fitterTC.Eval();
1823 fitterT2.Eval();
1824 Double_t rpolz0 = fitterT2.GetParameter(3);
1825 Double_t rpolz1 = fitterT2.GetParameter(4);
1826 //
6c94f330 1827 // linear fitter - not possible to make boundaries
69b55c55 1828 // non accept non possible z and dzdx combination
1829 //
6c94f330 1830 Bool_t acceptablez =kTRUE;
1831 for (Int_t iLayer=0; iLayer<6;iLayer++){
1832 if (cseed[iLayer].IsOK()){
1833 Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
1834 if (TMath::Abs(cseed[iLayer].fZProb-zT2)>padlength[iLayer]*0.5+1)
69b55c55 1835 acceptablez = kFALSE;
1836 }
1837 }
6c94f330 1838 if (!acceptablez){
69b55c55 1839 fitterT2.FixParameter(3,zmf);
1840 fitterT2.FixParameter(4,dzmf);
1841 fitterT2.Eval();
1842 fitterT2.ReleaseParameter(3);
1843 fitterT2.ReleaseParameter(4);
1844 rpolz0 = fitterT2.GetParameter(3);
1845 rpolz1 = fitterT2.GetParameter(4);
1846 }
6c94f330 1847 //
1848 Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
1849 Double_t chi2TC = fitterTC.GetChisquare()/Float_t(npointsT);
1850 //
69b55c55 1851 Double_t polz1c = fitterTC.GetParameter(2);
6c94f330 1852 Double_t polz0c = polz1c*xref2;
1853 //
1854 Double_t aC = fitterTC.GetParameter(0);
1855 Double_t bC = fitterTC.GetParameter(1);
1856 Double_t cC = aC/TMath::Sqrt(bC*bC+1.); // curvature
1857 //
1858 Double_t aR = fitterT2.GetParameter(0);
1859 Double_t bR = fitterT2.GetParameter(1);
1860 Double_t dR = fitterT2.GetParameter(2);
1861 Double_t cR = 1+bR*bR-dR*aR;
1862 Double_t dca = 0.;
1863 if (cR>0){
1864 dca = -dR/(TMath::Sqrt(1+bR*bR-dR*aR)+TMath::Sqrt(1+bR*bR));
1865 cR = aR/TMath::Sqrt(cR);
69b55c55 1866 }
6c94f330 1867 //
1868 Double_t chi2ZT2=0, chi2ZTC=0;
1869 for (Int_t iLayer=0; iLayer<6;iLayer++){
1870 if (cseed[iLayer].IsOK()){
1871 Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
1872 Double_t zTC = polz0c+polz1c*(xcl[iLayer] - xref2);
1873 chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz-zT2);
1874 chi2ZTC += TMath::Abs(cseed[iLayer].fMeanz-zTC);
69b55c55 1875 }
1876 }
6c94f330 1877 chi2ZT2/=TMath::Max((nlayers-3.),1.);
1878 chi2ZTC/=TMath::Max((nlayers-3.),1.);
1879 //
1880 //
1881 //
1882 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1883 Float_t sumdaf = 0;
1884 for (Int_t iLayer=0;iLayer<6;iLayer++){
1885 if (cseed[iLayer].IsOK())
1886 sumdaf += TMath::Abs((cseed[iLayer].fYfit[1]-cseed[iLayer].fYref[1])/cseed[iLayer].fSigmaY2);
69b55c55 1887 }
6c94f330 1888 sumdaf /= Float_t (nlayers-2.);
69b55c55 1889 //
6c94f330 1890 // likelihoods for full track
69b55c55 1891 //
6c94f330 1892 Double_t likezf = TMath::Exp(-chi2ZF*0.14);
1893 Double_t likechi2C = TMath::Exp(-chi2TC*0.677);
1894 Double_t likechi2TR = TMath::Exp(-chi2TR*0.78);
1895 Double_t likeaf = TMath::Exp(-sumdaf*3.23);
1896 seedquality2[registered] = likezf*likechi2TR*likeaf;
69b55c55 1897// Bool_t isGold = kFALSE;
1898//
6c94f330 1899// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
1900// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
69b55c55 1901// if (isGold &&nusedf<10){
1902// for (Int_t jLayer=0;jLayer<6;jLayer++){
6c94f330 1903// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
69b55c55 1904// seed[index][jLayer].UseClusters(); //sign gold
1905// }
1906// }
6c94f330 1907 //
1908 //
1909 //
1910 Int_t index0=0;
1911 if (!cseed[0].IsOK()){
69b55c55 1912 index0 = 1;
6c94f330 1913 if (!cseed[1].IsOK()) index0 = 2;
69b55c55 1914 }
1915 seedparams[registered][0] = cseed[index0].fX0;
1916 seedparams[registered][1] = cseed[index0].fYref[0];
1917 seedparams[registered][2] = cseed[index0].fZref[0];
c6f438c0 1918 seedparams[registered][5] = cR;
6c94f330 1919 seedparams[registered][3] = cseed[index0].fX0*cR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
1920 seedparams[registered][4] = cseed[index0].fZref[1]/
1921 TMath::Sqrt(1+cseed[index0].fYref[1]*cseed[index0].fYref[1]);
69b55c55 1922 seedparams[registered][6] = ns;
6c94f330 1923 //
1924 //
1925 Int_t labels[12], outlab[24];
1926 Int_t nlab=0;
1927 for (Int_t iLayer=0;iLayer<6;iLayer++){
c6f438c0 1928 if (!cseed[iLayer].IsOK()) continue;
6c94f330 1929 if (cseed[iLayer].fLabels[0]>=0) {
69b55c55 1930 labels[nlab] = cseed[iLayer].fLabels[0];
1931 nlab++;
1932 }
6c94f330 1933 if (cseed[iLayer].fLabels[1]>=0) {
69b55c55 1934 labels[nlab] = cseed[iLayer].fLabels[1];
1935 nlab++;
1936 }
1937 }
1938 Freq(nlab,labels,outlab,kFALSE);
6c94f330 1939 Int_t label = outlab[0];
1940 Int_t frequency = outlab[1];
1941 for (Int_t iLayer=0;iLayer<6;iLayer++){
69b55c55 1942 cseed[iLayer].fFreq = frequency;
c6f438c0 1943 cseed[iLayer].fC = cR;
6c94f330 1944 cseed[iLayer].fCC = cC;
69b55c55 1945 cseed[iLayer].fChi2 = chi2TR;
1946 cseed[iLayer].fChi2Z = chi2ZF;
1947 }
6c94f330 1948 //
1949 if (1||(!isFake)){ //debugging print
69b55c55 1950 Float_t zvertex = GetZ();
6c94f330 1951 TTreeSRedirector& cstream = *fDebugStreamer;
1952 if (AliTRDReconstructor::StreamLevel()>0)
1953 cstream<<"Seeds1"<<
1954 "isFake="<<isFake<<
1955 "Vertex="<<zvertex<<
1956 "Rieman2.="<<&rieman2<<
1957 "Rieman.="<<&rieman<<
1958 "Xref="<<xref<<
1959 "X0="<<xcl[0]<<
1960 "X1="<<xcl[1]<<
1961 "X2="<<xcl[2]<<
1962 "X3="<<xcl[3]<<
1963 "X4="<<xcl[4]<<
1964 "X5="<<xcl[5]<<
1965 "Chi2R="<<chi2R<<
1966 "Chi2Z="<<chi2Z<<
1967 "Chi2RF="<<chi2RF<< //chi2 of trackletes on full track
1968 "Chi2ZF="<<chi2ZF<< //chi2 z on tracklets on full track
1969 "Chi2ZT2="<<chi2ZT2<< //chi2 z on tracklets on full track - rieman tilt
1970 "Chi2ZTC="<<chi2ZTC<< //chi2 z on tracklets on full track - rieman tilt const
1971 //
1972 "Chi2TR="<<chi2TR<< //chi2 without vertex constrain
1973 "Chi2TC="<<chi2TC<< //chi2 with vertex constrain
1974 "C="<<curv<< // non constrained - no tilt correction
1975 "DR="<<dR<< // DR parameter - tilt correction
1976 "DCA="<<dca<< // DCA - tilt correction
1977 "CR="<<cR<< // non constrained curvature - tilt correction
1978 "CC="<<cC<< // constrained curvature
1979 "Polz0="<<polz0c<<
1980 "Polz1="<<polz1c<<
1981 "RPolz0="<<rpolz0<<
1982 "RPolz1="<<rpolz1<<
1983 "Ncl="<<nclusters<<
1984 "Nlayers="<<nlayers<<
1985 "NUsedS="<<nusedCl<<
1986 "NUsed="<<nusedf<<
1987 "Findable="<<findable<<
1988 "Like="<<like<<
1989 "LikePrim="<<likePrim<<
1990 "Likechi2C="<<likechi2C<<
1991 "Likechi2TR="<<likechi2TR<<
1992 "Likezf="<<likezf<<
1993 "LikeF="<<seedquality2[registered]<<
1994 "S0.="<<&cseed[0]<<
1995 "S1.="<<&cseed[1]<<
1996 "S2.="<<&cseed[2]<<
1997 "S3.="<<&cseed[3]<<
1998 "S4.="<<&cseed[4]<<
1999 "S5.="<<&cseed[5]<<
2000 "SB0.="<<&seedb[0]<<
2001 "SB1.="<<&seedb[1]<<
2002 "SB2.="<<&seedb[2]<<
2003 "SB3.="<<&seedb[3]<<
2004 "SB4.="<<&seedb[4]<<
2005 "SB5.="<<&seedb[5]<<
2006 "Label="<<label<<
2007 "Freq="<<frequency<<
2008 "sLayer="<<sLayer<<
2009 "\n";
69b55c55 2010 }
6c94f330 2011 if (registered<kMaxSeed-1) {
69b55c55 2012 registered++;
2013 cseed = seed[registered];
2014 }
6c94f330 2015 }// end of loop over layer 1
2016 } // end of loop over layer 0
2017 } // end of loop over layer 3
2018 } // end of loop over seeding time bins
69b55c55 2019 //
6c94f330 2020 // choos best
69b55c55 2021 //
2022 TMath::Sort(registered,seedquality2,sort,kTRUE);
c6f438c0 2023 Bool_t signedseed[kMaxSeed];
6c94f330 2024 for (Int_t i=0;i<registered;i++){
2025 signedseed[i]= kFALSE;
69b55c55 2026 }
6c94f330 2027 for (Int_t iter=0; iter<5; iter++){
2028 for (Int_t iseed=0;iseed<registered;iseed++){
69b55c55 2029 Int_t index = sort[iseed];
2030 if (signedseed[index]) continue;
2031 Int_t labelsall[1000];
6c94f330 2032 Int_t nlabelsall=0;
2033 Int_t naccepted=0;;
2034 Int_t sLayer = seedlayer[index];
2035 Int_t ncl = 0;
2036 Int_t nused = 0;
2037 Int_t nlayers =0;
69b55c55 2038 Int_t findable = 0;
6c94f330 2039 for (Int_t jLayer=0;jLayer<6;jLayer++){
2040 if (TMath::Abs(seed[index][jLayer].fYref[0]/xcl[jLayer])<0.15)
69b55c55 2041 findable++;
6c94f330 2042 if (seed[index][jLayer].IsOK()){
69b55c55 2043 seed[index][jLayer].UpdateUsed();
6c94f330 2044 ncl +=seed[index][jLayer].fN2;
2045 nused +=seed[index][jLayer].fNUsed;
69b55c55 2046 nlayers++;
6c94f330 2047 //cooking label
2048 for (Int_t itime=0;itime<25;itime++){
2049 if (seed[index][jLayer].fUsable[itime]){
69b55c55 2050 naccepted++;
6c94f330 2051 for (Int_t ilab=0;ilab<3;ilab++){
69b55c55 2052 Int_t tindex = seed[index][jLayer].fClusters[itime]->GetLabel(ilab);
6c94f330 2053 if (tindex>=0){
69b55c55 2054 labelsall[nlabelsall] = tindex;
2055 nlabelsall++;
2056 }
2057 }
2058 }
2059 }
2060 }
2061 }
6c94f330 2062 //
2063 if (nused>30) continue;
2064 //
2065 if (iter==0){
2066 if (nlayers<6) continue;
2067 if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue; // gold
69b55c55 2068 }
6c94f330 2069 //
2070 if (iter==1){
2071 if (nlayers<findable) continue;
2072 if (TMath::Log(0.000000001+seedquality2[index])<-4.) continue; //
7ad19338 2073 }
6c94f330 2074 //
2075 //
2076 if (iter==2){
2077 if (nlayers==findable || nlayers==6) continue;
2078 if (TMath::Log(0.000000001+seedquality2[index])<-6.) continue;
69b55c55 2079 }
6c94f330 2080 //
2081 if (iter==3){
2082 if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue;
69b55c55 2083 }
6c94f330 2084 //
2085 if (iter==4){
2086 if (TMath::Log(0.000000001+seedquality2[index])-nused/(nlayers-3.)<-15.) continue;
69b55c55 2087 }
6c94f330 2088 //
69b55c55 2089 signedseed[index] = kTRUE;
6c94f330 2090 //
2091 Int_t labels[1000], outlab[1000];
2092 Int_t nlab=0;
2093 for (Int_t iLayer=0;iLayer<6;iLayer++){
2094 if (seed[index][iLayer].IsOK()){
2095 if (seed[index][iLayer].fLabels[0]>=0) {
69b55c55 2096 labels[nlab] = seed[index][iLayer].fLabels[0];
2097 nlab++;
2098 }
6c94f330 2099 if (seed[index][iLayer].fLabels[1]>=0) {
69b55c55 2100 labels[nlab] = seed[index][iLayer].fLabels[1];
2101 nlab++;
2102 }
2103 }
7ad19338 2104 }
69b55c55 2105 Freq(nlab,labels,outlab,kFALSE);
6c94f330 2106 Int_t label = outlab[0];
2107 Int_t frequency = outlab[1];
69b55c55 2108 Freq(nlabelsall,labelsall,outlab,kFALSE);
6c94f330 2109 Int_t label1 = outlab[0];
2110 Int_t label2 = outlab[2];
2111 Float_t fakeratio = (naccepted-outlab[1])/Float_t(naccepted);
2112 Float_t ratio = Float_t(nused)/Float_t(ncl);
2113 if (ratio<0.25){
2114 for (Int_t jLayer=0;jLayer<6;jLayer++){
2115 if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.2 )
2116 seed[index][jLayer].UseClusters(); //sign gold
69b55c55 2117 }
7ad19338 2118 }
6c94f330 2119 //
69b55c55 2120 Int_t eventNr = esd->GetEventNumber();
6c94f330 2121 TTreeSRedirector& cstream = *fDebugStreamer;
69b55c55 2122 //
6c94f330 2123 // register seed
69b55c55 2124 //
6c94f330 2125 AliTRDtrack * track = RegisterSeed(seed[index],seedparams[index]);
2126 AliTRDtrack dummy;
2127 if (!track) track=&dummy;
2128 else{
69b55c55 2129 AliESDtrack esdtrack;
6c94f330 2130 esdtrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
69b55c55 2131 esdtrack.SetLabel(label);
2132 esd->AddTrack(&esdtrack);
6c94f330 2133 TTreeSRedirector& cstream = *fDebugStreamer;
2134 if (AliTRDReconstructor::StreamLevel()>0)
2135 cstream<<"Tracks"<<
2136 "EventNr="<<eventNr<<
2137 "ESD.="<<&esdtrack<<
2138 "trd.="<<track<<
2139 "trdback.="<<track<<
2140 "\n";
7ad19338 2141 }
6c94f330 2142 if (AliTRDReconstructor::StreamLevel()>0)
2143 cstream<<"Seeds2"<<
2144 "Iter="<<iter<<
2145 "Track.="<<track<<
2146 "Like="<<seedquality[index]<<
2147 "LikeF="<<seedquality2[index]<<
2148 "S0.="<<&seed[index][0]<<
2149 "S1.="<<&seed[index][1]<<
2150 "S2.="<<&seed[index][2]<<
2151 "S3.="<<&seed[index][3]<<
2152 "S4.="<<&seed[index][4]<<
2153 "S5.="<<&seed[index][5]<<
2154 "Label="<<label<<
2155 "Label1="<<label1<<
2156 "Label2="<<label2<<
2157 "FakeRatio="<<fakeratio<<
2158 "Freq="<<frequency<<
2159 "Ncl="<<ncl<<
2160 "Nlayers="<<nlayers<<
2161 "Findable="<<findable<<
2162 "NUsed="<<nused<<
2163 "sLayer="<<sLayer<<
2164 "EventNr="<<eventNr<<
2165 "\n";
7ad19338 2166 }
2167 }
6c94f330 2168 } // end of loop over sectors
75bd7f81 2169
69b55c55 2170 delete [] pseed;
75bd7f81 2171
69b55c55 2172}
2173
5443e65e 2174//_____________________________________________________________________________
b7a0917f 2175Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const
5443e65e 2176{
2177 //
a819a5f7 2178 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2179 // from the file. The names of the cluster tree and branches
2180 // should match the ones used in AliTRDclusterizer::WriteClusters()
2181 //
75bd7f81 2182
6c94f330 2183 Int_t nsize = Int_t(ClusterTree->GetTotBytes()/(sizeof(AliTRDcluster)));
2184 TObjArray *clusterArray = new TObjArray(nsize+1000);
5443e65e 2185
c630aafd 2186 TBranch *branch=ClusterTree->GetBranch("TRDcluster");
2187 if (!branch) {
6c94f330 2188 Error("ReadClusters","Can't get the branch !");
c630aafd 2189 return 1;
2190 }
029cd327 2191 branch->SetAddress(&clusterArray);
5443e65e 2192
6c94f330 2193 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
2194 // printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
2195
a819a5f7 2196 // Loop through all entries in the tree
6c94f330 2197 Int_t nbytes = 0;
a819a5f7 2198 AliTRDcluster *c = 0;
6c94f330 2199 // printf("\n");
a819a5f7 2200 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2201
2202 // Import the tree
5443e65e 2203 nbytes += ClusterTree->GetEvent(iEntry);
2204
a819a5f7 2205 // Get the number of points in the detector
029cd327 2206 Int_t nCluster = clusterArray->GetEntriesFast();
6c94f330 2207// printf("\r Read %d clusters from entry %d", nCluster, iEntry);
5443e65e 2208
a819a5f7 2209 // Loop through all TRD digits
2210 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
6c94f330 2211 c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
4f1c04d3 2212 AliTRDcluster *co = c;
a819a5f7 2213 array->AddLast(co);
6c94f330 2214 // delete clusterArray->RemoveAt(iCluster);
4f1c04d3 2215 clusterArray->RemoveAt(iCluster);
a819a5f7 2216 }
2217 }
6c94f330 2218// cout<<"Allocated"<<nsize<<"\tLoaded"<<array->GetEntriesFast()<<"\n";
a819a5f7 2219
029cd327 2220 delete clusterArray;
5443e65e 2221
c630aafd 2222 return 0;
75bd7f81 2223
a819a5f7 2224}
2225
75bd7f81 2226//_____________________________________________________________________________
6c94f330 2227Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint& p) const
3551db50 2228{
2229 //
2230 // Get track space point with index i
2231 // Origin: C.Cheshkov
2232 //
2233
6c94f330 2234 AliTRDcluster *cl = (AliTRDcluster*)fClusters->UncheckedAt(index);
2235 Int_t idet = cl->GetDetector();
2236 Int_t isector = fGeom->GetSector(idet);
2237 Int_t ichamber= fGeom->GetChamber(idet);
2238 Int_t iplan = fGeom->GetPlane(idet);
3551db50 2239 Double_t local[3];
6c94f330 2240 local[0]=GetX(isector,iplan,cl->GetLocalTimeBin());
2241 local[1]=cl->GetY();
2242 local[2]=cl->GetZ();
3551db50 2243 Double_t global[3];
2244 fGeom->RotateBack(idet,local,global);
2245 p.SetXYZ(global[0],global[1],global[2]);
2246 AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
2247 switch (iplan) {
2248 case 0:
2249 iLayer = AliAlignObj::kTRD1;
2250 break;
2251 case 1:
2252 iLayer = AliAlignObj::kTRD2;
2253 break;
2254 case 2:
2255 iLayer = AliAlignObj::kTRD3;
2256 break;
2257 case 3:
2258 iLayer = AliAlignObj::kTRD4;
2259 break;
2260 case 4:
2261 iLayer = AliAlignObj::kTRD5;
2262 break;
2263 case 5:
2264 iLayer = AliAlignObj::kTRD6;
2265 break;
2266 };
6c94f330 2267 Int_t modId = isector*fGeom->Ncham()+ichamber;
3551db50 2268 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
2269 p.SetVolumeID(volid);
2270
2271 return kTRUE;
2272
2273}
2274
75bd7f81 2275//_____________________________________________________________________________
029cd327 2276void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
2277{
2278 //
2279 // This cooks a label. Mmmmh, smells good...
2280 //
46d29e70 2281
6c94f330 2282 Int_t label=123456789, index, i, j;
2283 Int_t ncl=pt->GetNumberOfClusters();
2284 const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
5443e65e 2285
029cd327 2286 Bool_t labelAdded;
46d29e70 2287
6c94f330 2288 // Int_t s[kRange][2];
2289 Int_t **s = new Int_t* [kRange];
2290 for (i=0; i<kRange; i++) {
d1b06c24 2291 s[i] = new Int_t[2];
2292 }
6c94f330 2293 for (i=0; i<kRange; i++) {
2294 s[i][0]=-1;
2295 s[i][1]=0;
46d29e70 2296 }
2297
6c94f330 2298 Int_t t0,t1,t2;
2299 for (i=0; i<ncl; i++) {
2300 index=pt->GetClusterIndex(i);
2301 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2302 t0=c->GetLabel(0);
2303 t1=c->GetLabel(1);
2304 t2=c->GetLabel(2);
46d29e70 2305 }
2306
6c94f330 2307 for (i=0; i<ncl; i++) {
2308 index=pt->GetClusterIndex(i);
2309 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2310 for (Int_t k=0; k<3; k++) {
2311 label=c->GetLabel(k);
2312 labelAdded=kFALSE; j=0;
46d29e70 2313 if (label >= 0) {
6c94f330 2314 while ( (!labelAdded) && ( j < kRange ) ) {
2315 if (s[j][0]==label || s[j][1]==0) {
2316 s[j][0]=label;
2317 s[j][1]=s[j][1]+1;
2318 labelAdded=kTRUE;
a9814c09 2319 }
2320 j++;
2321 }
46d29e70 2322 }
2323 }
2324 }
2325
6c94f330 2326 Int_t max=0;
2327 label = -123456789;
46d29e70 2328
6c94f330 2329 for (i=0; i<kRange; i++) {
2330 if (s[i][1]>max) {
2331 max=s[i][1]; label=s[i][0];
46d29e70 2332 }
2333 }
5443e65e 2334
6c94f330 2335 for (i=0; i<kRange; i++) {
2336 delete []s[i];
5443e65e 2337 }
2338
6c94f330 2339 delete []s;
5443e65e 2340
6c94f330 2341 if ((1.- Float_t(max)/ncl) > wrong) label=-label;
5443e65e 2342
2343 pt->SetLabel(label);
2344
46d29e70 2345}
2346
75bd7f81 2347//_____________________________________________________________________________
029cd327 2348void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
2349{
2350 //
2351 // Use clusters, but don't abuse them!
2352 //
75bd7f81 2353
6c94f330 2354 const Float_t kmaxchi2 =18;
2355 const Float_t kmincl =10;
2356 AliTRDtrack * track = (AliTRDtrack*)t;
2357 //
2358 Int_t ncl=t->GetNumberOfClusters();
2359 for (Int_t i=from; i<ncl; i++) {
5443e65e 2360 Int_t index = t->GetClusterIndex(i);
6c94f330 2361 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2362 //
69b55c55 2363 Int_t iplane = fGeom->GetPlane(c->GetDetector());
6c94f330 2364 if (track->fTracklets[iplane].GetChi2()>kmaxchi2) continue;
2365 if (track->fTracklets[iplane].GetN()<kmincl) continue;
69b55c55 2366 if (!(c->IsUsed())) c->Use();
5443e65e 2367 }
5443e65e 2368
75bd7f81 2369}
5443e65e 2370
75bd7f81 2371//_____________________________________________________________________________
029cd327 2372Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
5443e65e 2373{
75bd7f81 2374 //
5443e65e 2375 // Parametrised "expected" error of the cluster reconstruction in Y
75bd7f81 2376 //
5443e65e 2377
2378 Double_t s = 0.08 * 0.08;
2379 return s;
75bd7f81 2380
5443e65e 2381}
2382
75bd7f81 2383//_____________________________________________________________________________
029cd327 2384Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
0a29d0f1 2385{
75bd7f81 2386 //
5443e65e 2387 // Parametrised "expected" error of the cluster reconstruction in Z
75bd7f81 2388 //
5443e65e 2389
6c94f330 2390 Double_t s = 9 * 9 /12.;
5443e65e 2391 return s;
75bd7f81 2392
5443e65e 2393}
2394
75bd7f81 2395//_____________________________________________________________________________
029cd327 2396Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
5443e65e 2397{
2398 //
029cd327 2399 // Returns radial position which corresponds to time bin <localTB>
5443e65e 2400 // in tracking sector <sector> and plane <plane>
2401 //
2402
6c94f330 2403 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
2404 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
5443e65e 2405 return fTrSec[sector]->GetLayer(pl)->GetX();
2406
2407}
2408
75bd7f81 2409//_____________________________________________________________________________
2410AliTRDtracker::AliTRDpropagationLayer
2411 ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
2412 , Double_t radLength, Int_t tbIndex, Int_t plane)
ad670fba 2413 :fN(0)
2414 ,fSec(0)
2415 ,fClusters(NULL)
2416 ,fIndex(NULL)
2417 ,fX(x)
2418 ,fdX(dx)
2419 ,fRho(rho)
2420 ,fX0(radLength)
2421 ,fTimeBinIndex(tbIndex)
2422 ,fPlane(plane)
2423 ,fYmax(0)
2424 ,fYmaxSensitive(0)
2425 ,fHole(kFALSE)
2426 ,fHoleZc(0)
2427 ,fHoleZmax(0)
2428 ,fHoleYc(0)
2429 ,fHoleYmax(0)
2430 ,fHoleRho(0)
2431 ,fHoleX0(0)
5443e65e 2432{
0a29d0f1 2433 //
5443e65e 2434 // AliTRDpropagationLayer constructor
0a29d0f1 2435 //
46d29e70 2436
ad670fba 2437 for (Int_t i = 0; i < (Int_t) kZones; i++) {
2438 fZc[i] = 0;
2439 fZmax[i] = 0;
a819a5f7 2440 }
5443e65e 2441
ad670fba 2442 if (fTimeBinIndex >= 0) {
029cd327 2443 fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
ad670fba 2444 fIndex = new UInt_t[kMaxClusterPerTimeBin];
a819a5f7 2445 }
46d29e70 2446
ad670fba 2447 for (Int_t i = 0; i < 5; i++) {
2448 fIsHole[i] = kFALSE;
2449 }
5443e65e 2450
2451}
2452
75bd7f81 2453//_____________________________________________________________________________
2454void AliTRDtracker::AliTRDpropagationLayer
2455 ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
2456 , Double_t radLength, Double_t Yc, Double_t Zc)
5443e65e 2457{
2458 //
2459 // Sets hole in the layer
2460 //
75bd7f81 2461
6c94f330 2462 fHole = kTRUE;
2463 fHoleZc = Zc;
5443e65e 2464 fHoleZmax = Zmax;
6c94f330 2465 fHoleYc = Yc;
5443e65e 2466 fHoleYmax = Ymax;
6c94f330 2467 fHoleRho = rho;
2468 fHoleX0 = radLength;
75bd7f81 2469
5443e65e 2470}
46d29e70 2471
75bd7f81 2472//_____________________________________________________________________________
2473AliTRDtracker::AliTRDtrackingSector
ad670fba 2474 ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2475 :fN(0)
2476 ,fGeom(geo)
2477 ,fGeomSector(gs)
5443e65e 2478{
2479 //
2480 // AliTRDtrackingSector Constructor
2481 //
75bd7f81 2482
ad670fba 2483 AliTRDpadPlane *padPlane = 0;
2484 AliTRDpropagationLayer *ppl = 0;
a5cadd36 2485
ad670fba 2486 // Get holes description from geometry
3c625a9b 2487 Bool_t holes[AliTRDgeometry::kNcham];
ad670fba 2488 for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3c625a9b 2489 holes[icham] = fGeom->IsHole(0,icham,gs);
3c625a9b 2490 }
3c625a9b 2491
ad670fba 2492 for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
2493 fTimeBinIndex[i] = -1;
2494 }
5443e65e 2495
ad670fba 2496 Double_t x;
2497 Double_t dx;
2498 Double_t rho;
2499 Double_t radLength;
5443e65e 2500
ad670fba 2501 // Add layers for each of the planes
2502 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
a305677e 2503 //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
5443e65e 2504
ad670fba 2505 const Int_t kNchambers = AliTRDgeometry::Ncham();
a305677e 2506 Int_t tbIndex;
ad670fba 2507 Double_t ymax = 0;
2508 Double_t ymaxsensitive = 0;
2509 Double_t *zc = new Double_t[kNchambers];
2510 Double_t *zmax = new Double_t[kNchambers];
3c625a9b 2511 Double_t *zmaxsensitive = new Double_t[kNchambers];
5443e65e 2512
ad670fba 2513 AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
2514 if (!commonParam) {
030b4415 2515 AliErrorGeneral("AliTRDtrackingSector::Ctor"
2516 ,"Could not get common parameters\n");
3551db50 2517 return;
2518 }
2519
ad670fba 2520 for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2521
2522 ymax = fGeom->GetChamberWidth(plane) / 2.0;
2523 padPlane = commonParam->GetPadPlane(plane,0);
2524 ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
2525
2526 for (Int_t ch = 0; ch < kNchambers; ch++) {
2527 zmax[ch] = fGeom->GetChamberLength(plane,ch) / 2.0;
2528 Float_t pad = padPlane->GetRowSize(1);
2529 Float_t row0 = commonParam->GetRow0(plane,ch,0);
2530 Int_t nPads = commonParam->GetRowMax(plane,ch,0);
2531 zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
2532 zc[ch] = -(pad * nPads) / 2.0 + row0;
5443e65e 2533 }
2534
ad670fba 2535 dx = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
2536 / AliTRDcalibDB::Instance()->GetSamplingFrequency();
2537 rho = 0.00295 * 0.85; //????
2538 radLength = 11.0;
5443e65e 2539
3551db50 2540 Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
a305677e 2541 //Double_t xbottom = x0 - dxDrift;
ad670fba 2542 //Double_t xtop = x0 + dxAmp;
2543
59393e34 2544 Int_t nTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
ad670fba 2545 for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2546
2547 Double_t xlayer = iTime * dx - dxAmp;
2548 //if (xlayer<0) xlayer = dxAmp / 2.0;
59393e34 2549 x = x0 - xlayer;
ad670fba 2550
2551 tbIndex = CookTimeBinIndex(plane,iTime);
2552 ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
3c625a9b 2553 ppl->SetYmax(ymax,ymaxsensitive);
ad670fba 2554 ppl->SetZ(zc,zmax,zmaxsensitive);
3c625a9b 2555 ppl->SetHoles(holes);
59393e34 2556 InsertLayer(ppl);
ad670fba 2557
5443e65e 2558 }
ad670fba 2559
5443e65e 2560 }
2561
5443e65e 2562 MapTimeBinLayers();
ad670fba 2563
029cd327 2564 delete [] zc;
2565 delete [] zmax;
4f1c04d3 2566 delete [] zmaxsensitive;
5443e65e 2567
2568}
2569
ad670fba 2570//_____________________________________________________________________________
2571AliTRDtracker::AliTRDtrackingSector
2572 ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2573 :fN(0)
2574 ,fGeom(0)
2575 ,fGeomSector(0)
2576{
2577 //
2578 // Copy constructor
2579 //
2580
2581}
2582
75bd7f81 2583//_____________________________________________________________________________
2584Int_t AliTRDtracker::AliTRDtrackingSector
2585 ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
5443e65e 2586{
2587 //
6c94f330 2588 // depending on the digitization parameters calculates "global"
029cd327 2589 // time bin index for timebin <localTB> in plane <plane>
5443e65e 2590 //
59393e34 2591 //
75bd7f81 2592
59393e34 2593 Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
6c94f330 2594 Int_t gtb = (plane+1) * tbPerPlane - localTB -1;
2595 if (localTB<0) return -1;
2596 if (gtb<0) return -1;
75bd7f81 2597
5443e65e 2598 return gtb;
5443e65e 2599
75bd7f81 2600}
5443e65e 2601
75bd7f81 2602//_____________________________________________________________________________
2603void AliTRDtracker::AliTRDtrackingSector
2604 ::MapTimeBinLayers()
5443e65e 2605{
2606 //
2607 // For all sensitive time bins sets corresponding layer index
2608 // in the array fTimeBins
2609 //
2610
2611 Int_t index;
2612
6c94f330 2613 for(Int_t i = 0; i < fN; i++) {
5443e65e 2614 index = fLayers[i]->GetTimeBinIndex();
6c94f330 2615
2616 // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
5443e65e 2617
2618 if(index < 0) continue;
029cd327 2619 if(index >= (Int_t) kMaxTimeBinIndex) {
6c94f330 2620 printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2621 printf(" index %d exceeds allowed maximum of %d!\n",
2622 index, kMaxTimeBinIndex-1);
5443e65e 2623 continue;
2624 }
2625 fTimeBinIndex[index] = i;
2626 }
5443e65e 2627
75bd7f81 2628}
5443e65e 2629
75bd7f81 2630//_____________________________________________________________________________
2631Int_t AliTRDtracker::AliTRDtrackingSector
2632 ::GetLayerNumber(Double_t x) const
5443e65e 2633{
2634 //
2635 // Returns the number of time bin which in radial position is closest to <x>
2636 //
2637
6c94f330 2638 if(x >= fLayers[fN-1]->GetX()) return fN-1;
2639 if(x <= fLayers[0]->GetX()) return 0;
ad670fba 2640
6c94f330 2641 Int_t b=0, e=fN-1, m=(b+e)/2;
2642 for (; b<e; m=(b+e)/2) {
2643 if (x > fLayers[m]->GetX()) b=m+1;
2644 else e=m;
5443e65e 2645 }
6c94f330 2646 if(TMath::Abs(x - fLayers[m]->GetX()) >
2647 TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
75bd7f81 2648
6c94f330 2649 else return m;
5443e65e 2650
2651}
2652
75bd7f81 2653//_____________________________________________________________________________
2654Int_t AliTRDtracker::AliTRDtrackingSector
2655 ::GetInnerTimeBin() const
5443e65e 2656{
2657 //
2658 // Returns number of the innermost SENSITIVE propagation layer
2659 //
2660
2661 return GetLayerNumber(0);
5443e65e 2662
75bd7f81 2663}
5443e65e 2664
75bd7f81 2665//_____________________________________________________________________________
2666Int_t AliTRDtracker::AliTRDtrackingSector
2667 ::GetOuterTimeBin() const
5443e65e 2668{
2669 //
2670 // Returns number of the outermost SENSITIVE time bin
2671 //
2672
2673 return GetLayerNumber(GetNumberOfTimeBins() - 1);
46d29e70 2674
75bd7f81 2675}
5443e65e 2676
75bd7f81 2677//_____________________________________________________________________________
2678Int_t AliTRDtracker::AliTRDtrackingSector
2679 ::GetNumberOfTimeBins() const
5443e65e 2680{
2681 //
2682 // Returns number of SENSITIVE time bins
2683 //
2684
6c94f330 2685 Int_t tb, layer;
2686 for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
5443e65e 2687 layer = GetLayerNumber(tb);
6c94f330 2688 if(layer>=0) break;
5443e65e 2689 }
75bd7f81 2690
6c94f330 2691 return tb+1;
5443e65e 2692
75bd7f81 2693}
5443e65e 2694
75bd7f81 2695//_____________________________________________________________________________
2696void AliTRDtracker::AliTRDtrackingSector
2697 ::InsertLayer(AliTRDpropagationLayer* pl)
5443e65e 2698{
2699 //
2700 // Insert layer <pl> in fLayers array.
2701 // Layers are sorted according to X coordinate.
75bd7f81 2702 //
5443e65e 2703
6c94f330 2704 if ( fN == ((Int_t) kMaxLayersPerSector)) {
2705 printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
ad670fba 2706 return;
2707 }
6c94f330 2708 if (fN==0) {fLayers[fN++] = pl; return;}
2709 Int_t i=Find(pl->GetX());
ad670fba 2710
6c94f330 2711 memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2712 fLayers[i]=pl; fN++;
5443e65e 2713
2714}
2715
75bd7f81 2716//_____________________________________________________________________________
2717Int_t AliTRDtracker::AliTRDtrackingSector
2718 ::Find(Double_t x) const
5443e65e 2719{
2720 //
2721 // Returns index of the propagation layer nearest to X
2722 //
2723
6c94f330 2724 if (x <= fLayers[0]->GetX()) return 0;
2725 if (x > fLayers[fN-1]->GetX()) return fN;
2726 Int_t b=0, e=fN-1, m=(b+e)/2;
2727 for (; b<e; m=(b+e)/2) {
2728 if (x > fLayers[m]->GetX()) b=m+1;
2729 else e=m;
5443e65e 2730 }
7ad19338 2731
75bd7f81 2732 return m;
7ad19338 2733
75bd7f81 2734}
7ad19338 2735
75bd7f81 2736//_____________________________________________________________________________
2737void AliTRDtracker::AliTRDpropagationLayer
6c94f330 2738 ::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
3c625a9b 2739{
2740 //
6c94f330 2741 // set centers and the width of sectors
75bd7f81 2742 //
2743
6c94f330 2744 for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2745 fZc[icham] = center[icham];
2746 fZmax[icham] = w[icham];
3c625a9b 2747 fZmaxSensitive[icham] = wsensitive[icham];
6c94f330 2748 // printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
3c625a9b 2749 }
75bd7f81 2750
3c625a9b 2751}
5443e65e 2752
75bd7f81 2753//_____________________________________________________________________________
3c625a9b 2754void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
2755{
2756 //
6c94f330 2757 // set centers and the width of sectors
75bd7f81 2758 //
2759
3c625a9b 2760 fHole = kFALSE;
6c94f330 2761 for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
3c625a9b 2762 fIsHole[icham] = holes[icham];
6c94f330 2763 if (holes[icham]) fHole = kTRUE;
3c625a9b 2764 }
5443e65e 2765
75bd7f81 2766}
5443e65e 2767
75bd7f81 2768//_____________________________________________________________________________
2769void AliTRDtracker::AliTRDpropagationLayer
6c94f330 2770 ::InsertCluster(AliTRDcluster* c, UInt_t index)
75bd7f81 2771{
2772 //
2773 // Insert cluster in cluster array.
2774 // Clusters are sorted according to Y coordinate.
2775 //
5443e65e 2776
6c94f330 2777 if(fTimeBinIndex < 0) {
2778 printf("*** attempt to insert cluster into non-sensitive time bin!\n");
5443e65e 2779 return;
2780 }
ad670fba 2781
6c94f330 2782 if (fN== (Int_t) kMaxClusterPerTimeBin) {
2783 printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
ad670fba 2784 return;
2785 }
6c94f330 2786 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2787 Int_t i=Find(c->GetY());
2788 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2789 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
2790 fIndex[i]=index; fClusters[i]=c; fN++;
5443e65e 2791
75bd7f81 2792}
5443e65e 2793
75bd7f81 2794//_____________________________________________________________________________
2795Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
2796{
2797 //
2798 // Returns index of the cluster nearest in Y
2799 //
5443e65e 2800
6c94f330 2801 if (fN<=0) return 0;
2802 if (y <= fClusters[0]->GetY()) return 0;
2803 if (y > fClusters[fN-1]->GetY()) return fN;
2804 Int_t b=0, e=fN-1, m=(b+e)/2;
2805 for (; b<e; m=(b+e)/2) {
2806 if (y > fClusters[m]->GetY()) b=m+1;
2807 else e=m;
5443e65e 2808 }
75bd7f81 2809
5443e65e 2810 return m;
75bd7f81 2811
5443e65e 2812}
2813
75bd7f81 2814//_____________________________________________________________________________
2815Int_t AliTRDtracker::AliTRDpropagationLayer
2816 ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
2817 , Float_t maxroadz) const
7ad19338 2818{
2819 //
2820 // Returns index of the cluster nearest to the given y,z
2821 //
75bd7f81 2822
6c94f330 2823 Int_t index = -1;
2824 Int_t maxn = fN;
69b55c55 2825 Float_t mindist = maxroad;
6c94f330 2826 //
2827 for (Int_t i=Find(y-maxroad); i<maxn; i++) {
2828 AliTRDcluster* c=(AliTRDcluster*)(fClusters[i]);
69b55c55 2829 Float_t ycl = c->GetY();
6c94f330 2830 //
2831 if (ycl > y+maxroad) break;
2832 if (TMath::Abs(c->GetZ()-z) > maxroadz) continue;
2833 if (TMath::Abs(ycl-y)<mindist){
2834 mindist = TMath::Abs(ycl-y);
2835 index = fIndex[i];
2836 }
2837 }
75bd7f81 2838
7ad19338 2839 return index;
7ad19338 2840
75bd7f81 2841}
7ad19338 2842
75bd7f81 2843//_____________________________________________________________________________
2844Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c)
2845{
2846 //
2847 // Returns correction factor for tilted pads geometry
2848 //
5443e65e 2849
6c94f330 2850 Int_t det = c->GetDetector();
fd621f36 2851 Int_t plane = fGeom->GetPlane(det);
3551db50 2852 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
6c94f330 2853 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
b8dc2353 2854
6c94f330 2855 if(fNoTilt) h01 = 0;
75bd7f81 2856
fd621f36 2857 return h01;
5443e65e 2858
75bd7f81 2859}
c630aafd 2860
75bd7f81 2861//_____________________________________________________________________________
eab5961e 2862void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
2863{
75bd7f81 2864 //
6c94f330 2865 // *** ADDED TO GET MORE INFORMATION FOR TRD PID ---- PS
eab5961e 2866 // This is setting fdEdxPlane and fTimBinPlane
2867 // Sums up the charge in each plane for track TRDtrack and also get the
2868 // Time bin for Max. Cluster
2869 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
75bd7f81 2870 //
eab5961e 2871
6d45eaef 2872 Double_t clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
2873 Double_t maxclscharge[AliESDtrack::kNPlane];
2874 Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
2875 Int_t timebin[AliESDtrack::kNPlane];
2876
6c94f330 2877 //Initialization of cluster charge per plane.
6d45eaef 2878 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
2879 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
2880 clscharge[iPlane][iSlice] = 0.0;
6c94f330 2881 nCluster[iPlane][iSlice] = 0;
6d45eaef 2882 }
2883 }
eab5961e 2884
6c94f330 2885 //Initialization of cluster charge per plane.
f122c485 2886 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
6c94f330 2887 timebin[iPlane] = -1;
eab5961e 2888 maxclscharge[iPlane] = 0.0;
2889 }
2890
2891 // Loop through all clusters associated to track TRDtrack
2892 Int_t nClus = TRDtrack.GetNumberOfClusters(); // from Kalmantrack
2893 for (Int_t iClus = 0; iClus < nClus; iClus++) {
2894 Double_t charge = TRDtrack.GetClusterdQdl(iClus);
6c94f330 2895 Int_t index = TRDtrack.GetClusterIndex(iClus);
c6f438c0 2896 AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index);
2897 if (!pTRDcluster) continue;
2898 Int_t tb = pTRDcluster->GetLocalTimeBin();
6c94f330 2899 if (!tb) continue;
c6f438c0 2900 Int_t detector = pTRDcluster->GetDetector();
eab5961e 2901 Int_t iPlane = fGeom->GetPlane(detector);
6c94f330 2902 Int_t iSlice = tb*AliESDtrack::kNSlice/AliTRDtrack::kNtimeBins;
2903 clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice]+charge;
eab5961e 2904 if(charge > maxclscharge[iPlane]) {
2905 maxclscharge[iPlane] = charge;
6c94f330 2906 timebin[iPlane] = tb;
eab5961e 2907 }
6d45eaef 2908 nCluster[iPlane][iSlice]++;
6c94f330 2909 } // end of loop over cluster
eab5961e 2910
2911 // Setting the fdEdxPlane and fTimBinPlane variabales
c6f438c0 2912 Double_t totalCharge = 0;
6c94f330 2913
f122c485 2914 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
6d45eaef 2915 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
6c94f330 2916 if (nCluster[iPlane][iSlice]) clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
2917 TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice], iPlane, iSlice);
2918 totalCharge= totalCharge+clscharge[iPlane][iSlice];
bd50219c 2919 }
6c94f330 2920 TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
eab5961e 2921 }
6d45eaef 2922
eab5961e 2923 // Int_t i;
2924 // Int_t nc=TRDtrack.GetNumberOfClusters();
2925 // Float_t dedx=0;
2926 // for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
2927 // dedx /= nc;
2928 // for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
2929 // TRDtrack.SetPIDsignals(dedx, iPlane);
2930 // TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
2931 // }
2932
75bd7f81 2933}
c630aafd 2934
75bd7f81 2935//_____________________________________________________________________________
2936Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
6c94f330 2937 , AliTRDtrack * track
2938 , Int_t *clusters,AliTRDtracklet&tracklet)
4f1c04d3 2939{
6c94f330 2940 //
4f1c04d3 2941 //
75bd7f81 2942 // Try to find nearest clusters to the track in timebins from t0 to t1
4f1c04d3 2943 //
6c94f330 2944 //
2945 //
2946 // correction coeficients - depends on TRD parameters - to be changed according it
2947 //
ad670fba 2948
6c94f330 2949 Double_t x[100],yt[100],zt[100];
2950 Double_t xmean=0; //reference x
2951 Double_t dz[10][100],dy[10][100];
2952 Float_t zmean[100], nmean[100];
2953 Int_t clfound=0;
2954 Int_t indexes[10][100]; // indexes of the clusters in the road
2955 AliTRDcluster *cl[10][100]; // pointers to the clusters in the road
2956 Int_t best[10][100]; // index of best matching cluster
2957 //
2958 //
ad670fba 2959
6c94f330 2960 for (Int_t it=0;it<100; it++){
2961 x[it]=0;
2962 yt[it]=0;
2963 zt[it]=0;
2964 clusters[it]=-2;
2965 zmean[it]=0;
2966 nmean[it]=0;
2967 //
2968 for (Int_t ih=0;ih<10;ih++){
2969 indexes[ih][it]=-2; //reset indexes1
2970 cl[ih][it]=0;
2971 dz[ih][it]=-100;
2972 dy[ih][it]=-100;
2973 best[ih][it]=0;
2974 }
2975 }
2976 //
2977 Double_t x0 = track->GetX();
2978 Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
2979 Int_t nall=0;
2980 Int_t nfound=0;
2981 Double_t h01 =0;
2982 Int_t plane =-1;
2983 Int_t detector =-1;
2984 Float_t padlength=0;
7ad19338 2985 AliTRDtrack track2(*track);
6c94f330 2986 Float_t snpy = track->GetSnp();
2987 Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy));
2988 if (snpy<0) tany*=-1;
2989 //
2990 Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
2991 Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
2992 Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2);
2993 if (road>6.) road=6.;
4f1c04d3 2994
6c94f330 2995 //
2996 for (Int_t it=0;it<t1-t0;it++){
2997 Double_t maxChi2[2]={fgkMaxChi2,fgkMaxChi2};
2998 AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it+t0));
2999 if (timeBin==0) continue; // no indexes1
4f1c04d3 3000 Int_t maxn = timeBin;
6c94f330 3001 x[it] = timeBin.GetX();
7ad19338 3002 track2.PropagateTo(x[it]);
6c94f330 3003 yt[it] = track2.GetY();
3004 zt[it] = track2.GetZ();
7ad19338 3005
6c94f330 3006 Double_t y=yt[it],z=zt[it];
3007 Double_t chi2 =1000000;
4f1c04d3 3008 nall++;
3009 //
6c94f330 3010 // find 2 nearest cluster at given time bin
7ad19338 3011 //
6c94f330 3012 //
3013 for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
3014 AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
7ad19338 3015 h01 = GetTiltFactor(c);
6c94f330 3016 if (plane<0){
c6f438c0 3017 Int_t det = c->GetDetector();
6c94f330 3018 plane = fGeom->GetPlane(det);
3019 padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
7ad19338 3020 }
6c94f330 3021 // if (c->GetLocalTimeBin()==0) continue;
4f1c04d3 3022 if (c->GetY() > y+road) break;
6c94f330 3023 if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;
7ad19338 3024
6c94f330 3025 Double_t dist = TMath::Abs(c->GetZ()-z);
3026 if (dist> (0.5*padlength+6.*sigmaz)) continue; // 6 sigma boundary cut
3027 Double_t cost = 0;
3028 //
3029 if (dist> (0.5*padlength-sigmaz)){ // sigma boundary cost function
3030 cost = (dist-0.5*padlength)/(2.*sigmaz);
3031 if (cost>-1) cost= (cost+1.)*(cost+1.);
3032 else cost=0;
3033 }
3034 // Int_t label = TMath::Abs(track->GetLabel());
3035 // if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3036 chi2=track2.GetPredictedChi2(c,h01)+cost;
3037 //
7ad19338 3038 clfound++;
6c94f330 3039
7ad19338 3040 if (chi2 > maxChi2[1]) continue;
6c94f330 3041 detector = c->GetDetector();
3042 for (Int_t ih=2;ih<9; ih++){ //store the clusters in the road
3043 if (cl[ih][it]==0){
3044 cl[ih][it] = c;
3045 indexes[ih][it] =timeBin.GetIndex(i); // index - 9 - reserved for outliers
7ad19338 3046 break;
3047 }
4f1c04d3 3048 }
6c94f330 3049 //
3050 if (chi2 <maxChi2[0]){
7ad19338 3051 maxChi2[1] = maxChi2[0];
3052 maxChi2[0] = chi2;
3053 indexes[1][it] = indexes[0][it];
3054 cl[1][it] = cl[0][it];
3055 indexes[0][it] = timeBin.GetIndex(i);
3056 cl[0][it] = c;
3057 continue;
3058 }
6c94f330 3059 maxChi2[1]=chi2;
3060 cl[1][it] = c;
3061 indexes[1][it] =timeBin.GetIndex(i);
7ad19338 3062 }
6c94f330 3063 if (cl[0][it]){
7ad19338 3064 nfound++;
3065 xmean += x[it];
3066 }
4f1c04d3 3067 }
3068 //
6c94f330 3069 if (nfound<4) return 0;
3070 xmean /=Float_t(nfound); // middle x
3071 track2.PropagateTo(xmean); // propagate track to the center
3072 //
3073 // choose one of the variants
3074 //
3075 Int_t changes[10];
3076 Float_t sumz = 0;
3077 Float_t sum = 0;
3078 Double_t sumdy = 0;
3079 Double_t sumdy2 = 0;
3080 Double_t sumx = 0;
3081 Double_t sumxy = 0;
3082 Double_t sumx2 = 0;
3083 Double_t mpads = 0;
3084 //
3085 Int_t ngood[10];
3086 Int_t nbad[10];
7ad19338 3087 //
3088 Double_t meanz[10];
6c94f330 3089 Double_t moffset[10]; // mean offset
3090 Double_t mean[10]; // mean value
3091 Double_t angle[10]; // angle
3092 //
3093 Double_t smoffset[10]; // sigma of mean offset
3094 Double_t smean[10]; // sigma of mean value
3095 Double_t sangle[10]; // sigma of angle
3096 Double_t smeanangle[10]; // correlation
3097 //
7ad19338 3098 Double_t sigmas[10];
6c94f330 3099 Double_t tchi2s[10]; // chi2s for tracklet
ad670fba 3100
6c94f330 3101 for (Int_t it=0;it<10;it++) {
ad670fba 3102
6c94f330 3103 ngood[it] = 0;
3104 nbad[it] = 0;
3105 //
3106 meanz[it] = 0;
3107 moffset[it] = 0; // mean offset
3108 mean[it] = 0; // mean value
3109 angle[it] = 0; // angle
3110 //
3111 smoffset[it] = 1e10; // sigma of mean offset
3112 smean[it] = 1e10; // sigma of mean value
3113 sangle[it] = 1e10; // sigma of angle
3114 smeanangle[it] = 0; // correlation
3115 //
3116 sigmas[it] = 1e10;
3117 tchi2s[it] = 1e10; // chi2s for tracklet
75bd7f81 3118
3119 }
3120
7ad19338 3121 //
6c94f330 3122 // calculate zmean
7ad19338 3123 //
6c94f330 3124 for (Int_t it=0;it<t1-t0;it++){
7ad19338 3125 if (!cl[0][it]) continue;
6c94f330 3126 for (Int_t dt=-3;dt<=3;dt++){
3127 if (it+dt<0) continue;
3128 if (it+dt>t1-t0) continue;
7ad19338 3129 if (!cl[0][it+dt]) continue;
6c94f330 3130 zmean[it]+=cl[0][it+dt]->GetZ();
3131 nmean[it]+=1.;
7ad19338 3132 }
6c94f330 3133 zmean[it]/=nmean[it];
7ad19338 3134 }
6c94f330 3135 //
3136 for (Int_t it=0; it<t1-t0;it++){
3137 best[0][it]=0;
3138 for (Int_t ih=0;ih<10;ih++){
3139 dz[ih][it]=-100;
3140 dy[ih][it]=-100;
4f1c04d3 3141 if (!cl[ih][it]) continue;
6c94f330 3142 Double_t xcluster = cl[ih][it]->GetX();
3143 Double_t ytrack,ztrack;
3144 track2.GetProlongation(xcluster, ytrack, ztrack );
3145 dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // calculate distance from track in z
3146 dy[ih][it] = cl[ih][it]->GetY()+ dz[ih][it]*h01 -ytrack; // in y
7ad19338 3147 }
6c94f330 3148 // minimize changes
7ad19338 3149 if (!cl[0][it]) continue;
6c94f330 3150 if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it])
3151 if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){
3152 best[0][it]=1;
4f1c04d3 3153 }
7ad19338 3154 }
3155 //
6c94f330 3156 // iterative choosing of "best path"
7ad19338 3157 //
6c94f330 3158 //
3159 Int_t label = TMath::Abs(track->GetLabel());
3160 Int_t bestiter=0;
3161 //
3162 for (Int_t iter=0;iter<9;iter++){
3163 //
3164 changes[iter]= 0;
3165 sumz = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0;
3166 // linear fit
3167 for (Int_t it=0;it<t1-t0;it++){
7ad19338 3168 if (!cl[best[iter][it]][it]) continue;
6c94f330 3169 //calculates pad-row changes
3170 Double_t zbefore= cl[best[iter][it]][it]->GetZ();
3171 Double_t zafter = cl[best[iter][it]][it]->GetZ();
3172 for (Int_t itd = it-1; itd>=0;itd--) {
7ad19338 3173 if (cl[best[iter][itd]][itd]) {
6c94f330 3174 zbefore= cl[best[iter][itd]][itd]->GetZ();
7ad19338 3175 break;
3176 }
3177 }
6c94f330 3178 for (Int_t itd = it+1; itd<t1-t0;itd++) {
7ad19338 3179 if (cl[best[iter][itd]][itd]) {
6c94f330 3180 zafter= cl[best[iter][itd]][itd]->GetZ();
7ad19338 3181 break;
3182 }
3183 }
6c94f330 3184 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]++;
3185 //
3186 Double_t dx = x[it]-xmean; // distance to reference x
3187 sumz += cl[best[iter][it]][it]->GetZ();
4f1c04d3 3188 sum++;
6c94f330 3189 sumdy += dy[best[iter][it]][it];
3190 sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it];
3191 sumx += dx;
3192 sumx2 += dx*dx;
7ad19338 3193 sumxy += dx*dy[best[iter][it]][it];
6c94f330 3194 mpads += cl[best[iter][it]][it]->GetNPads();
3195 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){
7ad19338 3196 ngood[iter]++;
4f1c04d3 3197 }
6c94f330 3198 else{
7ad19338 3199 nbad[iter]++;
4f1c04d3 3200 }
3201 }
7ad19338 3202 //
6c94f330 3203 // calculates line parameters
7ad19338 3204 //
6c94f330 3205 Double_t det = sum*sumx2-sumx*sumx;
3206 angle[iter] = (sum*sumxy-sumx*sumdy)/det;
3207 mean[iter] = (sumx2*sumdy-sumx*sumxy)/det;
3208 meanz[iter] = sumz/sum;
3209 moffset[iter] = sumdy/sum;
3210 mpads /= sum; // mean number of pads
3211 //
3212 //
3213 Double_t sigma2 = 0; // normalized residuals - for line fit
3214 Double_t sigma1 = 0; // normalized residuals - constant fit
3215 //
3216 for (Int_t it=0;it<t1-t0;it++){
7ad19338 3217 if (!cl[best[iter][it]][it]) continue;
6c94f330 3218 Double_t dx = x[it]-xmean;
3219 Double_t ytr = mean[iter]+angle[iter]*dx;
3220 sigma2 += (dy[best[iter][it]][it]-ytr)*(dy[best[iter][it]][it]-ytr);
3221 sigma1 += (dy[best[iter][it]][it]-moffset[iter])*(dy[best[iter][it]][it]-moffset[iter]);
7ad19338 3222 sum++;
4f1c04d3 3223 }
6c94f330 3224 sigma2 /=(sum-2); // normalized residuals
3225 sigma1 /=(sum-1); // normalized residuals
7ad19338 3226 //
6c94f330 3227 smean[iter] = sigma2*(sumx2/det); // estimated error2 of mean
3228 sangle[iter] = sigma2*(sum/det); // estimated error2 of angle
3229 smeanangle[iter] = sigma2*(-sumx/det); // correlation
7ad19338 3230 //
6c94f330 3231 //
3232 sigmas[iter] = TMath::Sqrt(sigma1); //
3233 smoffset[iter]= (sigma1/sum)+0.01*0.01; // sigma of mean offset + unisochronity sigma
3234 //
3235 // iterative choosing of "better path"
3236 //
3237 for (Int_t it=0;it<t1-t0;it++){
7ad19338 3238 if (!cl[best[iter][it]][it]) continue;
6c94f330 3239 //
3240 Double_t sigmatr2 = smoffset[iter]+0.5*tany*tany; //add unisochronity + angular effect contribution
3241 Double_t sweight = 1./sigmatr2+1./track->GetSigmaY2();
3242 Double_t weighty = (moffset[iter]/sigmatr2)/sweight; // weighted mean
3243 Double_t sigmacl = TMath::Sqrt(sigma1*sigma1+track->GetSigmaY2()); //
3244 Double_t mindist=100000;
3245 Int_t ihbest=0;
3246 for (Int_t ih=0;ih<10;ih++){
7ad19338 3247 if (!cl[ih][it]) break;
6c94f330 3248 Double_t dist2 = (dy[ih][it]-weighty)/sigmacl;
3249 dist2*=dist2; //chi2 distance
3250 if (dist2<mindist){
7ad19338 3251 mindist = dist2;
6c94f330 3252 ihbest =ih;
7ad19338 3253 }
3254 }
6c94f330 3255 best[iter+1][it]=ihbest;
4f1c04d3 3256 }
4f1c04d3 3257 //
6c94f330 3258 // update best hypothesy if better chi2 according tracklet position and angle
7ad19338 3259 //
69b55c55 3260 Double_t sy2 = smean[iter] + track->GetSigmaY2();
6c94f330 3261 Double_t sa2 = sangle[iter] + track->GetSigmaSnp2();//track->fCee;
3262 Double_t say = track->GetSigmaSnpY();//track->fCey;
3263 // Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
3264 // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
3265
3266 Double_t detchi = sy2*sa2-say*say;
3267 Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi}; //inverse value of covariance matrix
4f1c04d3 3268
6c94f330 3269 Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+
3270 2.*mean[bestiter]*angle[bestiter]*invers[2];
3271 Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+
3272 2*mean[iter]*angle[iter]*invers[2];
3273 tchi2s[iter] =chi21;
3274 //
3275 if (changes[iter]<=changes[bestiter] && chi21<chi20) {
3276 bestiter =iter;
7ad19338 3277 }
3278 }
3279 //
6c94f330 3280 //set clusters
7ad19338 3281 //
6c94f330 3282 Double_t sigma2 = sigmas[0]; // choose as sigma from 0 iteration
3283 Short_t maxpos = -1;
3284 Float_t maxcharge = 0;
3285 Short_t maxpos4 = -1;
3286 Float_t maxcharge4 = 0;
3287 Short_t maxpos5 = -1;
3288 Float_t maxcharge5 = 0;
8979685e 3289
7ad19338 3290 //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3291 //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
3292
6c94f330 3293 Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0));
3294 Double_t expectederr = sigma2*sigma2+0.01*0.01;
3295 if (mpads>3.5) expectederr += (mpads-3.5)*0.04;
3296 if (changes[bestiter]>1) expectederr+= changes[bestiter]*0.01;
3297 expectederr+=(0.03*(tany-exB)*(tany-exB))*15;
3298 // if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
7ad19338 3299 //expectederr+=10000;
6c94f330 3300 for (Int_t it=0;it<t1-t0;it++){
7ad19338 3301 if (!cl[best[bestiter][it]][it]) continue;
6c94f330 3302 cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // set cluster error
3303 if (!cl[best[bestiter][it]][it]->IsUsed()){
3304 cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY());
3305 // cl[best[bestiter][it]][it]->Use();
69b55c55 3306 }
6c94f330 3307 //
3308 // time bins with maximal charge
3309 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
69b55c55 3310 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
6c94f330 3311 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3312 }
6c94f330 3313
3314 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
3315 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
69b55c55 3316 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
6c94f330 3317 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3318 }
3319 }
6c94f330 3320 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
3321 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
69b55c55 3322 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
6c94f330 3323 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
69b55c55 3324 }
7ad19338 3325 }
6c94f330 3326 //
3327 // time bins with maximal charge
3328 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
8979685e 3329 maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
6c94f330 3330 maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3331 }
6c94f330 3332
3333 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
3334 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
8979685e 3335 maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
6c94f330 3336 maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3337 }
3338 }
6c94f330 3339 if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
3340 if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
8979685e 3341 maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
6c94f330 3342 maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
8979685e 3343 }
3344 }
7ad19338 3345 clusters[it+t0] = indexes[best[bestiter][it]][it];
6c94f330 3346 //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 3347 }
3348 //
6c94f330 3349 // set tracklet parameters
7ad19338 3350 //
6c94f330 3351 Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01;
3352 if (mpads>3.5) trackleterr2 += (mpads-3.5)*0.04;
3353 trackleterr2+= changes[bestiter]*0.01;
3354 trackleterr2*= TMath::Max(14.-nfound,1.);
3355 trackleterr2+= 0.2*(tany-exB)*(tany-exB);
3356 //
3357 tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2); //set tracklet parameters
7ad19338 3358 tracklet.SetTilt(h01);
3359 tracklet.SetP0(mean[bestiter]);
3360 tracklet.SetP1(angle[bestiter]);
3361 tracklet.SetN(nfound);
3362 tracklet.SetNCross(changes[bestiter]);
3363 tracklet.SetPlane(plane);
3364 tracklet.SetSigma2(expectederr);
3365 tracklet.SetChi2(tchi2s[bestiter]);
8979685e 3366 tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
6c94f330 3367 track->fTracklets[plane] = tracklet;
3368 track->fNWrong+=nbad[0];
7ad19338 3369 //
3370 // Debuging part
3371 //
69b55c55 3372 TClonesArray array0("AliTRDcluster");
3373 TClonesArray array1("AliTRDcluster");
6c94f330 3374 array0.ExpandCreateFast(t1-t0+1);
3375 array1.ExpandCreateFast(t1-t0+1);
3376 TTreeSRedirector& cstream = *fDebugStreamer;
7ad19338 3377 AliTRDcluster dummy;
3378 Double_t dy0[100];
8979685e 3379 Double_t dyb[100];
3380
6c94f330 3381 for (Int_t it=0;it<t1-t0;it++){
7ad19338 3382 dy0[it] = dy[0][it];
3383 dyb[it] = dy[best[bestiter][it]][it];
3384 if(cl[0][it]) {
3385 new(array0[it]) AliTRDcluster(*cl[0][it]);
3386 }
3387 else{
3388 new(array0[it]) AliTRDcluster(dummy);
3389 }
3390 if(cl[best[bestiter][it]][it]) {
3391 new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3392 }
3393 else{
3394 new(array1[it]) AliTRDcluster(dummy);
3395 }
4f1c04d3 3396 }
7ad19338 3397 TGraph graph0(t1-t0,x,dy0);
3398 TGraph graph1(t1-t0,x,dyb);
3399 TGraph graphy(t1-t0,x,yt);
3400 TGraph graphz(t1-t0,x,zt);
6c94f330 3401 //
3402 //
3403 if (AliTRDReconstructor::StreamLevel()>0)
3404 cstream<<"tracklet"<<
3405 "track.="<<track<< // track parameters
3406 "tany="<<tany<< // tangent of the local track angle
3407 "xmean="<<xmean<< // xmean - reference x of tracklet
3408 "tilt="<<h01<< // tilt angle
3409 "nall="<<nall<< // number of foundable clusters
3410 "nfound="<<nfound<< // number of found clusters
3411 "clfound="<<clfound<< // total number of found clusters in road
3412 "mpads="<<mpads<< // mean number of pads per cluster
3413 "plane="<<plane<< // plane number
3414 "detector="<<detector<< // detector number
3415 "road="<<road<< // the width of the used road
3416 "graph0.="<<&graph0<< // x - y = dy for closest cluster
3417 "graph1.="<<&graph1<< // x - y = dy for second closest cluster
3418 "graphy.="<<&graphy<< // y position of the track
3419 "graphz.="<<&graphz<< // z position of the track
3420 // "fCl.="<<&array0<< // closest cluster
3421 //"fCl2.="<<&array1<< // second closest cluster
3422 "maxpos="<<maxpos<< // maximal charge postion
3423 "maxcharge="<<maxcharge<< // maximal charge
3424 "maxpos4="<<maxpos4<< // maximal charge postion - after bin 4
3425 "maxcharge4="<<maxcharge4<< // maximal charge - after bin 4
3426 "maxpos5="<<maxpos5<< // maximal charge postion - after bin 5
3427 "maxcharge5="<<maxcharge5<< // maximal charge - after bin 5
3428 //
3429 "bestiter="<<bestiter<< // best iteration number
3430 "tracklet.="<<&tracklet<< // corrspond to the best iteration
3431 "tchi20="<<tchi2s[0]<< // chi2 of cluster in the 0 iteration
3432 "tchi2b="<<tchi2s[bestiter]<< // chi2 of cluster in the best iteration
3433 "sigmas0="<<sigmas[0]<< // residuals sigma
3434 "sigmasb="<<sigmas[bestiter]<< // residulas sigma
3435 //
3436 "ngood0="<<ngood[0]<< // number of good clusters in 0 iteration
3437 "nbad0="<<nbad[0]<< // number of bad clusters in 0 iteration
3438 "ngoodb="<<ngood[bestiter]<< // in best iteration
3439 "nbadb="<<nbad[bestiter]<< // in best iteration
3440 //
3441 "changes0="<<changes[0]<< // changes of pardrows in iteration number 0
3442 "changesb="<<changes[bestiter]<< // changes of pardrows in best iteration
3443 //
3444 "moffset0="<<moffset[0]<< // offset fixing angle in iter=0
3445 "smoffset0="<<smoffset[0]<< // sigma of offset fixing angle in iter=0
3446 "moffsetb="<<moffset[bestiter]<< // offset fixing angle in iter=best
3447 "smoffsetb="<<smoffset[bestiter]<< // sigma of offset fixing angle in iter=best
3448 //
3449 "mean0="<<mean[0]<< // mean dy in iter=0;
3450 "smean0="<<smean[0]<< // sigma of mean dy in iter=0
3451 "meanb="<<mean[bestiter]<< // mean dy in iter=best
3452 "smeanb="<<smean[bestiter]<< // sigma of mean dy in iter=best
3453 //
3454 "angle0="<<angle[0]<< // angle deviation in the iteration number 0
3455 "sangle0="<<sangle[0]<< // sigma of angular deviation in iteration number 0
3456 "angleb="<<angle[bestiter]<< // angle deviation in the best iteration
3457 "sangleb="<<sangle[bestiter]<< // sigma of angle deviation in the best iteration
3458 //
3459 "expectederr="<<expectederr<< // expected error of cluster position
3460 "\n";
3461 //
3462 //
4f1c04d3 3463 return nfound;
4f1c04d3 3464
75bd7f81 3465}
4f1c04d3 3466
75bd7f81 3467//_____________________________________________________________________________
3468Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
3469 , Int_t *outlist, Bool_t down)
69b55c55 3470{
3471 //
3472 // Sort eleements according occurancy
3473 // The size of output array has is 2*n
3474 //
75bd7f81 3475
6c94f330 3476 Int_t * sindexS = new Int_t[n]; // temp array for sorting
3477 Int_t * sindexF = new Int_t[2*n];
3478 for (Int_t i=0;i<n;i++) sindexF[i]=0;
3479 //
3480 TMath::Sort(n,inlist, sindexS, down);
3481 Int_t last = inlist[sindexS[0]];
3482 Int_t val = last;
3483 sindexF[0] = 1;
3484 sindexF[0+n] = last;
3485 Int_t countPos = 0;
3486 //
3487 // find frequency
3488 for(Int_t i=1;i<n; i++){
69b55c55 3489 val = inlist[sindexS[i]];
6c94f330 3490 if (last == val) sindexF[countPos]++;
3491 else{
69b55c55 3492 countPos++;
3493 sindexF[countPos+n] = val;
3494 sindexF[countPos]++;
6c94f330 3495 last =val;
69b55c55 3496 }
3497 }
6c94f330 3498 if (last==val) countPos++;
3499 // sort according frequency
3500 TMath::Sort(countPos, sindexF, sindexS, kTRUE);
3501 for (Int_t i=0;i<countPos;i++){
69b55c55 3502 outlist[2*i ] = sindexF[sindexS[i]+n];
3503 outlist[2*i+1] = sindexF[sindexS[i]];
3504 }
3505 delete [] sindexS;
3506 delete [] sindexF;
3507
3508 return countPos;
75bd7f81 3509
69b55c55 3510}
3511
75bd7f81 3512//_____________________________________________________________________________
6c94f330 3513AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
69b55c55 3514{
3515 //
75bd7f81 3516 // Register a seed
69b55c55 3517 //
6c94f330 3518 Double_t alpha=AliTRDgeometry::GetAlpha();
3519 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
69b55c55 3520 Double_t c[15];
6c94f330 3521 c[0] = 0.2;
3522 c[1] = 0 ; c[2] = 2;
3523 c[3] = 0 ; c[4] = 0; c[5] = 0.02;
3524 c[6] = 0 ; c[7] = 0; c[8] = 0; c[9] = 0.1;
3525 c[10] = 0 ; c[11] = 0; c[12] = 0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
3526 //
3527 Int_t index =0;
3528 AliTRDcluster *cl =0;
3529 for (Int_t ilayer=0;ilayer<6;ilayer++){
3530 if (seeds[ilayer].IsOK()){
3531 for (Int_t itime=22;itime>0;itime--){
3532 if (seeds[ilayer].fIndexes[itime]>0){
69b55c55 3533 index = seeds[ilayer].fIndexes[itime];
6c94f330 3534 cl = seeds[ilayer].fClusters[itime];
69b55c55 3535 break;
3536 }
3537 }
3538 }
6c94f330 3539 if (index>0) break;
69b55c55 3540 }
6c94f330 3541 if (cl==0) return 0;
3542 AliTRDtrack * track = new AliTRDtrack(cl,index,&params[1],c, params[0],params[6]*alpha+shift);
3543 track->PropagateTo(params[0]-5.);
69b55c55 3544 track->ResetCovariance(1);
6c94f330 3545 //
3546 Int_t rc=FollowBackProlongation(*track);
3547 if (rc<30) {
69b55c55 3548 delete track;
6c94f330 3549 track =0;
3550 }else{
69b55c55 3551 track->CookdEdx();
3552 CookdEdxTimBin(*track);
6c94f330 3553 CookLabel(track, 0.9);
69b55c55 3554 }
69b55c55 3555
75bd7f81 3556 return track;
69b55c55 3557
69b55c55 3558}