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