1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.24 2003/02/19 09:02:28 hristov
19 Track time measurement (S.Radomski)
21 Revision 1.23 2003/02/10 14:06:10 cblume
22 Add tracking without tilted pads as option
24 Revision 1.22 2003/01/30 15:19:58 cblume
27 Revision 1.21 2003/01/27 16:34:49 cblume
28 Update of tracking by Sergei and Chuncheng
30 Revision 1.20 2002/11/07 15:52:09 cblume
31 Update of tracking code for tilted pads
33 Revision 1.19 2002/10/22 15:53:08 alibrary
34 Introducing Riostream.h
36 Revision 1.18 2002/10/14 14:57:44 hristov
37 Merging the VirtualMC branch to the main development branch (HEAD)
39 Revision 1.14.6.2 2002/07/24 10:09:31 alibrary
42 Revision 1.17 2002/06/13 12:09:58 hristov
45 Revision 1.16 2002/06/12 09:54:36 cblume
46 Update of tracking code provided by Sergei
48 Revision 1.14 2001/11/14 10:50:46 cblume
49 Changes in digits IO. Add merging of summable digits
51 Revision 1.13 2001/05/30 12:17:47 hristov
52 Loop variables declared once
54 Revision 1.12 2001/05/28 17:07:58 hristov
55 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
57 Revision 1.8 2000/12/20 13:00:44 cblume
58 Modifications for the HP-compiler
60 Revision 1.7 2000/12/08 16:07:02 cblume
61 Update of the tracking by Sergei
63 Revision 1.6 2000/11/30 17:38:08 cblume
64 Changes to get in line with new STEER and EVGEN
66 Revision 1.5 2000/11/14 14:40:27 cblume
67 Correction for the Sun compiler (kTRUE and kFALSE)
69 Revision 1.4 2000/11/10 14:57:52 cblume
70 Changes in the geometry constants for the DEC compiler
72 Revision 1.3 2000/10/15 23:40:01 cblume
75 Revision 1.2 2000/10/06 16:49:46 cblume
78 Revision 1.1.2.2 2000/10/04 16:34:58 cblume
79 Replace include files by forward declarations
81 Revision 1.1.2.1 2000/09/22 14:47:52 cblume
86 #include <Riostream.h>
91 #include <TObjArray.h>
93 #include "AliTRDgeometry.h"
94 #include "AliTRDparameter.h"
95 #include "AliTRDgeometryDetail.h"
96 #include "AliTRDcluster.h"
97 #include "AliTRDtrack.h"
98 #include "../TPC/AliTPCtrack.h"
100 #include "AliTRDtracker.h"
102 ClassImp(AliTRDtracker)
104 const Float_t AliTRDtracker::fSeedDepth = 0.5;
105 const Float_t AliTRDtracker::fSeedStep = 0.10;
106 const Float_t AliTRDtracker::fSeedGap = 0.25;
108 const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.;
109 const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.;
110 const Float_t AliTRDtracker::fMaxSeedC = 0.0052;
111 const Float_t AliTRDtracker::fMaxSeedTan = 1.2;
112 const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.;
114 const Double_t AliTRDtracker::fSeedErrorSY = 0.2;
115 const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5;
116 const Double_t AliTRDtracker::fSeedErrorSZ = 0.1;
118 const Float_t AliTRDtracker::fMinClustersInSeed = 0.7;
120 const Float_t AliTRDtracker::fMinClustersInTrack = 0.5;
121 const Float_t AliTRDtracker::fMinFractionOfFoundClusters = 0.8;
123 const Float_t AliTRDtracker::fSkipDepth = 0.3;
124 const Float_t AliTRDtracker::fLabelFraction = 0.8;
125 const Float_t AliTRDtracker::fWideRoad = 20.;
127 const Double_t AliTRDtracker::fMaxChi2 = 12.;
130 //____________________________________________________________________
131 AliTRDtracker::AliTRDtracker(const TFile *geomfile)
139 fAddTRDseeds = kFALSE;
143 TDirectory *savedir=gDirectory;
144 TFile *in=(TFile*)geomfile;
146 printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
147 printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
152 fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
153 fPar = (AliTRDparameter*) in->Get("TRDparameter");
158 // fTzero = geo->GetT0();
159 printf("Found geometry version %d on file \n", fGeom->IsVersion());
162 printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n");
163 printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
164 fGeom = new AliTRDgeometryDetail();
165 fPar = new AliTRDparameter();
171 // fGeom->SetT0(fTzero);
174 fClusters = new TObjArray(2000);
176 fSeeds = new TObjArray(2000);
178 fTracks = new TObjArray(1000);
180 for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
181 Int_t tr_s = CookSectorIndex(geom_s);
182 fTrSec[tr_s] = new AliTRDtrackingSector(fGeom, geom_s, fPar);
185 Float_t tilt_angle = TMath::Abs(fPar->GetTiltingAngle());
186 if(tilt_angle < 0.1) {
193 if(fNoTilt && (tilt_angle > 0.1)) fSY2corr = fSY2corr + tilt_angle * 0.05;
196 // calculate max gap on track
198 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
199 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
201 Double_t dx = (Double_t) fPar->GetTimeBinSize();
202 Int_t tbAmp = fPar->GetTimeBefore();
203 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
204 if(kTRUE) maxAmp = 0; // intentional until we change the parameter class
205 Int_t tbDrift = fPar->GetTimeMax();
206 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
208 tbDrift = TMath::Min(tbDrift,maxDrift);
209 tbAmp = TMath::Min(tbAmp,maxAmp);
211 fTimeBinsPerPlane = tbAmp + tbDrift;
212 fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth);
218 //___________________________________________________________________
219 AliTRDtracker::~AliTRDtracker()
227 for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
228 delete fTrSec[geom_s];
232 //_____________________________________________________________________
233 inline Double_t f1trd(Double_t x1,Double_t y1,
234 Double_t x2,Double_t y2,
235 Double_t x3,Double_t y3)
238 // Initial approximation of the track curvature
240 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
241 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
242 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
243 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
244 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
246 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
248 return -xr*yr/sqrt(xr*xr+yr*yr);
251 //_____________________________________________________________________
252 inline Double_t f2trd(Double_t x1,Double_t y1,
253 Double_t x2,Double_t y2,
254 Double_t x3,Double_t y3)
257 // Initial approximation of the track curvature times X coordinate
258 // of the center of curvature
261 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
262 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
263 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
264 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
265 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
267 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
269 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
272 //_____________________________________________________________________
273 inline Double_t f3trd(Double_t x1,Double_t y1,
274 Double_t x2,Double_t y2,
275 Double_t z1,Double_t z2)
278 // Initial approximation of the tangent of the track dip angle
281 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
284 //___________________________________________________________________
285 Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
288 // Finds tracks within the TRD. File <inp> is expected to contain seeds
289 // at the outer part of the TRD. If <inp> is NULL, the seeds
290 // are found within the TRD if fAddTRDseeds is TRUE.
291 // The tracks are propagated to the innermost time bin
292 // of the TRD and stored in file <out>.
297 TDirectory *savedir=gDirectory;
301 if (!out->IsOpen()) {
302 cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
306 sprintf(tname,"seedTRDtoTPC_%d",GetEventNumber());
307 TTree tpc_tree(tname,"Tree with seeds from TRD at outer TPC pad row");
308 AliTPCtrack *iotrack=0;
309 tpc_tree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
311 sprintf(tname,"TreeT%d_TRD",GetEventNumber());
312 TTree trd_tree(tname,"TRD tracks at inner TRD time bin");
313 AliTRDtrack *iotrack_trd=0;
314 trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0);
316 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
317 Float_t foundMin = fMinClustersInTrack * timeBins;
320 TFile *in=(TFile*)inp;
322 cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
323 cerr<<" ... going for seeds finding inside the TRD\n";
327 sprintf(tname,"TRDb_%d",GetEventNumber());
328 TTree *seedTree=(TTree*)in->Get(tname);
330 cerr<<"AliTRDtracker::Clusters2Tracks(): ";
331 cerr<<"can't get a tree with track seeds !\n";
334 AliTRDtrack *seed=new AliTRDtrack;
335 seedTree->SetBranchAddress("tracks",&seed);
337 Int_t n=(Int_t)seedTree->GetEntries();
338 for (Int_t i=0; i<n; i++) {
339 seedTree->GetEvent(i);
340 seed->ResetCovariance();
341 AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
353 // find tracks from loaded seeds
355 Int_t nseed=fSeeds->GetEntriesFast();
357 Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
359 for (i=0; i<nseed; i++) {
360 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
361 FollowProlongation(t, innerTB);
362 if (t.GetNumberOfClusters() >= foundMin) {
364 CookLabel(pt, 1-fLabelFraction);
370 // cout<<found<<'\r';
372 if(PropagateToTPC(t)) {
373 AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
378 delete fSeeds->RemoveAt(i);
382 cout<<"Number of loaded seeds: "<<nseed<<endl;
383 cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
385 // after tracks from loaded seeds are found and the corresponding
386 // clusters are used, look for additional seeds from TRD
389 // Find tracks for the seeds in the TRD
390 Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
392 Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
393 Int_t gap = (Int_t) (timeBins * fSeedGap);
394 Int_t step = (Int_t) (timeBins * fSeedStep);
396 // make a first turn with tight cut on initial curvature
397 for(Int_t turn = 1; turn <= 2; turn++) {
399 nSteps = (Int_t) (fSeedDepth / (3*fSeedStep));
400 step = (Int_t) (timeBins * (3*fSeedStep));
402 for(Int_t i=0; i<nSteps; i++) {
403 Int_t outer=timeBins-1-i*step;
404 Int_t inner=outer-gap;
406 nseed=fSeeds->GetEntriesFast();
408 MakeSeeds(inner, outer, turn);
410 nseed=fSeeds->GetEntriesFast();
411 printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
414 for (Int_t i=0; i<nseed; i++) {
415 AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
416 FollowProlongation(t,innerTB);
417 if (t.GetNumberOfClusters() >= foundMin) {
419 CookLabel(pt, 1-fLabelFraction);
422 // cout<<found<<'\r';
425 if(PropagateToTPC(t)) {
426 AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
432 delete fSeeds->RemoveAt(i);
441 cout<<"Total number of found tracks: "<<found<<endl;
452 //_____________________________________________________________________________
453 Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
455 // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
456 // backpropagated by the TPC tracker. Each seed is first propagated
457 // to the TRD, and then its prolongation is searched in the TRD.
458 // If sufficiently long continuation of the track is found in the TRD
459 // the track is updated, otherwise it's stored as originaly defined
460 // by the TPC tracker.
465 TDirectory *savedir=gDirectory;
467 TFile *in=(TFile*)inp;
470 cerr<<"AliTRDtracker::PropagateBack(): ";
471 cerr<<"file with back propagated TPC tracks is not open !\n";
475 if (!out->IsOpen()) {
476 cerr<<"AliTRDtracker::PropagateBack(): ";
477 cerr<<"file for back propagated TRD tracks is not open !\n";
483 sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());
484 TTree *seedTree=(TTree*)in->Get(tname);
486 cerr<<"AliTRDtracker::PropagateBack(): ";
487 cerr<<"can't get a tree with seeds from TPC !\n";
488 cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
492 AliTPCtrack *seed=new AliTPCtrack;
493 seedTree->SetBranchAddress("tracks",&seed);
495 Int_t n=(Int_t)seedTree->GetEntries();
496 for (Int_t i=0; i<n; i++) {
497 seedTree->GetEvent(i);
498 Int_t lbl = seed->GetLabel();
499 AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
500 tr->SetSeedLabel(lbl);
510 AliTPCtrack *otrack=0;
512 sprintf(tname,"seedsTRDtoTOF1_%d",GetEventNumber());
513 TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
514 tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);
516 sprintf(tname,"seedsTRDtoTOF2_%d",GetEventNumber());
517 TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
518 tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);
520 sprintf(tname,"seedsTRDtoPHOS_%d",GetEventNumber());
521 TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
522 phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
524 sprintf(tname,"seedsTRDtoRICH_%d",GetEventNumber());
525 TTree richTree(tname,"Tracks back propagated through TPC and TRD");
526 richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
528 sprintf(tname,"TRDb_%d",GetEventNumber());
529 TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
530 AliTRDtrack *otrack_trd=0;
531 trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0);
535 Int_t nseed=fSeeds->GetEntriesFast();
537 // Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan();
538 Float_t foundMin = 0;
540 Int_t outermost_tb = fTrSec[0]->GetOuterTimeBin();
542 for (Int_t i=0; i<nseed; i++) {
544 AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
545 Int_t expectedClr = FollowBackProlongation(s);
546 Int_t foundClr = s.GetNumberOfClusters();
547 Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX());
549 // printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
550 // i, foundClr, expectedClr, foundMin);
552 if (foundClr >= foundMin) {
555 CookLabel(ps, 1-fLabelFraction);
559 // Propagate to outer reference plane [SR, GSI, 18.02.2003]
560 ps->PropagateTo(364.8);
564 // cout<<found<<'\r';
567 if(((expectedClr < 10) && (last_tb == outermost_tb)) ||
568 ((expectedClr >= 10) &&
569 (((Float_t) foundClr) / ((Float_t) expectedClr) >=
570 fMinFractionOfFoundClusters) && (last_tb == outermost_tb))) {
572 Double_t x_tof = 375.5;
574 if(PropagateToOuterPlane(s,x_tof)) {
575 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
582 if(PropagateToOuterPlane(s,x_tof)) {
583 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
588 Double_t x_phos = 460.;
590 if(PropagateToOuterPlane(s,x_phos)) {
591 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
596 Double_t x_rich = 490+1.267;
598 if(PropagateToOuterPlane(s,x_rich)) {
599 AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
617 cerr<<"Number of seeds: "<<nseed<<endl;
618 cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
627 //---------------------------------------------------------------------------
628 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
630 // Starting from current position on track=t this function tries
631 // to extrapolate the track up to timeBin=0 and to confirm prolongation
632 // if a close cluster is found. Returns the number of clusters
633 // expected to be found in sensitive layers
635 Float_t wIndex, wTB, wChi2;
636 Float_t wYrt, wYclosest, wYcorrect, wYwindow;
637 Float_t wZrt, wZclosest, wZcorrect, wZwindow;
638 Float_t wPx, wPy, wPz, wC;
640 Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
642 Int_t trackIndex = t.GetLabel();
644 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
646 Int_t try_again=fMaxGap;
648 Double_t alpha=t.GetAlpha();
650 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
651 if (alpha < 0. ) alpha += 2.*TMath::Pi();
653 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
655 Double_t rad_length, rho, x, dx, y, ymax, z;
657 Int_t expectedNumberOfClusters = 0;
658 Bool_t lookForCluster;
660 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
663 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
665 y = t.GetY(); z = t.GetZ();
667 // first propagate to the inner surface of the current time bin
668 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
669 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
670 if(!t.PropagateTo(x,rad_length,rho)) break;
672 ymax = x*TMath::Tan(0.5*alpha);
675 if (!t.Rotate(alpha)) break;
676 if(!t.PropagateTo(x,rad_length,rho)) break;
677 } else if (y <-ymax) {
679 if (!t.Rotate(-alpha)) break;
680 if(!t.PropagateTo(x,rad_length,rho)) break;
683 y = t.GetY(); z = t.GetZ();
685 // now propagate to the middle plane of the next time bin
686 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
687 x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
688 if(!t.PropagateTo(x,rad_length,rho)) break;
690 ymax = x*TMath::Tan(0.5*alpha);
693 if (!t.Rotate(alpha)) break;
694 if(!t.PropagateTo(x,rad_length,rho)) break;
695 } else if (y <-ymax) {
697 if (!t.Rotate(-alpha)) break;
698 if(!t.PropagateTo(x,rad_length,rho)) break;
704 expectedNumberOfClusters++;
705 wIndex = (Float_t) t.GetLabel();
708 AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr-1));
710 Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
711 Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
714 if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
715 else return expectedNumberOfClusters;
719 wYwindow = (Float_t) road;
720 t.GetPxPyPz(Px,Py,Pz);
724 wC = (Float_t) t.GetC();
725 wSigmaC2 = (Float_t) t.GetSigmaC2();
726 wSigmaTgl2 = (Float_t) t.GetSigmaTgl2();
727 wSigmaY2 = (Float_t) t.GetSigmaY2();
728 wSigmaZ2 = (Float_t) t.GetSigmaZ2();
731 if (road>fWideRoad) {
732 if (t.GetNumberOfClusters()>4)
733 cerr<<t.GetNumberOfClusters()
734 <<"FindProlongation warning: Too broad road !\n";
741 Double_t max_chi2=fMaxChi2;
743 wYclosest = 12345678;
744 wYcorrect = 12345678;
745 wZclosest = 12345678;
746 wZcorrect = 12345678;
747 wZwindow = TMath::Sqrt(2.25 * 12 * sz2);
749 // Find the closest correct cluster for debugging purposes
751 Float_t minDY = 1000000;
752 for (Int_t i=0; i<time_bin; i++) {
753 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
754 if((c->GetLabel(0) != trackIndex) &&
755 (c->GetLabel(1) != trackIndex) &&
756 (c->GetLabel(2) != trackIndex)) continue;
757 if(TMath::Abs(c->GetY() - y) > minDY) continue;
758 minDY = TMath::Abs(c->GetY() - y);
759 wYcorrect = c->GetY();
760 wZcorrect = c->GetZ();
762 Double_t h01 = GetTiltFactor(c);
763 wChi2 = t.GetPredictedChi2(c, h01);
767 // Now go for the real cluster search
771 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
772 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
773 if (c->GetY() > y+road) break;
774 if (c->IsUsed() > 0) continue;
775 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
777 Double_t h01 = GetTiltFactor(c);
778 Double_t chi2=t.GetPredictedChi2(c,h01);
780 if (chi2 > max_chi2) continue;
783 index=time_bin.GetIndex(i);
788 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
789 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
791 if (c->GetY() > y+road) break;
792 if (c->IsUsed() > 0) continue;
793 if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
795 Double_t h01 = GetTiltFactor(c);
796 Double_t chi2=t.GetPredictedChi2(c, h01);
798 if (chi2 > max_chi2) continue;
801 index=time_bin.GetIndex(i);
807 wYclosest = cl->GetY();
808 wZclosest = cl->GetZ();
809 Double_t h01 = GetTiltFactor(cl);
811 t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
812 if(!t.Update(cl,max_chi2,index,h01)) {
813 if(!try_again--) return 0;
815 else try_again=fMaxGap;
818 if (try_again==0) break;
823 if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
825 printf(" %f", wIndex); //1
826 printf(" %f", wTB); //2
827 printf(" %f", wYrt); //3
828 printf(" %f", wYclosest); //4
829 printf(" %f", wYcorrect); //5
830 printf(" %f", wYwindow); //6
831 printf(" %f", wZrt); //7
832 printf(" %f", wZclosest); //8
833 printf(" %f", wZcorrect); //9
834 printf(" %f", wZwindow); //10
835 printf(" %f", wPx); //11
836 printf(" %f", wPy); //12
837 printf(" %f", wPz); //13
838 printf(" %f", wSigmaC2*1000000); //14
839 printf(" %f", wSigmaTgl2*1000); //15
840 printf(" %f", wSigmaY2); //16
841 // printf(" %f", wSigmaZ2); //17
842 printf(" %f", wChi2); //17
843 printf(" %f", wC); //18
850 return expectedNumberOfClusters;
855 //___________________________________________________________________
857 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
859 // Starting from current radial position of track <t> this function
860 // extrapolates the track up to outer timebin and in the sensitive
861 // layers confirms prolongation if a close cluster is found.
862 // Returns the number of clusters expected to be found in sensitive layers
864 Float_t wIndex, wTB, wChi2;
865 Float_t wYrt, wYclosest, wYcorrect, wYwindow;
866 Float_t wZrt, wZclosest, wZcorrect, wZwindow;
867 Float_t wPx, wPy, wPz, wC;
869 Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
871 Int_t trackIndex = t.GetLabel();
873 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
875 Int_t try_again=fMaxGap;
877 Double_t alpha=t.GetAlpha();
879 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
880 if (alpha < 0. ) alpha += 2.*TMath::Pi();
882 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
884 Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
886 Double_t rad_length, rho, x, dx, y, ymax, z;
888 Bool_t lookForCluster;
890 Int_t expectedNumberOfClusters = 0;
893 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
896 for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB; nr++) {
898 y = t.GetY(); z = t.GetZ();
900 // first propagate to the outer surface of the current time bin
902 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
903 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
905 if(!t.PropagateTo(x,rad_length,rho)) break;
908 ymax = x*TMath::Tan(0.5*alpha);
912 if (!t.Rotate(alpha)) break;
913 if(!t.PropagateTo(x,rad_length,rho)) break;
914 } else if (y <-ymax) {
916 if (!t.Rotate(-alpha)) break;
917 if(!t.PropagateTo(x,rad_length,rho)) break;
919 y = t.GetY(); z = t.GetZ();
921 // now propagate to the middle plane of the next time bin
922 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
924 x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
926 if(!t.PropagateTo(x,rad_length,rho)) break;
930 ymax = x*TMath::Tan(0.5*alpha);
932 if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
936 if (!t.Rotate(alpha)) break;
937 if(!t.PropagateTo(x,rad_length,rho)) break;
938 } else if (y <-ymax) {
940 if (!t.Rotate(-alpha)) break;
941 if(!t.PropagateTo(x,rad_length,rho)) break;
944 // printf("label %d, pl %d, lookForCluster %d \n",
945 // trackIndex, nr+1, lookForCluster);
948 expectedNumberOfClusters++;
950 wIndex = (Float_t) t.GetLabel();
951 wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
953 AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr+1));
954 Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
955 Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
956 if((t.GetSigmaY2() + sy2) < 0) break;
957 Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2);
958 Double_t y=t.GetY(), z=t.GetZ();
962 wYwindow = (Float_t) road;
963 t.GetPxPyPz(Px,Py,Pz);
967 wC = (Float_t) t.GetC();
968 wSigmaC2 = (Float_t) t.GetSigmaC2();
969 wSigmaTgl2 = (Float_t) t.GetSigmaTgl2();
970 wSigmaY2 = (Float_t) t.GetSigmaY2();
971 wSigmaZ2 = (Float_t) t.GetSigmaZ2();
974 if (road>fWideRoad) {
975 if (t.GetNumberOfClusters()>4)
976 cerr<<t.GetNumberOfClusters()
977 <<"FindProlongation warning: Too broad road !\n";
984 Double_t max_chi2=fMaxChi2;
986 wYclosest = 12345678;
987 wYcorrect = 12345678;
988 wZclosest = 12345678;
989 wZcorrect = 12345678;
990 wZwindow = TMath::Sqrt(2.25 * 12 * sz2);
992 // Find the closest correct cluster for debugging purposes
994 Float_t minDY = 1000000;
995 for (Int_t i=0; i<time_bin; i++) {
996 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
997 if((c->GetLabel(0) != trackIndex) &&
998 (c->GetLabel(1) != trackIndex) &&
999 (c->GetLabel(2) != trackIndex)) continue;
1000 if(TMath::Abs(c->GetY() - y) > minDY) continue;
1001 minDY = TMath::Abs(c->GetY() - y);
1002 wYcorrect = c->GetY();
1003 wZcorrect = c->GetZ();
1005 Double_t h01 = GetTiltFactor(c);
1006 wChi2 = t.GetPredictedChi2(c, h01);
1010 // Now go for the real cluster search
1014 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
1015 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1016 if (c->GetY() > y+road) break;
1017 // if (c->IsUsed() > 0) continue;
1018 if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1020 Double_t h01 = GetTiltFactor(c);
1021 Double_t chi2=t.GetPredictedChi2(c,h01);
1023 if (chi2 > max_chi2) continue;
1026 index=time_bin.GetIndex(i);
1031 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
1032 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1034 if (c->GetY() > y+road) break;
1035 // if (c->IsUsed() > 0) continue;
1036 if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1038 Double_t h01 = GetTiltFactor(c);
1039 Double_t chi2=t.GetPredictedChi2(c,h01);
1041 if (chi2 > max_chi2) continue;
1044 index=time_bin.GetIndex(i);
1049 wYclosest = cl->GetY();
1050 wZclosest = cl->GetZ();
1052 t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1053 Double_t h01 = GetTiltFactor(cl);
1054 if(!t.Update(cl,max_chi2,index,h01)) {
1055 if(!try_again--) return 0;
1057 else try_again=fMaxGap;
1060 if (try_again==0) break;
1065 if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1067 printf(" %f", wIndex); //1
1068 printf(" %f", wTB); //2
1069 printf(" %f", wYrt); //3
1070 printf(" %f", wYclosest); //4
1071 printf(" %f", wYcorrect); //5
1072 printf(" %f", wYwindow); //6
1073 printf(" %f", wZrt); //7
1074 printf(" %f", wZclosest); //8
1075 printf(" %f", wZcorrect); //9
1076 printf(" %f", wZwindow); //10
1077 printf(" %f", wPx); //11
1078 printf(" %f", wPy); //12
1079 printf(" %f", wPz); //13
1080 printf(" %f", wSigmaC2*1000000); //14
1081 printf(" %f", wSigmaTgl2*1000); //15
1082 printf(" %f", wSigmaY2); //16
1083 // printf(" %f", wSigmaZ2); //17
1084 printf(" %f", wChi2); //17
1085 printf(" %f", wC); //18
1092 return expectedNumberOfClusters;
1095 //___________________________________________________________________
1097 Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1099 // Starting from current radial position of track <t> this function
1100 // extrapolates the track up to radial position <xToGo>.
1101 // Returns 1 if track reaches the plane, and 0 otherwise
1103 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1105 Double_t alpha=t.GetAlpha();
1107 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1108 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1110 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1112 Bool_t lookForCluster;
1113 Double_t rad_length, rho, x, dx, y, ymax, z;
1117 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1119 Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1121 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
1123 y = t.GetY(); z = t.GetZ();
1125 // first propagate to the outer surface of the current time bin
1126 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1127 x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1128 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1130 ymax = x*TMath::Tan(0.5*alpha);
1133 if (!t.Rotate(alpha)) return 0;
1134 } else if (y <-ymax) {
1136 if (!t.Rotate(-alpha)) return 0;
1138 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1140 y = t.GetY(); z = t.GetZ();
1142 // now propagate to the middle plane of the next time bin
1143 fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1144 x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1145 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1147 ymax = x*TMath::Tan(0.5*alpha);
1150 if (!t.Rotate(alpha)) return 0;
1151 } else if (y <-ymax) {
1153 if (!t.Rotate(-alpha)) return 0;
1155 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1160 //___________________________________________________________________
1162 Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1164 // Starting from current radial position of track <t> this function
1165 // extrapolates the track up to radial position of the outermost
1166 // padrow of the TPC.
1167 // Returns 1 if track reaches the TPC, and 0 otherwise
1169 Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1171 Double_t alpha=t.GetAlpha();
1173 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1174 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1176 Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1178 Bool_t lookForCluster;
1179 Double_t rad_length, rho, x, dx, y, ymax, z;
1183 alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
1185 Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1187 for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1189 y = t.GetY(); z = t.GetZ();
1191 // first propagate to the outer surface of the current time bin
1192 fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1193 x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1194 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1196 ymax = x*TMath::Tan(0.5*alpha);
1199 if (!t.Rotate(alpha)) return 0;
1200 } else if (y <-ymax) {
1202 if (!t.Rotate(-alpha)) return 0;
1204 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1206 y = t.GetY(); z = t.GetZ();
1208 // now propagate to the middle plane of the next time bin
1209 fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1210 x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1211 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1213 ymax = x*TMath::Tan(0.5*alpha);
1216 if (!t.Rotate(alpha)) return 0;
1217 } else if (y <-ymax) {
1219 if (!t.Rotate(-alpha)) return 0;
1221 if(!t.PropagateTo(x,rad_length,rho)) return 0;
1227 //_____________________________________________________________________________
1228 void AliTRDtracker::LoadEvent()
1230 // Fills clusters into TRD tracking_sectors
1231 // Note that the numbering scheme for the TRD tracking_sectors
1232 // differs from that of TRD sectors
1234 ReadClusters(fClusters);
1235 Int_t ncl=fClusters->GetEntriesFast();
1236 cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1240 // printf("\r %d left ",ncl);
1241 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1242 Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
1243 Int_t sector=fGeom->GetSector(detector);
1244 Int_t plane=fGeom->GetPlane(detector);
1246 Int_t tracking_sector = CookSectorIndex(sector);
1248 Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin);
1249 if(gtb < 0) continue;
1250 Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb);
1253 fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index);
1259 //_____________________________________________________________________________
1260 void AliTRDtracker::UnloadEvent()
1263 // Clears the arrays of clusters and tracks. Resets sectors and timebins
1268 nentr = fClusters->GetEntriesFast();
1269 for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1271 nentr = fSeeds->GetEntriesFast();
1272 for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1274 nentr = fTracks->GetEntriesFast();
1275 for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1277 Int_t nsec = AliTRDgeometry::kNsect;
1279 for (i = 0; i < nsec; i++) {
1280 for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1281 fTrSec[i]->GetLayer(pl)->Clear();
1287 //__________________________________________________________________________
1288 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1290 // Creates track seeds using clusters in timeBins=i1,i2
1293 cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1297 Double_t x[5], c[15];
1298 Int_t max_sec=AliTRDgeometry::kNsect;
1300 Double_t alpha=AliTRDgeometry::GetAlpha();
1301 Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1302 Double_t cs=cos(alpha), sn=sin(alpha);
1303 Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1306 Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1307 Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1309 Double_t x1 =fTrSec[0]->GetX(i1);
1310 Double_t xx2=fTrSec[0]->GetX(i2);
1312 for (Int_t ns=0; ns<max_sec; ns++) {
1314 Int_t nl2 = *(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1315 Int_t nl=(*fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1316 Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1317 Int_t nu=(*fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1318 Int_t nu2=(*fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1320 AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1322 for (Int_t is=0; is < r1; is++) {
1323 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1325 for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1327 const AliTRDcluster *cl;
1328 Double_t x2, y2, z2;
1329 Double_t x3=0., y3=0.;
1332 if(turn != 2) continue;
1333 AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1335 y2=cl->GetY(); z2=cl->GetZ();
1340 else if (js<nl2+nl) {
1341 if(turn != 1) continue;
1342 AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1344 y2=cl->GetY(); z2=cl->GetZ();
1349 else if (js<nl2+nl+nm) {
1350 if(turn != 1) continue;
1351 AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1353 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1355 else if (js<nl2+nl+nm+nu) {
1356 if(turn != 1) continue;
1357 AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1358 cl=r2[js-nl2-nl-nm];
1359 y2=cl->GetY(); z2=cl->GetZ();
1365 if(turn != 2) continue;
1366 AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1367 cl=r2[js-nl2-nl-nm-nu];
1368 y2=cl->GetY(); z2=cl->GetZ();
1374 if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
1376 Double_t zz=z1 - z1/x1*(x1-x2);
1378 if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
1380 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1381 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
1385 x[4]=f1trd(x1,y1,x2,y2,x3,y3);
1387 if (TMath::Abs(x[4]) > fMaxSeedC) continue;
1389 x[2]=f2trd(x1,y1,x2,y2,x3,y3);
1391 if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1393 x[3]=f3trd(x1,y1,x2,y2,z1,z2);
1395 if (TMath::Abs(x[3]) > fMaxSeedTan) continue;
1397 Double_t a=asin(x[2]);
1398 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1400 if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
1402 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
1403 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2();
1404 Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
1407 Double_t h01 = GetTiltFactor(r1[is]);
1408 Double_t xu_factor = 100.;
1414 sy1=sy1+sz1*h01*h01;
1415 Double_t syz=sz1*(-h01);
1416 // end of tilt changes
1418 Double_t f40=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1419 Double_t f42=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1420 Double_t f43=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1421 Double_t f20=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1422 Double_t f22=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1423 Double_t f23=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1424 Double_t f30=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1425 Double_t f31=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1426 Double_t f32=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1427 Double_t f34=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1431 // c[1]=0.; c[2]=sz1;
1432 c[1]=syz; c[2]=sz1*xu_factor;
1433 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1434 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
1435 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1436 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1437 c[13]=f30*sy1*f40+f32*sy2*f42;
1438 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1440 UInt_t index=r1.GetIndex(is);
1442 AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
1444 Int_t rc=FollowProlongation(*track, i2);
1447 (track->GetNumberOfClusters() <
1448 (outer-inner)*fMinClustersInSeed)) delete track;
1450 fSeeds->AddLast(track); fNseeds++;
1451 // cerr<<"\r found seed "<<fNseeds;
1458 //_____________________________________________________________________________
1459 void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp)
1462 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
1463 // from the file. The names of the cluster tree and branches
1464 // should match the ones used in AliTRDclusterizer::WriteClusters()
1467 TDirectory *savedir=gDirectory;
1470 TFile *in=(TFile*)inp;
1471 if (!in->IsOpen()) {
1472 cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
1480 Char_t treeName[12];
1481 sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
1482 TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1484 TObjArray *ClusterArray = new TObjArray(400);
1486 ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
1488 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1489 printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
1491 // Loop through all entries in the tree
1493 AliTRDcluster *c = 0;
1496 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1499 nbytes += ClusterTree->GetEvent(iEntry);
1501 // Get the number of points in the detector
1502 Int_t nCluster = ClusterArray->GetEntriesFast();
1503 // printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1505 // Loop through all TRD digits
1506 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1507 c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1508 AliTRDcluster *co = new AliTRDcluster(*c);
1509 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1510 Int_t ltb = co->GetLocalTimeBin();
1511 if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
1512 else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1514 delete ClusterArray->RemoveAt(iCluster);
1518 delete ClusterArray;
1523 //______________________________________________________________________
1524 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
1527 // Reads AliTRDclusters from file <filename>. The names of the cluster
1528 // tree and branches should match the ones used in
1529 // AliTRDclusterizer::WriteClusters()
1530 // if <array> == 0, clusters are added into AliTRDtracker fCluster array
1533 TDirectory *savedir=gDirectory;
1535 TFile *file = TFile::Open(filename);
1536 if (!file->IsOpen()) {
1537 cerr<<"Can't open file with TRD clusters"<<endl;
1541 Char_t treeName[12];
1542 sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
1543 TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1546 cerr<<"AliTRDtracker::ReadClusters(): ";
1547 cerr<<"can't get a tree with clusters !\n";
1551 TObjArray *ClusterArray = new TObjArray(400);
1553 ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
1555 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1556 cout<<"found "<<nEntries<<" in ClusterTree"<<endl;
1558 // Loop through all entries in the tree
1560 AliTRDcluster *c = 0;
1564 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1567 nbytes += ClusterTree->GetEvent(iEntry);
1569 // Get the number of points in the detector
1570 Int_t nCluster = ClusterArray->GetEntriesFast();
1571 printf("\n Read %d clusters from entry %d", nCluster, iEntry);
1573 // Loop through all TRD digits
1574 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1575 c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1576 AliTRDcluster *co = new AliTRDcluster(*c);
1577 co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1578 Int_t ltb = co->GetLocalTimeBin();
1579 if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
1580 else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1582 delete ClusterArray->RemoveAt(iCluster);
1587 delete ClusterArray;
1593 //__________________________________________________________________
1594 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const {
1596 Int_t label=123456789, index, i, j;
1597 Int_t ncl=pt->GetNumberOfClusters();
1598 const Int_t range = fTrSec[0]->GetOuterTimeBin()+1;
1602 // Int_t s[range][2];
1603 Int_t **s = new Int_t* [range];
1604 for (i=0; i<range; i++) {
1605 s[i] = new Int_t[2];
1607 for (i=0; i<range; i++) {
1613 for (i=0; i<ncl; i++) {
1614 index=pt->GetClusterIndex(i);
1615 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1621 for (i=0; i<ncl; i++) {
1622 index=pt->GetClusterIndex(i);
1623 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1624 for (Int_t k=0; k<3; k++) {
1625 label=c->GetLabel(k);
1626 label_added=kFALSE; j=0;
1628 while ( (!label_added) && ( j < range ) ) {
1629 if (s[j][0]==label || s[j][1]==0) {
1643 for (i=0; i<range; i++) {
1645 max=s[i][1]; label=s[i][0];
1649 for (i=0; i<range; i++) {
1655 if ((1.- Float_t(max)/ncl) > wrong) label=-label;
1657 pt->SetLabel(label);
1662 //__________________________________________________________________
1663 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const {
1664 Int_t ncl=t->GetNumberOfClusters();
1665 for (Int_t i=from; i<ncl; i++) {
1666 Int_t index = t->GetClusterIndex(i);
1667 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1673 //_____________________________________________________________________
1674 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
1676 // Parametrised "expected" error of the cluster reconstruction in Y
1678 Double_t s = 0.08 * 0.08;
1682 //_____________________________________________________________________
1683 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
1685 // Parametrised "expected" error of the cluster reconstruction in Z
1687 Double_t s = 9 * 9 /12.;
1692 //_____________________________________________________________________
1693 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const
1696 // Returns radial position which corresponds to time bin <local_tb>
1697 // in tracking sector <sector> and plane <plane>
1700 Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, local_tb);
1701 Int_t pl = fTrSec[sector]->GetLayerNumber(index);
1702 return fTrSec[sector]->GetLayer(pl)->GetX();
1707 //_______________________________________________________
1708 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
1709 Double_t dx, Double_t rho, Double_t rad_length, Int_t tb_index)
1712 // AliTRDpropagationLayer constructor
1715 fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = rad_length;
1716 fClusters = NULL; fIndex = NULL; fTimeBinIndex = tb_index;
1719 for(Int_t i=0; i < (Int_t) kZONES; i++) {
1720 fZc[i]=0; fZmax[i] = 0;
1725 if(fTimeBinIndex >= 0) {
1726 fClusters = new AliTRDcluster*[kMAX_CLUSTER_PER_TIME_BIN];
1727 fIndex = new UInt_t[kMAX_CLUSTER_PER_TIME_BIN];
1740 //_______________________________________________________
1741 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
1742 Double_t Zmax, Double_t Ymax, Double_t rho,
1743 Double_t rad_length, Double_t Yc, Double_t Zc)
1746 // Sets hole in the layer
1755 fHoleX0 = rad_length;
1759 //_______________________________________________________
1760 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
1763 // AliTRDtrackingSector Constructor
1772 for(UInt_t i=0; i < kMAX_TIME_BIN_INDEX; i++) fTimeBinIndex[i] = -1;
1775 AliTRDpropagationLayer* ppl;
1777 Double_t x, xin, xout, dx, rho, rad_length;
1780 // set time bins in the gas of the TPC
1782 xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
1783 rho = 0.9e-3; rad_length = 28.94;
1785 for(Int_t i=0; i<steps; i++) {
1786 x = xin + i*dx + dx/2;
1787 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1791 // set time bins in the outer field cage vessel
1793 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
1794 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1797 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
1798 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1801 dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; rad_length = 41.28; // Nomex
1802 steps = 5; dx = (xout - xin)/steps;
1803 for(Int_t i=0; i<steps; i++) {
1804 x = xin + i*dx + dx/2;
1805 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1809 dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
1810 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1813 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
1814 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1818 // set time bins in CO2
1820 xin = xout; xout = 275.0;
1821 steps = 50; dx = (xout - xin)/steps;
1822 rho = 1.977e-3; rad_length = 36.2;
1824 for(Int_t i=0; i<steps; i++) {
1825 x = xin + i*dx + dx/2;
1826 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1830 // set time bins in the outer containment vessel
1832 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; rad_length = 24.01; // Al
1833 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1836 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
1837 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1840 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
1841 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1844 dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; rad_length = 41.28; // Nomex
1845 steps = 10; dx = (xout - xin)/steps;
1846 for(Int_t i=0; i<steps; i++) {
1847 x = xin + i*dx + dx/2;
1848 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1852 dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
1853 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1856 dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
1857 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1860 dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; rad_length = 24.01; // Al
1861 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1864 Double_t xtrd = (Double_t) fGeom->Rmin();
1866 // add layers between TPC and TRD (Air temporarily)
1867 xin = xout; xout = xtrd;
1868 steps = 50; dx = (xout - xin)/steps;
1869 rho = 1.2e-3; rad_length = 36.66;
1871 for(Int_t i=0; i<steps; i++) {
1872 x = xin + i*dx + dx/2;
1873 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1878 Double_t alpha=AliTRDgeometry::GetAlpha();
1880 // add layers for each of the planes
1882 Double_t dxRo = (Double_t) fGeom->CroHght(); // Rohacell
1883 Double_t dxSpace = (Double_t) fGeom->Cspace(); // Spacing between planes
1884 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
1885 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
1886 Double_t dxRad = (Double_t) fGeom->CraHght(); // Radiator
1887 Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
1888 Double_t dxPlane = dxTEC + dxSpace;
1891 const Int_t nChambers = AliTRDgeometry::Ncham();
1892 Double_t Ymax = 0, holeYmax = 0;
1893 Double_t * Zc = new Double_t[nChambers];
1894 Double_t * Zmax = new Double_t[nChambers];
1895 Double_t holeZmax = 1000.; // the whole sector is missing
1897 for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
1900 xin = xtrd + plane * dxPlane; xout = xin + dxRad;
1901 steps = 12; dx = (xout - xin)/steps; rho = 0.074; rad_length = 40.6;
1902 for(Int_t i=0; i<steps; i++) {
1903 x = xin + i*dx + dx/2;
1904 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1905 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1906 holeYmax = x*TMath::Tan(0.5*alpha);
1907 ppl->SetHole(holeYmax, holeZmax);
1909 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1910 holeYmax = x*TMath::Tan(0.5*alpha);
1911 ppl->SetHole(holeYmax, holeZmax);
1916 Ymax = fGeom->GetChamberWidth(plane)/2;
1917 for(Int_t ch = 0; ch < nChambers; ch++) {
1918 Zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
1919 Float_t pad = fPar->GetRowPadSize(plane,ch,0);
1920 Float_t row0 = fPar->GetRow0(plane,ch,0);
1921 Int_t nPads = fPar->GetRowMax(plane,ch,0);
1922 Zc[ch] = (pad * nPads)/2 + row0 - pad/2;
1925 dx = fPar->GetTimeBinSize();
1926 rho = 0.00295 * 0.85; rad_length = 11.0;
1928 Double_t x0 = (Double_t) fPar->GetTime0(plane);
1929 Double_t xbottom = x0 - dxDrift;
1930 Double_t xtop = x0 + dxAmp;
1932 // Amplification region
1934 steps = (Int_t) (dxAmp/dx);
1936 for(tb = 0; tb < steps; tb++) {
1937 x = x0 + tb * dx + dx/2;
1938 tb_index = CookTimeBinIndex(plane, -tb-1);
1939 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
1941 for(Int_t ch = 0; ch < nChambers; ch++) {
1942 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1944 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1945 holeYmax = x*TMath::Tan(0.5*alpha);
1946 ppl->SetHole(holeYmax, holeZmax);
1948 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1949 holeYmax = x*TMath::Tan(0.5*alpha);
1950 ppl->SetHole(holeYmax, holeZmax);
1954 tb_index = CookTimeBinIndex(plane, -steps);
1955 x = (x + dx/2 + xtop)/2;
1957 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
1959 for(Int_t ch = 0; ch < nChambers; ch++) {
1960 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1962 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1963 holeYmax = x*TMath::Tan(0.5*alpha);
1964 ppl->SetHole(holeYmax, holeZmax);
1966 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1967 holeYmax = x*TMath::Tan(0.5*alpha);
1968 ppl->SetHole(holeYmax, holeZmax);
1973 dx = fPar->GetTimeBinSize();
1974 steps = (Int_t) (dxDrift/dx);
1976 for(tb = 0; tb < steps; tb++) {
1977 x = x0 - tb * dx - dx/2;
1978 tb_index = CookTimeBinIndex(plane, tb);
1980 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
1982 for(Int_t ch = 0; ch < nChambers; ch++) {
1983 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1985 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1986 holeYmax = x*TMath::Tan(0.5*alpha);
1987 ppl->SetHole(holeYmax, holeZmax);
1989 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1990 holeYmax = x*TMath::Tan(0.5*alpha);
1991 ppl->SetHole(holeYmax, holeZmax);
1995 tb_index = CookTimeBinIndex(plane, steps);
1996 x = (x - dx/2 + xbottom)/2;
1998 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
2000 for(Int_t ch = 0; ch < nChambers; ch++) {
2001 ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
2003 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2004 holeYmax = x*TMath::Tan(0.5*alpha);
2005 ppl->SetHole(holeYmax, holeZmax);
2007 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2008 holeYmax = x*TMath::Tan(0.5*alpha);
2009 ppl->SetHole(holeYmax, holeZmax);
2014 xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; rad_length = 33.0;
2015 ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2016 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2017 holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
2018 ppl->SetHole(holeYmax, holeZmax);
2020 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2021 holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
2022 ppl->SetHole(holeYmax, holeZmax);
2027 xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
2028 steps = 5; dx = (xout - xin)/steps; rho = 0.074; rad_length = 40.6;
2029 for(Int_t i=0; i<steps; i++) {
2030 x = xin + i*dx + dx/2;
2031 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2032 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2033 holeYmax = x*TMath::Tan(0.5*alpha);
2034 ppl->SetHole(holeYmax, holeZmax);
2036 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2037 holeYmax = x*TMath::Tan(0.5*alpha);
2038 ppl->SetHole(holeYmax, holeZmax);
2043 // Space between the chambers, air
2044 xin = xout; xout = xtrd + (plane + 1) * dxPlane;
2045 steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; rad_length = 36.66;
2046 for(Int_t i=0; i<steps; i++) {
2047 x = xin + i*dx + dx/2;
2048 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2049 if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2050 holeYmax = x*TMath::Tan(0.5*alpha);
2051 ppl->SetHole(holeYmax, holeZmax);
2053 if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2054 holeYmax = x*TMath::Tan(0.5*alpha);
2055 ppl->SetHole(holeYmax, holeZmax);
2061 // Space between the TRD and RICH
2062 Double_t xRICH = 500.;
2063 xin = xout; xout = xRICH;
2064 steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; rad_length = 36.66;
2065 for(Int_t i=0; i<steps; i++) {
2066 x = xin + i*dx + dx/2;
2067 ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2077 //______________________________________________________
2079 Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t local_tb) const
2082 // depending on the digitization parameters calculates "global"
2083 // time bin index for timebin <local_tb> in plane <plane>
2086 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
2087 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2088 Double_t dx = (Double_t) fPar->GetTimeBinSize();
2090 Int_t tbAmp = fPar->GetTimeBefore();
2091 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
2092 if(kTRUE) maxAmp = 0; // intentional until we change parameter class
2093 Int_t tbDrift = fPar->GetTimeMax();
2094 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2096 Int_t tb_per_plane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
2098 Int_t gtb = (plane+1) * tb_per_plane - local_tb - 1 - TMath::Min(tbAmp,maxAmp);
2100 if((local_tb < 0) &&
2101 (TMath::Abs(local_tb) > TMath::Min(tbAmp,maxAmp))) return -1;
2102 if(local_tb >= TMath::Min(tbDrift,maxDrift)) return -1;
2109 //______________________________________________________
2111 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2114 // For all sensitive time bins sets corresponding layer index
2115 // in the array fTimeBins
2120 for(Int_t i = 0; i < fN; i++) {
2121 index = fLayers[i]->GetTimeBinIndex();
2123 // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2125 if(index < 0) continue;
2126 if(index >= (Int_t) kMAX_TIME_BIN_INDEX) {
2127 printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2128 printf(" index %d exceeds allowed maximum of %d!\n",
2129 index, kMAX_TIME_BIN_INDEX-1);
2132 fTimeBinIndex[index] = i;
2135 Double_t x1, dx1, x2, dx2, gap;
2137 for(Int_t i = 0; i < fN-1; i++) {
2138 x1 = fLayers[i]->GetX();
2139 dx1 = fLayers[i]->GetdX();
2140 x2 = fLayers[i+1]->GetX();
2141 dx2 = fLayers[i+1]->GetdX();
2142 gap = (x2 - dx2/2) - (x1 + dx1/2);
2144 printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2145 printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2148 printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2149 printf(" (%f - %f) - (%f + %f) = %f\n",
2150 x2, dx2/2, x1, dx1, gap);
2156 //______________________________________________________
2159 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2162 // Returns the number of time bin which in radial position is closest to <x>
2165 if(x >= fLayers[fN-1]->GetX()) return fN-1;
2166 if(x <= fLayers[0]->GetX()) return 0;
2168 Int_t b=0, e=fN-1, m=(b+e)/2;
2169 for (; b<e; m=(b+e)/2) {
2170 if (x > fLayers[m]->GetX()) b=m+1;
2173 if(TMath::Abs(x - fLayers[m]->GetX()) >
2174 TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2179 //______________________________________________________
2181 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2184 // Returns number of the innermost SENSITIVE propagation layer
2187 return GetLayerNumber(0);
2190 //______________________________________________________
2192 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2195 // Returns number of the outermost SENSITIVE time bin
2198 return GetLayerNumber(GetNumberOfTimeBins() - 1);
2201 //______________________________________________________
2203 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2206 // Returns number of SENSITIVE time bins
2210 for(tb = kMAX_TIME_BIN_INDEX-1; tb >=0; tb--) {
2211 layer = GetLayerNumber(tb);
2217 //______________________________________________________
2219 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2222 // Insert layer <pl> in fLayers array.
2223 // Layers are sorted according to X coordinate.
2225 if ( fN == ((Int_t) kMAX_LAYERS_PER_SECTOR)) {
2226 printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2229 if (fN==0) {fLayers[fN++] = pl; return;}
2230 Int_t i=Find(pl->GetX());
2232 memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2233 fLayers[i]=pl; fN++;
2237 //______________________________________________________
2239 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
2242 // Returns index of the propagation layer nearest to X
2245 if (x <= fLayers[0]->GetX()) return 0;
2246 if (x > fLayers[fN-1]->GetX()) return fN;
2247 Int_t b=0, e=fN-1, m=(b+e)/2;
2248 for (; b<e; m=(b+e)/2) {
2249 if (x > fLayers[m]->GetX()) b=m+1;
2255 //______________________________________________________
2257 void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
2258 Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &rad_length,
2259 Bool_t &lookForCluster) const
2262 // Returns radial step <dx>, density <rho>, rad. length <rad_length>,
2263 // and sensitivity <lookForCluster> in point <y,z>
2269 lookForCluster = kFALSE;
2271 // check dead regions
2272 if(fTimeBinIndex >= 0) {
2273 for(Int_t ch = 0; ch < (Int_t) kZONES; ch++) {
2274 if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
2275 lookForCluster = kTRUE;
2276 // else { rho = 1.7; rad_length = 33.0; } // G10
2278 if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE;
2279 if(!lookForCluster) {
2280 // rho = 1.7; rad_length = 33.0; // G10
2285 if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) &&
2286 (TMath::Abs(z - fHoleZc) < fHoleZmax)) {
2287 lookForCluster = kFALSE;
2289 rad_length = fHoleX0;
2295 //______________________________________________________
2297 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
2300 // Insert cluster in cluster array.
2301 // Clusters are sorted according to Y coordinate.
2303 if(fTimeBinIndex < 0) {
2304 printf("*** attempt to insert cluster into non-sensitive time bin!\n");
2308 if (fN== (Int_t) kMAX_CLUSTER_PER_TIME_BIN) {
2309 printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
2312 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2313 Int_t i=Find(c->GetY());
2314 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2315 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
2316 fIndex[i]=index; fClusters[i]=c; fN++;
2319 //______________________________________________________
2321 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
2323 // Returns index of the cluster nearest in Y
2325 if (y <= fClusters[0]->GetY()) return 0;
2326 if (y > fClusters[fN-1]->GetY()) return fN;
2327 Int_t b=0, e=fN-1, m=(b+e)/2;
2328 for (; b<e; m=(b+e)/2) {
2329 if (y > fClusters[m]->GetY()) b=m+1;
2335 //---------------------------------------------------------
2337 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
2339 // Returns correction factor for tilted pads geometry
2342 Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
2343 Int_t det = c->GetDetector();
2344 Int_t plane = fGeom->GetPlane(det);
2346 if((plane == 1) || (plane == 3) || (plane == 5)) h01=-h01;
2348 if(fNoTilt) h01 = 0;