]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Services added to ITS
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
e674b993 18Revision 1.14 1999/09/29 09:24:33 fca
19Introduction of the Copyright and cvs Log
20
4c039060 21*/
22
fe4da5cc 23///////////////////////////////////////////////////////////////////////////////
24// //
25// Time Projection Chamber //
26// This class contains the basic functions for the Time Projection Chamber //
27// detector. Functions specific to one particular geometry are //
28// contained in the derived classes //
29// //
30//Begin_Html
31/*
1439f98e 32<img src="picts/AliTPCClass.gif">
fe4da5cc 33*/
34//End_Html
35// //
8c555625 36// //
fe4da5cc 37///////////////////////////////////////////////////////////////////////////////
38
39#include <TMath.h>
40#include <TRandom.h>
41#include <TVector.h>
8c555625 42#include <TMatrix.h>
fe4da5cc 43#include <TGeometry.h>
44#include <TNode.h>
45#include <TTUBS.h>
46#include <TObjectTable.h>
1578254f 47#include "TParticle.h"
fe4da5cc 48#include "AliTPC.h"
49#include "AliRun.h"
50#include <iostream.h>
51#include <fstream.h>
52#include "AliMC.h"
53
8c555625 54//MI change
55#include "AliTPCParam.h"
56#include "AliTPCD.h"
57#include "AliTPCPRF2D.h"
58#include "AliTPCRF1D.h"
59
60
1283eee5 61
fe4da5cc 62ClassImp(AliTPC)
63
64//_____________________________________________________________________________
65AliTPC::AliTPC()
66{
67 //
68 // Default constructor
69 //
70 fIshunt = 0;
71 fClusters = 0;
72 fHits = 0;
73 fDigits = 0;
74 fTracks = 0;
75 fNsectors = 0;
76 fNtracks = 0;
77 fNclusters= 0;
3c0f9266 78
8c555625 79 fDigParam= new AliTPCD();
80 fDigits = fDigParam->GetArray();
fe4da5cc 81}
82
83//_____________________________________________________________________________
84AliTPC::AliTPC(const char *name, const char *title)
85 : AliDetector(name,title)
86{
87 //
88 // Standard constructor
89 //
90
91 //
92 // Initialise arrays of hits and digits
93 fHits = new TClonesArray("AliTPChit", 176);
8c555625 94 // fDigits = new TClonesArray("AliTPCdigit",10000);
95 //MI change
96 fDigParam= new AliTPCD;
97 fDigits = fDigParam->GetArray();
3c0f9266 98
99 AliTPCParam *fTPCParam = &(fDigParam->GetParam());
100
fe4da5cc 101 //
102 // Initialise counters
3c0f9266 103 //
fe4da5cc 104 fClusters = 0;
105 fTracks = 0;
3c0f9266 106 fNsectors = fTPCParam->GetNSector();
fe4da5cc 107 fNtracks = 0;
108 fNclusters= 0;
109 fDigitsIndex = new Int_t[fNsectors+1];
fe4da5cc 110 //
111 fIshunt = 0;
112 //
113 // Initialise color attributes
114 SetMarkerColor(kYellow);
115}
116
117//_____________________________________________________________________________
118AliTPC::~AliTPC()
119{
120 //
121 // TPC destructor
122 //
123 fIshunt = 0;
124 delete fHits;
125 delete fDigits;
126 delete fClusters;
127 delete fTracks;
8c555625 128 delete fDigParam;
fe4da5cc 129 if (fDigitsIndex) delete [] fDigitsIndex;
fe4da5cc 130}
131
132//_____________________________________________________________________________
133void AliTPC::AddCluster(Float_t *hits, Int_t *tracks)
134{
135 //
136 // Add a simulated cluster to the list
137 //
138 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
139 TClonesArray &lclusters = *fClusters;
140 new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks);
141}
142
143//_____________________________________________________________________________
144void AliTPC::AddCluster(const AliTPCcluster &c)
145{
146 //
147 // Add a simulated cluster copy to the list
148 //
149 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
150 TClonesArray &lclusters = *fClusters;
151 new(lclusters[fNclusters++]) AliTPCcluster(c);
152}
153
154//_____________________________________________________________________________
155void AliTPC::AddDigit(Int_t *tracks, Int_t *digits)
156{
157 //
158 // Add a TPC digit to the list
159 //
8c555625 160 // TClonesArray &ldigits = *fDigits;
161 //MI change
162 TClonesArray &ldigits = *fDigParam->GetArray();
fe4da5cc 163 new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits);
164}
165
166//_____________________________________________________________________________
167void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
168{
169 //
170 // Add a hit to the list
171 //
172 TClonesArray &lhits = *fHits;
173 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
174}
175
176//_____________________________________________________________________________
177void AliTPC::AddTrack(Float_t *hits)
178{
179 //
180 // Add a track to the list of tracks
181 //
182 TClonesArray &ltracks = *fTracks;
183 new(ltracks[fNtracks++]) AliTPCtrack(hits);
184}
185
186//_____________________________________________________________________________
187void AliTPC::AddTrack(const AliTPCtrack& t)
188{
189 //
190 // Add a track copy to the list of tracks
191 //
192 if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000);
193 TClonesArray &ltracks = *fTracks;
194 new(ltracks[fNtracks++]) AliTPCtrack(t);
195}
196
197//_____________________________________________________________________________
198void AliTPC::BuildGeometry()
199{
200 //
201 // Build TPC ROOT TNode geometry for the event display
202 //
203 TNode *Node, *Top;
204 TTUBS *tubs;
205 Int_t i;
206 const int kColorTPC=19;
1283eee5 207 char name[5], title[25];
fe4da5cc 208 const Double_t kDegrad=TMath::Pi()/180;
1283eee5 209 const Double_t kRaddeg=180./TMath::Pi();
210
211 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
212
213 Float_t InnerOpenAngle = fTPCParam->GetInnerAngle();
214 Float_t OuterOpenAngle = fTPCParam->GetOuterAngle();
215
216 Float_t InnerAngleShift = fTPCParam->GetInnerAngleShift();
217 Float_t OuterAngleShift = fTPCParam->GetOuterAngleShift();
218
219 Int_t nLo = fTPCParam->GetNInnerSector()/2;
220 Int_t nHi = fTPCParam->GetNOuterSector()/2;
221
222 const Double_t loAng = (Double_t)TMath::Nint(InnerOpenAngle*kRaddeg);
223 const Double_t hiAng = (Double_t)TMath::Nint(OuterOpenAngle*kRaddeg);
224 const Double_t loAngSh = (Double_t)TMath::Nint(InnerAngleShift*kRaddeg);
225 const Double_t hiAngSh = (Double_t)TMath::Nint(OuterAngleShift*kRaddeg);
226
227
fe4da5cc 228 const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad);
229 const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad);
1283eee5 230
231 Double_t rl,ru;
232
233
fe4da5cc 234 //
235 // Get ALICE top node
fe4da5cc 236 //
1283eee5 237
238 Top=gAlice->GetGeometry()->GetNode("alice");
239
240 // inner sectors
241
242 rl = fTPCParam->GetInSecLowEdge();
243 ru = fTPCParam->GetInSecUpEdge();
244
245
fe4da5cc 246 for(i=0;i<nLo;i++) {
247 sprintf(name,"LS%2.2d",i);
1283eee5 248 name[4]='\0';
249 sprintf(title,"TPC low sector %3d",i);
250 title[24]='\0';
251
252 tubs = new TTUBS(name,title,"void",rl*loCorr,ru*loCorr,250.,
253 loAng*(i-0.5)+loAngSh,loAng*(i+0.5)+loAngSh);
fe4da5cc 254 tubs->SetNumberOfDivisions(1);
255 Top->cd();
256 Node = new TNode(name,title,name,0,0,0,"");
257 Node->SetLineColor(kColorTPC);
258 fNodes->Add(Node);
259 }
1283eee5 260
fe4da5cc 261 // Outer sectors
1283eee5 262
263 rl = fTPCParam->GetOuSecLowEdge();
264 ru = fTPCParam->GetOuSecUpEdge();
265
fe4da5cc 266 for(i=0;i<nHi;i++) {
267 sprintf(name,"US%2.2d",i);
1283eee5 268 name[4]='\0';
fe4da5cc 269 sprintf(title,"TPC upper sector %d",i);
1283eee5 270 title[24]='\0';
271 tubs = new TTUBS(name,title,"void",rl*hiCorr,ru*hiCorr,250,
272 hiAng*(i-0.5)+hiAngSh,hiAng*(i+0.5)+hiAngSh);
fe4da5cc 273 tubs->SetNumberOfDivisions(1);
274 Top->cd();
275 Node = new TNode(name,title,name,0,0,0,"");
276 Node->SetLineColor(kColorTPC);
277 fNodes->Add(Node);
278 }
1283eee5 279}
280
281
282
fe4da5cc 283//_____________________________________________________________________________
284Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
285{
286 //
287 // Calculate distance from TPC to mouse on the display
288 // Dummy procedure
289 //
290 return 9999;
291}
292
fe4da5cc 293//_____________________________________________________________________________
294static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
295{
296 //
3c0f9266 297 // Parametrised error of the cluster reconstruction (pad direction)
fe4da5cc 298 //
299 pt=TMath::Abs(pt)*1000.;
300 Double_t x=r/pt;
301 tgl=TMath::Abs(tgl);
302 Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x;
303 if (s<0.4e-3) s=0.4e-3;
3c0f9266 304 s*=1.3; //Iouri Belikov
fe4da5cc 305 return s;
306}
307
308//_____________________________________________________________________________
309static Double_t SigmaZ2(Double_t r, Double_t tgl)
310{
311 //
3c0f9266 312 // Parametrised error of the cluster reconstruction (drift direction)
fe4da5cc 313 //
314 tgl=TMath::Abs(tgl);
315 Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl;
316 if (s<0.4e-3) s=0.4e-3;
3c0f9266 317 s*=1.3; //Iouri Belikov
fe4da5cc 318 return s;
319}
320
321//_____________________________________________________________________________
3c0f9266 322inline Double_t f1(Double_t x1,Double_t y1,
fe4da5cc 323 Double_t x2,Double_t y2,
324 Double_t x3,Double_t y3)
325{
3c0f9266 326 //-----------------------------------------------------------------
327 // Initial approximation of the track curvature
fe4da5cc 328 //
3c0f9266 329 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
330 //-----------------------------------------------------------------
fe4da5cc 331 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
332 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
333 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
334 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
335 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
8c555625 336
fe4da5cc 337 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
8c555625 338
fe4da5cc 339 return -xr*yr/sqrt(xr*xr+yr*yr);
340}
341
342
343//_____________________________________________________________________________
3c0f9266 344inline Double_t f2(Double_t x1,Double_t y1,
fe4da5cc 345 Double_t x2,Double_t y2,
346 Double_t x3,Double_t y3)
347{
3c0f9266 348 //-----------------------------------------------------------------
349 // Initial approximation of the track curvature times center of curvature
fe4da5cc 350 //
3c0f9266 351 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
352 //-----------------------------------------------------------------
fe4da5cc 353 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
354 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
355 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
356 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
357 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
8c555625 358
fe4da5cc 359 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
360
361 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
362}
363
364//_____________________________________________________________________________
3c0f9266 365inline Double_t f3(Double_t x1,Double_t y1,
fe4da5cc 366 Double_t x2,Double_t y2,
367 Double_t z1,Double_t z2)
368{
3c0f9266 369 //-----------------------------------------------------------------
370 // Initial approximation of the tangent of the track dip angle
fe4da5cc 371 //
3c0f9266 372 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
373 //-----------------------------------------------------------------
fe4da5cc 374 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
375}
376
377//_____________________________________________________________________________
3c0f9266 378static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec,
379 int s, int rf=0)
fe4da5cc 380{
3c0f9266 381 //-----------------------------------------------------------------
382 // This function tries to find a track prolongation.
fe4da5cc 383 //
3c0f9266 384 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
385 //-----------------------------------------------------------------
386 const int ROWS_TO_SKIP=int(0.5*sec->GetNRows());
387 const Float_t MAX_CHI2=12.;
fe4da5cc 388 int try_again=ROWS_TO_SKIP;
389 Double_t alpha=sec->GetAlpha();
3c0f9266 390 int ns=int(2*TMath::Pi()/alpha+0.5);
8c555625 391
3c0f9266 392 for (int nr=sec->GetRowNumber(t.GetX())-1; nr>=rf; nr--) {
393 Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
394 if (!t.PropagateTo(x)) return 0;
8c555625 395
3c0f9266 396 AliTPCcluster *cl=0;
fe4da5cc 397 Double_t max_chi2=MAX_CHI2;
398 const AliTPCRow& row=sec[s][nr];
399 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
400 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
3c0f9266 401 Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
8c555625 402
fe4da5cc 403 if (road>30) {
3c0f9266 404 if (t>4) cerr<<t<<" FindProlongation warning: Too broad road !\n";
405 return 0;
fe4da5cc 406 }
8c555625 407
408 if (row) {
fe4da5cc 409 for (int i=row.Find(y-road); i<row; i++) {
410 AliTPCcluster* c=(AliTPCcluster*)(row[i]);
411 if (c->fY > y+road) break;
412 if (c->IsUsed()) continue;
3c0f9266 413 if ((c->fZ - z)*(c->fZ - z) > 25.*(t.GetSigmaZ2() + sz2)) continue;
fe4da5cc 414 Double_t chi2=t.GetPredictedChi2(c);
415 if (chi2 > max_chi2) continue;
416 max_chi2=chi2;
417 cl=c;
418 }
8c555625 419 }
fe4da5cc 420 if (cl) {
421 t.Update(cl,max_chi2);
3c0f9266 422 Double_t ll=TMath::Sqrt((1+t.GetTgl()*t.GetTgl())/
423 (1-(t.GetC()*x-t.GetEta())*(t.GetC()*x-t.GetEta())));
424 cl->fdEdX = cl->fQ/ll;
fe4da5cc 425 try_again=ROWS_TO_SKIP;
426 } else {
8c555625 427 if (try_again==0) break;
428 if (y > ymax) {
429 s = (s+1) % ns;
3c0f9266 430 if (!t.Rotate(alpha)) return 0;
8c555625 431 } else if (y <-ymax) {
3c0f9266 432 s = (s-1+ns) % ns;
433 if (!t.Rotate(-alpha)) return 0;
fe4da5cc 434 }
8c555625 435 try_again--;
fe4da5cc 436 }
437 }
8c555625 438
3c0f9266 439 return 1;
fe4da5cc 440}
441
8c555625 442
fe4da5cc 443//_____________________________________________________________________________
3c0f9266 444static void MakeSeeds(TObjArray& seeds,const AliTPCSector *sec, int max_sec,
445int i1, int i2)
fe4da5cc 446{
3c0f9266 447 //-----------------------------------------------------------------
448 // This function creates track seeds.
fe4da5cc 449 //
3c0f9266 450 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
451 //-----------------------------------------------------------------
fe4da5cc 452 TMatrix C(5,5); TVector x(5);
3c0f9266 453 double alpha=sec->GetAlpha(), shift=sec->GetAlphaShift();
454 double cs=cos(alpha), sn=sin(alpha);
fe4da5cc 455 for (int ns=0; ns<max_sec; ns++) {
456 int nl=sec[(ns-1+max_sec)%max_sec][i2];
457 int nm=sec[ns][i2];
458 int nu=sec[(ns+1)%max_sec][i2];
fe4da5cc 459 const AliTPCRow& r1=sec[ns][i1];
460 for (int is=0; is < r1; is++) {
3c0f9266 461 double x1=sec->GetX(i1), y1=r1[is]->fY, z1=r1[is]->fZ;
fe4da5cc 462 for (int js=0; js < nl+nm+nu; js++) {
463 const AliTPCcluster *cl;
8c555625 464 int ks;
3c0f9266 465 double x2=sec->GetX(i2), y2, z2, tmp;
466
fe4da5cc 467 if (js<nl) {
8c555625 468 ks=(ns-1+max_sec)%max_sec;
fe4da5cc 469 const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2];
470 cl=r2[js];
3c0f9266 471 y2=cl->fY; z2=cl->fZ;
472 tmp= x2*cs+y2*sn;
473 y2 =-x2*sn+y2*cs; x2=tmp;
fe4da5cc 474 } else
475 if (js<nl+nm) {
8c555625 476 ks=ns;
fe4da5cc 477 const AliTPCRow& r2=sec[ns][i2];
478 cl=r2[js-nl];
3c0f9266 479 y2=cl->fY; z2=cl->fZ;
fe4da5cc 480 } else {
8c555625 481 ks=(ns+1)%max_sec;
fe4da5cc 482 const AliTPCRow& r2=sec[(ns+1)%max_sec][i2];
483 cl=r2[js-nl-nm];
3c0f9266 484 y2=cl->fY; z2=cl->fZ;
485 tmp=x2*cs-y2*sn;
486 y2 =x2*sn+y2*cs; x2=tmp;
fe4da5cc 487 }
3c0f9266 488
489 double d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
490 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
491
fe4da5cc 492 x(0)=y1;
493 x(1)=z1;
494 x(2)=f1(x1,y1,x2,y2,0.,0.);
495 x(3)=f2(x1,y1,x2,y2,0.,0.);
496 x(4)=f3(x1,y1,x2,y2,z1,z2);
497
498 if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue;
499
500 if (TMath::Abs(x(4)) > 1.2) continue;
8c555625 501
fe4da5cc 502 Double_t a=asin(x(3));
fe4da5cc 503 Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3)));
504 if (TMath::Abs(zv)>33.) continue;
8c555625 505
fe4da5cc 506 TMatrix X(6,6); X=0.;
507 X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2;
508 X(2,2)=cl->fSigmaY2; X(3,3)=cl->fSigmaZ2;
509 X(4,4)=3./12.; X(5,5)=3./12.;
510 TMatrix F(5,6); F.UnitMatrix();
511 Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1));
512 F(2,0)=(f1(x1,y1+sy,x2,y2,0.,0.)-x(2))/sy;
513 F(2,2)=(f1(x1,y1,x2,y2+sy,0.,0.)-x(2))/sy;
514 F(2,4)=(f1(x1,y1,x2,y2,0.,0.+sy)-x(2))/sy;
515 F(3,0)=(f2(x1,y1+sy,x2,y2,0.,0.)-x(3))/sy;
516 F(3,2)=(f2(x1,y1,x2,y2+sy,0.,0.)-x(3))/sy;
517 F(3,4)=(f2(x1,y1,x2,y2,0.,0.+sy)-x(3))/sy;
518 F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy;
519 F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz;
520 F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy;
521 F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz;
522 F(4,4)=0;
523 F(3,3)=0;
524
525 TMatrix t(F,TMatrix::kMult,X);
526 C.Mult(t,TMatrix(TMatrix::kTransposed,F));
8c555625 527
3c0f9266 528 AliTPCtrack *track=new AliTPCtrack(r1[is], x, C, x1, ns*alpha+shift);
529 int rc=FindProlongation(*track,sec,ns,i2);
8c555625 530 if (rc<0 || *track<(i1-i2)/2) delete track;
531 else seeds.AddLast(track);
fe4da5cc 532 }
533 }
534 }
535}
536
537//_____________________________________________________________________________
3c0f9266 538AliTPCParam *AliTPCSector::param;
fe4da5cc 539void AliTPC::Clusters2Tracks()
540{
3c0f9266 541 //-----------------------------------------------------------------
542 // This is a track finder.
fe4da5cc 543 //
3c0f9266 544 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
545 //-----------------------------------------------------------------
fe4da5cc 546 if (!fClusters) return;
8c555625 547
548 AliTPCParam *p=&fDigParam->GetParam();
3c0f9266 549 AliTPCSector::SetParam(p);
8c555625 550
3c0f9266 551 const int nis=p->GetNInnerSector()/2;
552 AliTPCSSector *ssec=new AliTPCSSector[nis];
553 int nrow_low=ssec->GetNRows();
8c555625 554
3c0f9266 555 const int nos=p->GetNOuterSector()/2;
556 AliTPCLSector *lsec=new AliTPCLSector[nos];
557 int nrow_up=lsec->GetNRows();
8c555625 558
fe4da5cc 559 int ncl=fClusters->GetEntriesFast();
560 while (ncl--) {
561 AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl);
3c0f9266 562 Int_t sec=c->fSector, row=c->fPadRow;
563 if (sec<nis*2) {
564 ssec[sec%nis][row].InsertCluster(c);
fe4da5cc 565 } else {
3c0f9266 566 sec -= nis*2;
567 lsec[sec%nos][row].InsertCluster(c);
fe4da5cc 568 }
569 }
570
fe4da5cc 571 TObjArray seeds(20000);
3c0f9266 572
573 int nrows=nrow_low+nrow_up;
574 int gap=int(0.125*nrows), shift=int(0.5*gap);
575 MakeSeeds(seeds, lsec, nos, nrow_up-1, nrow_up-1-gap);
576 MakeSeeds(seeds, lsec, nos, nrow_up-1-shift, nrow_up-1-shift-gap);
8c555625 577
fe4da5cc 578 seeds.Sort();
579
580 int found=0;
581 int nseed=seeds.GetEntriesFast();
582
583 for (int s=0; s<nseed; s++) {
584 AliTPCtrack& t=*((AliTPCtrack*)seeds.UncheckedAt(s));
3c0f9266 585 double alpha=t.GetAlpha();
fe4da5cc 586 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
587 if (alpha < 0. ) alpha += 2.*TMath::Pi();
3c0f9266 588 int ns=int(alpha/lsec->GetAlpha())%nos;
589
590 if (!FindProlongation(t,lsec,ns)) continue;
591
592 alpha=t.GetAlpha() + 0.5*ssec->GetAlpha() - ssec->GetAlphaShift();
593 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
594 if (alpha < 0. ) alpha += 2.*TMath::Pi();
595 ns=int(alpha/ssec->GetAlpha())%nis; //index of the inner sector needed
596
597 alpha=ns*ssec->GetAlpha() - t.GetAlpha();
598 if (!t.Rotate(alpha)) continue;
599
600 if (!FindProlongation(t,ssec,ns)) continue;
fe4da5cc 601
3c0f9266 602 if (t < int(0.4*nrows)) continue;
fe4da5cc 603
604 AddTrack(t);
605 t.UseClusters();
606 cerr<<found++<<'\r';
607 }
3c0f9266 608
609 delete[] ssec;
610 delete[] lsec;
fe4da5cc 611}
612
613//_____________________________________________________________________________
614void AliTPC::CreateMaterials()
615{
8c555625 616 //-----------------------------------------------
fe4da5cc 617 // Create Materials for for TPC
8c555625 618 //-----------------------------------------------
619
620 //-----------------------------------------------------------------
621 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
622 //-----------------------------------------------------------------
1283eee5 623
fe4da5cc 624 Int_t ISXFLD=gAlice->Field()->Integ();
625 Float_t SXMGMX=gAlice->Field()->Max();
1283eee5 626
627 Float_t amat[5]; // atomic numbers
628 Float_t zmat[5]; // z
629 Float_t wmat[5]; // proportions
630
631 Float_t density;
632
633 // ********************* Gases *******************
634
635 //--------------------------------------------------------------
636 // pure gases
637 //--------------------------------------------------------------
638
639 // Ne
640
641
642 Float_t a_ne = 20.18;
643 Float_t z_ne = 10.;
fe4da5cc 644
1283eee5 645 density = 0.0009;
646
647 AliMaterial(20,"Ne",a_ne,z_ne,density,999.,999.);
648
649 // Ar
650
651 Float_t a_ar = 39.948;
652 Float_t z_ar = 18.;
653
654 density = 0.001782;
655
656 AliMaterial(21,"Ar",a_ar,z_ar,density,999.,999.);
657
658 Float_t a_pure[2];
fe4da5cc 659
1283eee5 660 a_pure[0] = a_ne;
661 a_pure[1] = a_ar;
fe4da5cc 662
1283eee5 663
664 //--------------------------------------------------------------
665 // gases - compounds
666 //--------------------------------------------------------------
667
668 Float_t amol[3];
669
670 // CO2
671
672 amat[0]=12.011;
673 amat[1]=15.9994;
674
675 zmat[0]=6.;
676 zmat[1]=8.;
677
678 wmat[0]=1.;
679 wmat[1]=2.;
680
681 density=0.001977;
682
683 amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
684
685 AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
686
687 // CF4
688
689 amat[0]=12.011;
690 amat[1]=18.998;
691
692 zmat[0]=6.;
693 zmat[1]=9.;
694
695 wmat[0]=1.;
696 wmat[1]=4.;
697
698 density=0.003034;
699
700 amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
701
702 AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
703
704 // CH4
705
706 amat[0]=12.011;
707 amat[1]=1.;
708
709 zmat[0]=6.;
710 zmat[1]=1.;
711
712 wmat[0]=1.;
713 wmat[1]=4.;
714
715 density=0.000717;
716
717 amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
718
719 AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
720
721 //----------------------------------------------------------------
722 // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
723 //----------------------------------------------------------------
724
fe4da5cc 725 char namate[21];
1283eee5 726
727 density = 0.;
728 Float_t am=0;
729 Int_t nc;
730
731 Float_t a,z,rho,absl,X0,buf[1];
732 Int_t nbuf;
733
734 for(nc = 0;nc<fNoComp;nc++)
735 {
736
737 // retrive material constants
738
739 gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
740
741 amat[nc] = a;
742 zmat[nc] = z;
743
744 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
745
746 am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc]);
747 density += fMixtProp[nc]*rho; // density of the mixture
748
749 }
750
751 // mixture proportions by weight!
752
753 for(nc = 0;nc<fNoComp;nc++)
754 {
755
756 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
757
758 wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc])/am;
759
760 }
fe4da5cc 761
1283eee5 762 AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
763 AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
764 AliMixture(33,"Drift gas 3",amat,zmat,density,fNoComp,wmat);
765
766 AliMedium(2, "Drift gas 1", 31, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
767 AliMedium(3, "Drift gas 2", 32, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
768 AliMedium(4, "Drift gas 3", 33, 1, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
769
770 // Air
771
772 AliMaterial(24, "Air", 14.61, 7.3, .001205, 30420., 67500.);
773
774 AliMedium(24, "Air", 24, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
775
776 //----------------------------------------------------------------------
777 // solid materials
778 //----------------------------------------------------------------------
779
780 // Al
781
782 AliMaterial(30, "Al", 26.98, 13., 2.7, 8.9, 37.2);
783
784 AliMedium(0, "Al",30, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
785
786 // Si
787
788 AliMaterial(31, "Si", 28.086, 14.,2.33, 9.36, 999.);
789
790 AliMedium(7, "Al",31, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
fe4da5cc 791
1283eee5 792
793 // Mylar C5H4O2
794
795 amat[0]=12.011;
796 amat[1]=1.;
797 amat[2]=15.9994;
798
799 zmat[0]=6.;
800 zmat[1]=1.;
801 zmat[2]=8.;
802
803 wmat[0]=5.;
804 wmat[1]=4.;
805 wmat[2]=2.;
806
807 density = 1.39;
fe4da5cc 808
1283eee5 809 AliMixture(32, "Mylar",amat,zmat,density,-3,wmat);
810
811 AliMedium(5, "Mylar",32, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
812
813
814
815
816 // Carbon (normal)
817
818 AliMaterial(33,"C normal",12.011,6.,2.265,18.8,999.);
819
820 AliMedium(6,"C normal",33,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
821
822 // G10 for inner and outr field cage
823 // G10 is 60% SiO2 + 40% epoxy, right now I use A and Z for SiO2
824
825 Float_t rhoFactor;
826
827 amat[0]=28.086;
828 amat[1]=15.9994;
829
830 zmat[0]=14.;
831 zmat[1]=8.;
832
833 wmat[0]=1.;
834 wmat[1]=2.;
835
836 density = 1.7;
fe4da5cc 837
1283eee5 838
839 AliMixture(34,"G10 aux.",amat,zmat,density,-2,wmat);
840
841
842 gMC->Gfmate((*fIdmate)[34],namate,a,z,rho,X0,absl,buf,nbuf);
843
844 Float_t thickX0 = 0.0052; // field cage in X0 units
fe4da5cc 845
1283eee5 846 Float_t thick = 2.; // in cm
847
848 X0=19.4; // G10
849
850 rhoFactor = X0*thickX0/thick;
851 density = rho*rhoFactor;
852
853 AliMaterial(35,"G10-fc",a,z,density,999.,999.);
854
855 AliMedium(8,"G10-fc",35,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
856
857 thickX0 = 0.0027; // inner vessel (eta <0.9)
858 thick=0.5;
859 rhoFactor = X0*thickX0/thick;
860 density = rho*rhoFactor;
861
862 AliMaterial(36,"G10-iv",a,z,density,999.,999.);
863
864 AliMedium(9,"G10-iv",36,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
865
866 // Carbon fibre
fe4da5cc 867
1283eee5 868 gMC->Gfmate((*fIdmate)[33],namate,a,z,rho,X0,absl,buf,nbuf);
869
870 thickX0 = 0.0133; // outer vessel
871 thick=3.0;
872 rhoFactor = X0*thickX0/thick;
873 density = rho*rhoFactor;
874
875
876 AliMaterial(37,"C-ov",a,z,density,999.,999.);
877
878 AliMedium(10,"C-ov",37,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
879
880 thickX0=0.015; // inner vessel (cone, eta > 0.9)
881 thick=1.5;
882 rhoFactor = X0*thickX0/thick;
883 density = rho*rhoFactor;
884
885 AliMaterial(38,"C-ivc",a,z,density,999.,999.);
886
887 AliMedium(11,"C-ivc",38,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
888
889 //
890
891 AliMedium(12,"CO2",10,0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
892
893
894
fe4da5cc 895}
896
897//_____________________________________________________________________________
898struct Bin {
899 const AliTPCdigit *dig;
900 int idx;
901 Bin() {dig=0; idx=-1;}
902};
903
904struct PreCluster : public AliTPCcluster {
3c0f9266 905 const AliTPCdigit* summit; //pointer to the maximum digit of this precluster
906 int idx; //index in AliTPC::fClusters
907 int npeaks; //number of peaks in this precluster
908 int ndigits; //number of digits in this precluster
909 PreCluster();
fe4da5cc 910};
3c0f9266 911PreCluster::PreCluster() : AliTPCcluster() {npeaks=ndigits=0;}
fe4da5cc 912
913//_____________________________________________________________________________
3c0f9266 914static void FindPreCluster(int i,int j,int maxj,Bin *bins,PreCluster &c)
fe4da5cc 915{
3c0f9266 916 //-----------------------------------------------------------------
917 // This function looks for "preclusters".
fe4da5cc 918 //
3c0f9266 919 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
920 //-----------------------------------------------------------------
921 Bin& b=bins[i*maxj+j];
922 double q=(double)TMath::Abs(b.dig->fSignal);
8c555625 923
fe4da5cc 924 if (b.idx >= 0 && b.idx != c.idx) {
925 c.idx=b.idx;
926 c.npeaks++;
927 }
928
929 if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
930
931 c.fY += i*q;
932 c.fZ += j*q;
933 c.fSigmaY2 += i*i*q;
934 c.fSigmaZ2 += j*j*q;
935 c.fQ += q;
3c0f9266 936 c.ndigits++;
8c555625 937
fe4da5cc 938 b.dig = 0; b.idx = c.idx;
939
3c0f9266 940 if (bins[(i-1)*maxj+j].dig) FindPreCluster(i-1,j,maxj,bins,c);
941 if (bins[i*maxj+(j-1)].dig) FindPreCluster(i,j-1,maxj,bins,c);
942 if (bins[(i+1)*maxj+j].dig) FindPreCluster(i+1,j,maxj,bins,c);
943 if (bins[i*maxj+(j+1)].dig) FindPreCluster(i,j+1,maxj,bins,c);
fe4da5cc 944
945}
946
947//_____________________________________________________________________________
948void AliTPC::Digits2Clusters()
949{
3c0f9266 950 //-----------------------------------------------------------------
951 // This is a simple cluster finder.
fe4da5cc 952 //
3c0f9266 953 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
954 //-----------------------------------------------------------------
955 AliTPCParam *par = &(fDigParam->GetParam());
8c555625 956
3c0f9266 957 int inp=par->GetNPads(0, par->GetNRowLow()-1);
958 int onp=par->GetNPads(par->GetNSector()-1,par->GetNRowUp() -1);
959 const int MAXY=(inp>onp) ? inp+2 : onp+2;
960 const int MAXTBKT=int((z_end+6*par->GetZSigma())/par->GetZWidth())+1;
961 const int MAXZ=MAXTBKT+2;
962 const int THRESHOLD=20;
fe4da5cc 963
8c555625 964 TTree *t=(TTree*)gDirectory->Get("TreeD0_Param1");
965 t->GetBranch("Digits")->SetAddress(&fDigits);
fe4da5cc 966 Int_t sectors_by_rows=(Int_t)t->GetEntries();
967
968 int ncls=0;
969
970 for (Int_t n=0; n<sectors_by_rows; n++) {
971 if (!t->GetEvent(n)) continue;
3c0f9266 972 Bin *bins=new Bin[MAXY*MAXZ];
fe4da5cc 973 AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
3c0f9266 974 int sec=dig->fSector, row=dig->fPadRow;
975 int ndigits=fDigits->GetEntriesFast();
fe4da5cc 976
3c0f9266 977 int npads, sign;
978 {
979 int nis=par->GetNInnerSector(), nos=par->GetNOuterSector();
980 if (sec < nis) {
981 npads = par->GetNPadsLow(row);
982 sign = (sec < nis/2) ? 1 : -1;
983 } else {
984 npads = par->GetNPadsUp(row);
985 sign = ((sec-nis) < nos/2) ? 1 : -1;
986 }
fe4da5cc 987 }
3c0f9266 988
fe4da5cc 989 int ndig;
990 for (ndig=0; ndig<ndigits; ndig++) {
991 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
3c0f9266 992 int i=dig->fPad+1, j=dig->fTime+1;
993 if (i > npads) {
994 cerr<<"AliTPC::Digits2Clusters error: pad number is out of range ! ";
995 cerr<<i<<' '<<npads<<endl;
996 continue;
997 }
998 if (j > MAXTBKT) {
999 cerr<<"AliTPC::Digits2Clusters error: time bucket is out of range ! ";
1000 cerr<<j<<' '<<MAXTBKT<<endl;
1001 continue;
1002 }
1003 if (dig->fSignal >= THRESHOLD) bins[i*MAXZ+j].dig=dig;
1004 if (i==1 || i==npads || j==1 || j==MAXTBKT) dig->fSignal*=-1;
fe4da5cc 1005 }
3c0f9266 1006
fe4da5cc 1007 int ncl=0;
1008 int i,j;
8c555625 1009
3c0f9266 1010 for (i=1; i<MAXY-1; i++) {
1011 for (j=1; j<MAXZ-1; j++) {
1012 if (bins[i*MAXZ+j].dig == 0) continue;
1013 PreCluster c; c.summit=bins[i*MAXZ+j].dig; c.idx=ncls;
1014 FindPreCluster(i, j, MAXZ, bins, c);
fe4da5cc 1015 c.fY /= c.fQ;
1016 c.fZ /= c.fQ;
8c555625 1017
1018 double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
1019 c.fSigmaY2 = s2 + 1./12.;
3c0f9266 1020 c.fSigmaY2 *= par->GetPadPitchWidth()*par->GetPadPitchWidth();
1021 if (s2 != 0.) c.fSigmaY2 *= 0.17;
8c555625 1022
1023 s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
1024 c.fSigmaZ2 = s2 + 1./12.;
3c0f9266 1025 c.fSigmaZ2 *= par->GetZWidth()*par->GetZWidth();
1026 if (s2 != 0.) c.fSigmaZ2 *= 0.41;
1027
1028 c.fY = (c.fY - 0.5 - 0.5*npads)*par->GetPadPitchWidth();
1029 c.fZ = par->GetZWidth()*c.fZ;
1030 c.fZ -= 3.*par->GetZSigma(); // PASA delay
1031 c.fZ = sign*(z_end - c.fZ);
1032
1033 c.fSector=sec;
1034 c.fPadRow=row;
fe4da5cc 1035 c.fTracks[0]=c.summit->fTracks[0];
1036 c.fTracks[1]=c.summit->fTracks[1];
1037 c.fTracks[2]=c.summit->fTracks[2];
8c555625 1038
3c0f9266 1039 if (c.summit->fSignal<0) {
fe4da5cc 1040 c.fSigmaY2 *= 25.;
1041 c.fSigmaZ2 *= 4.;
1042 }
1043
1044 AddCluster(c); ncls++; ncl++;
1045 }
1046 }
3c0f9266 1047
fe4da5cc 1048 for (ndig=0; ndig<ndigits; ndig++) {
1049 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
3c0f9266 1050 int i=dig->fPad+1, j=dig->fTime+1;
1051 if (i > npads) {
1052 cerr<<"AliTPC::Digits2Clusters error: pad number is out of range ! ";
1053 cerr<<i<<' '<<npads<<endl;
1054 continue;
1055 }
1056 if (j > MAXTBKT) {
1057 cerr<<"AliTPC::Digits2Clusters error: time bucket is out of range ! ";
1058 cerr<<j<<' '<<MAXTBKT<<endl;
1059 continue;
1060 }
1061 if (TMath::Abs(dig->fSignal)>=par->GetZeroSup()) bins[i*MAXZ+j].dig=dig;
fe4da5cc 1062 }
1063
3c0f9266 1064 for (i=1; i<MAXY-1; i++) {
1065 for (j=1; j<MAXZ-1; j++) {
1066 if (bins[i*MAXZ+j].dig == 0) continue;
1067 PreCluster c; c.summit=bins[i*MAXZ+j].dig; c.idx=ncls;
1068 FindPreCluster(i, j, MAXZ, bins, c);
1069 if (c.ndigits < 2) continue; //noise cluster
8c555625 1070 if (c.npeaks>1) continue; //overlapped cluster
fe4da5cc 1071 c.fY /= c.fQ;
1072 c.fZ /= c.fQ;
8c555625 1073
1074 double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
1075 c.fSigmaY2 = s2 + 1./12.;
3c0f9266 1076 c.fSigmaY2 *= par->GetPadPitchWidth()*par->GetPadPitchWidth();
1077 if (s2 != 0.) c.fSigmaY2 *= 0.04;
8c555625 1078
1079 s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
1080 c.fSigmaZ2 = s2 + 1./12.;
3c0f9266 1081 c.fSigmaZ2 *= par->GetZWidth()*par->GetZWidth();
1082 if (s2 != 0.) c.fSigmaZ2 *= 0.10;
1083
1084 c.fY = (c.fY - 0.5 - 0.5*npads)*par->GetPadPitchWidth();
1085 c.fZ = par->GetZWidth()*c.fZ;
1086 c.fZ -= 3.*par->GetZSigma(); // PASA delay
1087 c.fZ = sign*(z_end - c.fZ);
1088
1089 c.fSector=sec;
1090 c.fPadRow=row;
fe4da5cc 1091 c.fTracks[0]=c.summit->fTracks[0];
1092 c.fTracks[1]=c.summit->fTracks[1];
1093 c.fTracks[2]=c.summit->fTracks[2];
1094
3c0f9266 1095 if (c.summit->fSignal<0) {
fe4da5cc 1096 c.fSigmaY2 *= 25.;
1097 c.fSigmaZ2 *= 4.;
1098 }
3c0f9266 1099
fe4da5cc 1100 if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;}
1101 else {
1102 new ((*fClusters)[c.idx]) AliTPCcluster(c);
1103 }
1104 }
1105 }
1106
1107 cerr<<"sector, row, digits, clusters: "
3c0f9266 1108 <<sec<<' '<<row<<' '<<ndigits<<' '<<ncl<<" \r";
fe4da5cc 1109
1110 fDigits->Clear();
3c0f9266 1111
1112 delete[] bins;
fe4da5cc 1113 }
1114}
1115
1116//_____________________________________________________________________________
1117void AliTPC::ElDiff(Float_t *xyz)
1118{
8c555625 1119 //--------------------------------------------------
fe4da5cc 1120 // calculates the diffusion of a single electron
8c555625 1121 //--------------------------------------------------
1122
1123 //-----------------------------------------------------------------
1124 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1125 //-----------------------------------------------------------------
1126 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
fe4da5cc 1127 Float_t driftl;
8c555625 1128
fe4da5cc 1129 Float_t z0=xyz[2];
8c555625 1130
fe4da5cc 1131 driftl=z_end-TMath::Abs(xyz[2]);
8c555625 1132
fe4da5cc 1133 if(driftl<0.01) driftl=0.01;
8c555625 1134
1135 // check the attachment
1136
fe4da5cc 1137 driftl=TMath::Sqrt(driftl);
8c555625 1138
1139 // Float_t sig_t = driftl*diff_t;
1140 //Float_t sig_l = driftl*diff_l;
1141 Float_t sig_t = driftl*fTPCParam->GetDiffT();
1142 Float_t sig_l = driftl*fTPCParam->GetDiffL();
fe4da5cc 1143 xyz[0]=gRandom->Gaus(xyz[0],sig_t);
1144 xyz[1]=gRandom->Gaus(xyz[1],sig_t);
1145 xyz[2]=gRandom->Gaus(xyz[2],sig_l);
8c555625 1146
fe4da5cc 1147 if (TMath::Abs(xyz[2])>z_end){
8c555625 1148 xyz[2]=TMath::Sign(z_end,z0);
fe4da5cc 1149 }
1150 if(xyz[2]*z0 < 0.){
8c555625 1151 Float_t eps = 0.0001;
1152 xyz[2]=TMath::Sign(eps,z0);
fe4da5cc 1153 }
1154}
1155
1156//_____________________________________________________________________________
1157void AliTPC::Hits2Clusters()
1158{
8c555625 1159 //--------------------------------------------------------
fe4da5cc 1160 // TPC simple cluster generator from hits
1161 // obtained from the TPC Fast Simulator
8c555625 1162 // The point errors are taken from the parametrization
1163 //--------------------------------------------------------
1164
1165 //-----------------------------------------------------------------
1166 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1167 //-----------------------------------------------------------------
3c0f9266 1168
8c555625 1169 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
fe4da5cc 1170 Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
1171 //
1578254f 1172 TParticle *particle; // pointer to a given particle
fe4da5cc 1173 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1174 TClonesArray *Particles; //pointer to the particle list
1175 Int_t sector,nhits;
1176 Int_t ipart;
1177 Float_t xyz[5];
1178 Float_t pl,pt,tanth,rpad,ratio;
fe4da5cc 1179 Float_t cph,sph;
1180
1181 //---------------------------------------------------------------
1182 // Get the access to the tracks
1183 //---------------------------------------------------------------
1184
1185 TTree *TH = gAlice->TreeH();
1186 Stat_t ntracks = TH->GetEntries();
3c0f9266 1187 Particles=gAlice->Particles();
fe4da5cc 1188
1189 //------------------------------------------------------------
1190 // Loop over all sectors (72 sectors)
3c0f9266 1191 // Sectors 0-35 are lower sectors, 0-17 z>0, 18-35 z<0
1192 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
fe4da5cc 1193 //
3c0f9266 1194 // First cluster for sector 0 starts at "0"
fe4da5cc 1195 //------------------------------------------------------------
1196
1197
fe4da5cc 1198 //
3c0f9266 1199 int Nsectors=fDigParam->GetParam().GetNSector();
1200 for(Int_t isec=0; isec<Nsectors; isec++){
8c555625 1201 //MI change
1202 fTPCParam->AdjustAngles(isec,cph,sph);
fe4da5cc 1203
1204 //------------------------------------------------------------
1205 // Loop over tracks
1206 //------------------------------------------------------------
1207
1208 for(Int_t track=0;track<ntracks;track++){
1209 ResetHits();
1210 TH->GetEvent(track);
1211 //
1212 // Get number of the TPC hits and a pointer
1213 // to the particles
1214 //
1215 nhits=fHits->GetEntriesFast();
fe4da5cc 1216 //
1217 // Loop over hits
1218 //
1219 for(Int_t hit=0;hit<nhits;hit++){
1220 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1221 sector=tpcHit->fSector; // sector number
1222 if(sector != isec) continue; //terminate iteration
1223 ipart=tpcHit->fTrack;
1578254f 1224 particle=(TParticle*)Particles->UncheckedAt(ipart);
1225 pl=particle->Pz();
1226 pt=particle->Pt();
fe4da5cc 1227 if(pt < 1.e-9) pt=1.e-9;
1228 tanth=pl/pt;
1229 tanth = TMath::Abs(tanth);
1230 rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
1231 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1232
1233 // space-point resolutions
1234
1235 sigma_rphi=SigmaY2(rpad,tanth,pt);
1236 sigma_z =SigmaZ2(rpad,tanth );
1237
1238 // cluster widths
1239
1240 cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
1241 cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
1242
1243 // temporary protection
1244
1245 if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
1246 if(sigma_z < 0.) sigma_z=0.4e-3;
1247 if(cl_rphi < 0.) cl_rphi=2.5e-3;
1248 if(cl_z < 0.) cl_z=2.5e-5;
1249
1250 //
3c0f9266 1251 // smearing --> rotate sectors firstly,
fe4da5cc 1252 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1253 //
1254 Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
1255 Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
1256 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi)); // y
3c0f9266 1257 Double_t alpha=(sector < fTPCParam->GetNInnerSector()) ?
1258 fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
8c555625 1259 if (TMath::Abs(xyz[0]/xprim) > TMath::Tan(0.5*alpha)) xyz[0]=yprim;
fe4da5cc 1260 xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z
8c555625 1261 if (TMath::Abs(xyz[1]) > 250) xyz[1]=tpcHit->fZ;
3c0f9266 1262 xyz[2]=tpcHit->fQ+1;// q; let it be not equal to zero (Y.Belikov)
fe4da5cc 1263 xyz[3]=sigma_rphi; // fSigmaY2
1264 xyz[4]=sigma_z; // fSigmaZ2
1283eee5 1265
fe4da5cc 1266 // and finally add the cluster
1283eee5 1267 Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, tpcHit->fPadRow};
fe4da5cc 1268 AddCluster(xyz,tracks);
1269
1270 } // end of loop over hits
1271 } // end of loop over tracks
1272
fe4da5cc 1273
1274 } // end of loop over sectors
1275
fe4da5cc 1276
1277}
1278
8c555625 1279
1280void AliTPC::Hits2Digits()
1281{
1282
1283 //----------------------------------------------------
1284 // Loop over all sectors (72 sectors)
3c0f9266 1285 // Sectors 0-35 are lower sectors, 0-17 z>0, 18-35 z<0
1286 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
8c555625 1287 //----
3c0f9266 1288 int Nsectors=fDigParam->GetParam().GetNSector();
1289 for(Int_t isec=0;isec<Nsectors;isec++) Hits2DigitsSector(isec);
8c555625 1290}
1291
1292
fe4da5cc 1293//_____________________________________________________________________________
8c555625 1294void AliTPC::Hits2DigitsSector(Int_t isec)
fe4da5cc 1295{
8c555625 1296 //-------------------------------------------------------------------
fe4da5cc 1297 // TPC conversion from hits to digits.
8c555625 1298 //-------------------------------------------------------------------
1299
1300 //-----------------------------------------------------------------
1301 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1302 //-----------------------------------------------------------------
1303
fe4da5cc 1304 //-------------------------------------------------------
8c555625 1305 // Get the access to the track hits
fe4da5cc 1306 //-------------------------------------------------------
8c555625 1307
1308 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1309 TTree *TH = gAlice->TreeH(); // pointer to the hits tree
fe4da5cc 1310 Stat_t ntracks = TH->GetEntries();
8c555625 1311
1312 if( ntracks > 0){
1313
1314 //-------------------------------------------
1315 // Only if there are any tracks...
1316 //-------------------------------------------
1317
fe4da5cc 1318
8c555625 1319 // TObjArrays for three neighouring pad-rows
1320
1321 TObjArray **rowTriplet = new TObjArray* [3];
fe4da5cc 1322
8c555625 1323 // TObjArray-s for each pad-row
1324
1325 TObjArray **row;
1326
1327 for (Int_t trip=0;trip<3;trip++){
1328 rowTriplet[trip]=new TObjArray;
fe4da5cc 1329 }
8c555625 1330
1331
fe4da5cc 1332
8c555625 1333 printf("*** Processing sector number %d ***\n",isec);
1334
1335 Int_t nrows =fTPCParam->GetNRow(isec);
1336
1337 row= new TObjArray* [nrows];
fe4da5cc 1338
8c555625 1339 MakeSector(isec,nrows,TH,ntracks,row);
1340
1341 //--------------------------------------------------------
1342 // Digitize this sector, row by row
1343 // row[i] is the pointer to the TObjArray of TVectors,
1344 // each one containing electrons accepted on this
1345 // row, assigned into tracks
1346 //--------------------------------------------------------
1347
1348 Int_t i;
1349
1350 for (i=0;i<nrows;i++){
1351
1352 // Triplets for i = 0 and i=1 are identical!
1353 // The same for i = nrows-1 and nrows!
1354
1355 if(i != 1 && i != nrows-1){
1356 MakeTriplet(i,rowTriplet,row);
1357 }
1358
1359 DigitizeRow(i,isec,rowTriplet);
1360
1361 fDigParam->Fill();
1362
1363 Int_t ndig=fDigParam->GetArray()->GetEntriesFast();
1364
1365 printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig);
1366
1367 ResetDigits(); // reset digits for this row after storing them
1368
1369 } // end of the sector digitization
1370
1371 // delete the last triplet
1372
1373 for (i=0;i<3;i++) rowTriplet[i]->Delete();
1374
1375 delete [] row; // delete the array of pointers to TObjArray-s
1376
1377 } // ntracks >0
1378} // end of Hits2Digits
1379//_____________________________________________________________________________
1380void AliTPC::MakeTriplet(Int_t row,
1381 TObjArray **rowTriplet, TObjArray **prow)
1382{
1383 //------------------------------------------------------------------
1384 // Makes the "triplet" of the neighbouring pad-row for the
1385 // digitization including the cross-talk between the pad-rows
1386 //------------------------------------------------------------------
1387
1388 //-----------------------------------------------------------------
1389 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1390 //-----------------------------------------------------------------
1391
1392 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1393 Float_t gasgain = fTPCParam->GetGasGain();
1394 Int_t nTracks[3];
1395
1396 Int_t nElements,nElectrons;
1397
1398 TVector *pv;
1399 TVector *tv;
1400
1401 //-------------------------------------------------------------------
1402 // pv is an "old" track, i.e. label + triplets of (x,y,z)
1403 // for each electron
1404 //
1405 //-------------------------------------------------------------------
1406
1407
1408 Int_t i1,i2;
1409 Int_t nel,nt;
1410
1411 if(row == 0 || row == 1){
1412
1413 // create entire triplet for the first AND the second row
1414
1415 nTracks[0] = prow[0]->GetEntries();
1416 nTracks[1] = prow[1]->GetEntries();
1417 nTracks[2] = prow[2]->GetEntries();
1418
1419 for(i1=0;i1<3;i1++){
1420 nt = nTracks[i1]; // number of tracks for this row
1421
1422 for(i2=0;i2<nt;i2++){
1423 pv = (TVector*)prow[i1]->At(i2);
1424 TVector &v1 = *pv;
1425 nElements = pv->GetNrows();
1426 nElectrons = (nElements-1)/3;
1427
1428 tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1429 TVector &v2 = *tv;
1430 v2(0)=v1(0); //track label
1431
1432 for(nel=0;nel<nElectrons;nel++){
1433 Int_t idx1 = nel*3;
1434 Int_t idx2 = nel*4;
1435 // Avalanche, including fluctuations
1436 Int_t aval = (Int_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1437 v2(idx2+1)= v1(idx1+1);
1438 v2(idx2+2)= v1(idx1+2);
1439 v2(idx2+3)= v1(idx1+3);
1440 v2(idx2+4)= (Float_t)aval; // in number of electrons!
1441 } // end of loop over electrons
1442 //
1443 // Add this track to a row
1444 //
1445
1446 rowTriplet[i1]->Add(tv);
1447
1448
1449 } // end of loop over tracks for this row
1450
1451 prow[i1]->Delete(); // remove "old" tracks
1452 delete prow[i1]; // delete TObjArray for this row
1453 prow[i1]=0; // set pointer to NULL
1454
1455 } // end of loop over row triplets
1456
1457
1458 }
1459 else{
1460
1461 rowTriplet[0]->Delete(); // remove old lower row
1462
1463 nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row
1464 nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row
1465 nTracks[2]=prow[row+1]->GetEntries(); // next row
1466
1467
1468 //-------------------------------------------
1469 // shift new tracks downwards
1470 //-------------------------------------------
1471
1472 for(i1=0;i1<nTracks[0];i1++){
1473 pv=(TVector*)rowTriplet[1]->At(i1);
1474 rowTriplet[0]->Add(pv);
1475 }
1476
1477 rowTriplet[1]->Clear(); // leave tracks on the heap!!!
1478
1479 for(i1=0;i1<nTracks[1];i1++){
1480 pv=(TVector*)rowTriplet[2]->At(i1);
1481 rowTriplet[1]->Add(pv);
1482 }
1483
1484 rowTriplet[2]->Clear(); // leave tracks on the heap!!!
1485
1486 //---------------------------------------------
1487 // Create new upper row
1488 //---------------------------------------------
1489
1490
1491
1492 for(i1=0;i1<nTracks[2];i1++){
1493 pv = (TVector*)prow[row+1]->At(i1);
1494 TVector &v1 = *pv;
1495 nElements = pv->GetNrows();
1496 nElectrons = (nElements-1)/3;
1497
1498 tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1499 TVector &v2 = *tv;
1500 v2(0)=v1(0); //track label
1501
1502 for(nel=0;nel<nElectrons;nel++){
1503
1504 Int_t idx1 = nel*3;
1505 Int_t idx2 = nel*4;
1506 // Avalanche, including fluctuations
1507 Int_t aval = (Int_t)
1508 (-gasgain*TMath::Log(gRandom->Rndm()));
1509
1510 v2(idx2+1)= v1(idx1+1);
1511 v2(idx2+2)= v1(idx1+2);
1512 v2(idx2+3)= v1(idx1+3);
1513 v2(idx2+4)= (Float_t)aval; // in number of electrons!
1514 } // end of loop over electrons
1515
1516 rowTriplet[2]->Add(tv);
1517
1518 } // end of loop over tracks
1519
1520 prow[row+1]->Delete(); // delete tracks for this row
1521 delete prow[row+1]; // delete TObjArray for this row
1522 prow[row+1]=0; // set a pointer to NULL
1523
1524 }
1525
1526} // end of MakeTriplet
1527//_____________________________________________________________________________
1528void AliTPC::ExB(Float_t *xyz)
1529{
1530 //------------------------------------------------
1531 // ExB at the wires and wire number calulation
1532 //------------------------------------------------
1533
1534 //-----------------------------------------------------------------
1535 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1536 //-----------------------------------------------------------------
1537 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1538
1539 Float_t x1=xyz[0];
1540 fTPCParam->GetWire(x1); //calculate nearest wire position
1541 Float_t dx=xyz[0]-x1;
1542 xyz[1]+=dx*fTPCParam->GetOmegaTau();
1543
1544} // end of ExB
1545//_____________________________________________________________________________
1546void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet)
1547{
1548 //-----------------------------------------------------------
1549 // Single row digitization, coupling from the neighbouring
1550 // rows taken into account
1551 //-----------------------------------------------------------
1552
1553 //-----------------------------------------------------------------
1554 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1555 //-----------------------------------------------------------------
1556
1557
1558 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1559 Float_t chipgain= fTPCParam->GetChipGain();
1560 Float_t zerosup = fTPCParam->GetZeroSup();
1561 Int_t nrows =fTPCParam->GetNRow(isec);
3c0f9266 1562 const int MAXTBKT=
1563 int((z_end+6*fTPCParam->GetZSigma())/fTPCParam->GetZWidth())+1;
8c555625 1564
1565 Int_t nTracks[3];
1566 Int_t n_of_pads[3];
1567 Int_t IndexRange[4];
1568 Int_t i1;
1569 Int_t iFlag;
1570
1571 //
1572 // iFlag = 0 -> inner row, iFlag = 1 -> the middle one, iFlag = 2 -> the outer one
1573 //
1574
1575 nTracks[0]=rowTriplet[0]->GetEntries(); // lower row
1576 nTracks[1]=rowTriplet[1]->GetEntries(); // middle row
1577 nTracks[2]=rowTriplet[2]->GetEntries(); // upper row
1578
1579
1580 if(irow == 0){
1581 iFlag=0;
1582 n_of_pads[0]=fTPCParam->GetNPads(isec,0);
1583 n_of_pads[1]=fTPCParam->GetNPads(isec,1);
1584 }
1585 else if(irow == nrows-1){
1586 iFlag=2;
1587 n_of_pads[1]=fTPCParam->GetNPads(isec,irow-1);
1588 n_of_pads[2]=fTPCParam->GetNPads(isec,irow);
1589 }
1590 else {
1591 iFlag=1;
1592 for(i1=0;i1<3;i1++){
1593 n_of_pads[i1]=fTPCParam->GetNPads(isec,irow-1+i1);
1594 }
1595 }
1596
1597 //
1598 // Integrated signal for this row
1599 // and a single track signal
1600 //
1601
3c0f9266 1602 TMatrix *m1 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTBKT-1); // integrated
1603 TMatrix *m2 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTBKT-1); // single
8c555625 1604
1605 //
1606
1607 TMatrix &Total = *m1;
1608
1609 // Array of pointers to the label-signal list
1610
3c0f9266 1611 Int_t NofDigits = n_of_pads[iFlag]*MAXTBKT; // number of digits for this row
8c555625 1612
1613 Float_t **pList = new Float_t* [NofDigits];
1614
1615 Int_t lp;
1616
1617 for(lp=0;lp<NofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1618
1619 //
1620 // Straight signal and cross-talk, cross-talk is integrated over all
1621 // tracks and added to the signal at the very end
1622 //
1623
1624
1625 for (i1=0;i1<nTracks[iFlag];i1++){
1626
1627 m2->Zero(); // clear single track signal matrix
1628
1629 Float_t TrackLabel =
1630 GetSignal(rowTriplet[iFlag],i1,n_of_pads[iFlag],m2,m1,IndexRange);
1631
1632 GetList(TrackLabel,n_of_pads[iFlag],m2,IndexRange,pList);
1633
1634 }
1635
1636 //
1637 // Cross talk from the neighbouring pad-rows
1638 //
1639
3c0f9266 1640 TMatrix *m3 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTBKT-1); // cross-talk
8c555625 1641
1642 TMatrix &Cross = *m3;
1643
1644 if(iFlag == 0){
1645
1646 // cross-talk from the outer row only (first pad row)
1647
1648 GetCrossTalk(0,rowTriplet[1],nTracks[1],n_of_pads,m3);
1649
1650 }
1651 else if(iFlag == 2){
1652
1653 // cross-talk from the inner row only (last pad row)
1654
1655 GetCrossTalk(2,rowTriplet[1],nTracks[1],n_of_pads,m3);
1656
1657 }
1658 else{
1659
1660 // cross-talk from both inner and outer rows
1661
1662 GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner
1663 GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer
1664 }
1665
1666 Total += Cross; // add the cross-talk
1667
1668 //
1669 // Convert analog signal to ADC counts
1670 //
1671
1672 Int_t tracks[3];
1673 Int_t digits[5];
1674
1675
1283eee5 1676 for(Int_t ip=0;ip<n_of_pads[iFlag];ip++){
3c0f9266 1677 for(Int_t it=0;it<MAXTBKT;it++){
8c555625 1678
1679 Float_t q = Total(ip,it);
1680
1283eee5 1681 Int_t gi =it*n_of_pads[iFlag]+ip; // global index
8c555625 1682
1683 q = gRandom->Gaus(q,fTPCParam->GetNoise()); // apply noise
1684 q *= (q_el*1.e15); // convert to fC
1685 q *= chipgain; // convert to mV
1686 q *= (adc_sat/dyn_range); // convert to ADC counts
1687
1688 if(q <zerosup) continue; // do not fill zeros
1689 if(q > adc_sat) q = adc_sat; // saturation
1690
1691 //
1692 // "real" signal or electronic noise (list = -1)?
1693 //
1694
1695 for(Int_t j1=0;j1<3;j1++){
1696 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1697 }
1698
1699 digits[0]=isec;
1283eee5 1700 digits[1]=irow;
8c555625 1701 digits[2]=ip;
1702 digits[3]=it;
1703 digits[4]= (Int_t)q;
1704
1705 // Add this digit
1706
1707 AddDigit(tracks,digits);
1708
1709 } // end of loop over time buckets
1710 } // end of lop over pads
1711
1712 //
1713 // This row has been digitized, delete nonused stuff
1714 //
1715
1716 for(lp=0;lp<NofDigits;lp++){
1717 if(pList[lp]) delete [] pList[lp];
1718 }
1719
1720 delete [] pList;
1721
1722 delete m1;
1723 delete m2;
1724 delete m3;
1725
1726} // end of DigitizeRow
1727//_____________________________________________________________________________
1728Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, Int_t np, TMatrix *m1, TMatrix *m2,
1729 Int_t *IndexRange)
1730{
1731
1732 //---------------------------------------------------------------
1733 // Calculates 2-D signal (pad,time) for a single track,
1734 // returns a pointer to the signal matrix and the track label
1735 // No digitization is performed at this level!!!
1736 //---------------------------------------------------------------
1737
1738 //-----------------------------------------------------------------
1739 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1740 //-----------------------------------------------------------------
1741
1742 TVector *tv;
1743 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1744 AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
1745 AliTPCRF1D * fRF = &(fDigParam->GetRF());
3c0f9266 1746 const int MAXTBKT=
1747 int((z_end+6*fTPCParam->GetZSigma())/fTPCParam->GetZWidth())+1;
8c555625 1748
1749 //to make the code faster we put parameters to the stack
1750
1751 Float_t zwidth = fTPCParam->GetZWidth();
7ee081c2 1752 Float_t zwidthm1 =1./zwidth;
8c555625 1753
1754 tv = (TVector*)p1->At(ntr); // pointer to a track
1755 TVector &v = *tv;
1756
1757 Float_t label = v(0);
1758
1283eee5 1759 Int_t CentralPad = (np-1)/2;
8c555625 1760 Int_t PadNumber;
1761 Int_t nElectrons = (tv->GetNrows()-1)/4;
1762 Float_t range=((np-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); // pad range
8c555625 1763
1283eee5 1764 range -= 0.5; // dead zone, 5mm from the edge, according to H.G. Fischer
1765
8c555625 1766 Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1767
1768
1769 Float_t PadSignal[7]; // signal from a single electron
1770
1771 TMatrix &signal = *m1;
1772 TMatrix &total = *m2;
1773
1774
1775 IndexRange[0]=9999; // min pad
1776 IndexRange[1]=-1; // max pad
1777 IndexRange[2]=9999; //min time
1778 IndexRange[3]=-1; // max time
1779
1780 //
1781 // Loop over all electrons
1782 //
1783
1784 for(Int_t nel=0; nel<nElectrons; nel++){
1785 Int_t idx=nel*4;
1786 Float_t xwire = v(idx+1);
1787 Float_t y = v(idx+2);
1788 Float_t z = v(idx+3);
1789
1790
1791 Float_t absy=TMath::Abs(y);
1792
1793 if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
1794 PadNumber=CentralPad;
1795 }
1796 else if (absy < range){
8022212f 1797 PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth()+1.);
8c555625 1798 PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
1799 }
1800 else continue; // electron out of pad-range , lost at the sector edge
1801
1802 Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
1803
1804
1805 Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
a03540ae 1806 for (Int_t i=0;i<7;i++){
b6738b36 1807 PadSignal[i]=fPRF2D->GetPRF(-dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
a03540ae 1808 PadSignal[i] *= fTPCParam->GetPadCoupling();
1809 }
8c555625 1810
1283eee5 1811 Int_t LeftPad = TMath::Max(0,PadNumber-3);
1812 Int_t RightPad = TMath::Min(np-1,PadNumber+3);
8c555625 1813
1814 Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
1815 Int_t pmax=RightPad-PadNumber+3; // upper index
1816
1283eee5 1817 Float_t z_drift = z*zwidthm1;
8c555625 1818 Float_t z_offset = z_drift-(Int_t)z_drift;
3c0f9266 1819
1820 Int_t FirstBucket = (Int_t)z_drift; // numbering starts from "0"
8c555625 1821
1822
1823 // loop over time bins (4 bins is enough - 3 sigma truncated Gaussian)
3c0f9266 1824
8c555625 1825 for (Int_t i2=0;i2<4;i2++){
1826 Int_t TrueTime = FirstBucket+i2; // current time bucket
3c0f9266 1827 Float_t dz = (Float_t(i2)+1.-z_offset)*zwidth;
8c555625 1828 Float_t ampl = fRF->GetRF(dz);
3c0f9266 1829 if( (TrueTime>MAXTBKT-1) ) break; // beyond the time range
8c555625 1830
1831 IndexRange[2]=TMath::Min(IndexRange[2],TrueTime); // min time
1832 IndexRange[3]=TMath::Max(IndexRange[3],TrueTime); // max time
1833
1834 // loop over pads, from pmin to pmax
1835 for(Int_t i3=pmin;i3<=pmax;i3++){
1836 Int_t TruePad = LeftPad+i3-pmin;
1837 IndexRange[0]=TMath::Min(IndexRange[0],TruePad); // min pad
1838 IndexRange[1]=TMath::Max(IndexRange[1],TruePad); // max pad
1839 signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!!
1840 total(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!!
1841 } // end of pads loop
1842 } // end of time loop
1843 } // end of loop over electrons
1844
1845 return label; // returns track label when finished
1846}
1847
1848//_____________________________________________________________________________
1849void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange,
1850 Float_t **pList)
1851{
1852 //----------------------------------------------------------------------
1853 // Updates the list of tracks contributing to digits for a given row
1854 //----------------------------------------------------------------------
1855
1856 //-----------------------------------------------------------------
1857 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1858 //-----------------------------------------------------------------
1859
1860 TMatrix &signal = *m;
1861
1862 // lop over nonzero digits
1863
1864 for(Int_t it=IndexRange[2];it<IndexRange[3]+1;it++){
1865 for(Int_t ip=IndexRange[0];ip<IndexRange[1]+1;ip++){
1866
1867
1283eee5 1868 Int_t GlobalIndex = it*np+ip; // GlobalIndex starts from 0!
8c555625 1869
1870 if(!pList[GlobalIndex]){
1871
1872 //
1873 // Create new list (6 elements - 3 signals and 3 labels),
1874 // but only if the signal is > 0.
1875 //
1876
1877 if(signal(ip,it)>0.){
1878
1879 pList[GlobalIndex] = new Float_t [6];
1880
1881 // set list to -1
1882
1883 *pList[GlobalIndex] = -1.;
1884 *(pList[GlobalIndex]+1) = -1.;
1885 *(pList[GlobalIndex]+2) = -1.;
1886 *(pList[GlobalIndex]+3) = -1.;
1887 *(pList[GlobalIndex]+4) = -1.;
1888 *(pList[GlobalIndex]+5) = -1.;
1889
1890
1891 *pList[GlobalIndex] = label;
1892 *(pList[GlobalIndex]+3) = signal(ip,it);}
1893 }
1894 else{
1895
1896 // check the signal magnitude
1897
1898 Float_t highest = *(pList[GlobalIndex]+3);
1899 Float_t middle = *(pList[GlobalIndex]+4);
1900 Float_t lowest = *(pList[GlobalIndex]+5);
1901
1902 //
1903 // compare the new signal with already existing list
1904 //
1905
1906 if(signal(ip,it)<lowest) continue; // neglect this track
1907
1908 //
1909
1910 if (signal(ip,it)>highest){
1911 *(pList[GlobalIndex]+5) = middle;
1912 *(pList[GlobalIndex]+4) = highest;
1913 *(pList[GlobalIndex]+3) = signal(ip,it);
1914
1915 *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1916 *(pList[GlobalIndex]+1) = *pList[GlobalIndex];
1917 *pList[GlobalIndex] = label;
1918 }
1919 else if (signal(ip,it)>middle){
1920 *(pList[GlobalIndex]+5) = middle;
1921 *(pList[GlobalIndex]+4) = signal(ip,it);
1922
1923 *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1924 *(pList[GlobalIndex]+1) = label;
1925 }
1926 else{
1927 *(pList[GlobalIndex]+5) = signal(ip,it);
1928 *(pList[GlobalIndex]+2) = label;
1929 }
1930 }
1931
1932 } // end of loop over pads
1933 } // end of loop over time bins
1934
1935
1936
1937
1938}//end of GetList
1939//___________________________________________________________________
1940void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1941 Stat_t ntracks,TObjArray **row)
1942{
1943
1944 //-----------------------------------------------------------------
1945 // Prepares the sector digitization, creates the vectors of
1946 // tracks for each row of this sector. The track vector
1947 // contains the track label and the position of electrons.
1948 //-----------------------------------------------------------------
1949
1950 //-----------------------------------------------------------------
1951 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1952 //-----------------------------------------------------------------
1953
1954 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1955 Int_t i;
1956 Float_t xyz[3];
1957
1958 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1959
1960 //----------------------------------------------
1961 // Create TObjArray-s, one for each row,
1962 // each TObjArray will store the TVectors
1963 // of electrons, one TVector per each track.
1964 //----------------------------------------------
1965
1966 for(i=0; i<nrows; i++){
1967 row[i] = new TObjArray;
1968 }
1969 Int_t *n_of_electrons = new Int_t [nrows]; // electron counter for each row
1970 TVector **tr = new TVector* [nrows]; //pointers to the track vectors
1971
1972 //--------------------------------------------------------------------
1973 // Loop over tracks, the "track" contains the full history
1974 //--------------------------------------------------------------------
1975
1976 Int_t previousTrack,currentTrack;
1977 previousTrack = -1; // nothing to store so far!
1978
1979 for(Int_t track=0;track<ntracks;track++){
1980
1981 ResetHits();
1982
1983 TH->GetEvent(track); // get next track
1984 Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1985
1986 if(nhits == 0) continue; // no hits in the TPC for this track
1987
1988 //--------------------------------------------------------------
1989 // Loop over hits
1990 //--------------------------------------------------------------
1991
1992 for(Int_t hit=0;hit<nhits;hit++){
1993
1994 tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1995
1996 Int_t sector=tpcHit->fSector; // sector number
1997 if(sector != isec) continue; //terminate iteration
1998
1999 currentTrack = tpcHit->fTrack; // track number
2000 if(currentTrack != previousTrack){
2001
2002 // store already filled fTrack
2003
2004 for(i=0;i<nrows;i++){
2005 if(previousTrack != -1){
2006 if(n_of_electrons[i]>0){
2007 TVector &v = *tr[i];
2008 v(0) = previousTrack;
2009 tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary
2010 row[i]->Add(tr[i]);
2011 }
2012 else{
2013 delete tr[i]; // delete empty TVector
2014 tr[i]=0;
2015 }
2016 }
2017
2018 n_of_electrons[i]=0;
2019 tr[i] = new TVector(361); // TVectors for the next fTrack
2020
2021 } // end of loop over rows
2022
2023 previousTrack=currentTrack; // update track label
2024 }
2025
2026 Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2027
2028 //---------------------------------------------------
2029 // Calculate the electron attachment probability
2030 //---------------------------------------------------
2031
2032 Float_t time = 1.e6*(z_end-TMath::Abs(tpcHit->fZ))/fTPCParam->GetDriftV();
2033 // in microseconds!
2034 Float_t AttProb = fTPCParam->GetAttCoef()*
2035 fTPCParam->GetOxyCont()*time; // fraction!
2036
2037 //-----------------------------------------------
2038 // Loop over electrons
2039 //-----------------------------------------------
2040
2041 for(Int_t nel=0;nel<QI;nel++){
2042 // skip if electron lost due to the attachment
2043 if((gRandom->Rndm(0)) < AttProb) continue; // electron lost!
2044 xyz[0]=tpcHit->fX;
2045 xyz[1]=tpcHit->fY;
2046 xyz[2]=tpcHit->fZ;
2047 ElDiff(xyz); // Appply the diffusion
2048 Int_t row_number;
2049 fTPCParam->XYZtoCRXYZ(xyz,isec,row_number,3);
2050
2051 //transform position to local coordinates
2052 //option 3 means that we calculate x position relative to
2053 //nearest pad row
2054
2055 if ((row_number<0)||row_number>=fTPCParam->GetNRow(isec)) continue;
2056 ExB(xyz); // ExB effect at the sense wires
2057 n_of_electrons[row_number]++;
2058 //----------------------------------
2059 // Expand vector if necessary
2060 //----------------------------------
fe4da5cc 2061 if(n_of_electrons[row_number]>120){
2062 Int_t range = tr[row_number]->GetNrows();
8c555625 2063 if(n_of_electrons[row_number] > (range-1)/3){
2064 tr[row_number]->ResizeTo(range+150); // Add 50 electrons
fe4da5cc 2065 }
2066 }
2067
fe4da5cc 2068 TVector &v = *tr[row_number];
8c555625 2069 Int_t idx = 3*n_of_electrons[row_number]-2;
2070
2071 v(idx)= xyz[0]; // X
2072 v(idx+1)=xyz[1]; // Y (along the pad-row)
2073 v(idx+2)=xyz[2]; // Z
2074
2075 } // end of loop over electrons
2076
2077 } // end of loop over hits
2078 } // end of loop over tracks
2079
2080 //
2081 // store remaining track (the last one) if not empty
2082 //
2083
2084 for(i=0;i<nrows;i++){
2085 if(n_of_electrons[i]>0){
2086 TVector &v = *tr[i];
2087 v(0) = previousTrack;
2088 tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary
2089 row[i]->Add(tr[i]);
fe4da5cc 2090 }
2091 else{
8c555625 2092 delete tr[i];
2093 tr[i]=0;
2094 }
2095 }
2096
2097 delete [] tr;
2098 delete [] n_of_electrons;
2099
2100} // end of MakeSector
2101//_____________________________________________________________________________
2102void AliTPC::GetCrossTalk (Int_t iFlag,TObjArray *p,Int_t ntracks,Int_t *npads,
2103 TMatrix *m)
2104{
2105
2106 //-------------------------------------------------------------
2107 // Calculates the cross-talk from one row (inner or outer)
2108 //-------------------------------------------------------------
2109
2110 //-----------------------------------------------------------------
2111 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2112 //-----------------------------------------------------------------
2113
2114 //
2115 // iFlag=2 & 3 --> cross-talk from the inner row
2116 // iFlag=0 & 4 --> cross-talk from the outer row
2117 //
2118 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
2119 AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
2120 AliTPCRF1D * fRF = &(fDigParam->GetRF());
3c0f9266 2121 const int MAXTBKT=
2122 int((z_end+6*fTPCParam->GetZSigma())/fTPCParam->GetZWidth())+1;
8c555625 2123
2124 //to make code faster
2125
2126 Float_t zwidth = fTPCParam->GetZWidth();
2127 Float_t zwidthm1 =1/fTPCParam->GetZWidth();
2128
2129 Int_t PadNumber;
2130 Float_t xwire;
2131
2132 Int_t nPadsSignal; // for this pads the signal is calculated
2133 Float_t range; // sense wire range
2134 Int_t nPadsDiff;
2135
2136 Float_t IneffFactor=0.5; // gain inefficiency close to the edges
2137
2138 if(iFlag == 0){
2139 // 1-->0
2140 nPadsSignal = npads[1];
2141 range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2142 nPadsDiff = (npads[1]-npads[0])/2;
2143 }
2144 else if (iFlag == 2){
2145 // 1-->2
2146 nPadsSignal = npads[2];
2147 range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2148 nPadsDiff = 0;
2149 }
2150 else if (iFlag == 3){
2151 // 0-->1
2152 nPadsSignal = npads[1];
2153 range = ((npads[0]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2154 nPadsDiff = 0;
2155 }
2156 else{
2157 // 2-->1
2158 nPadsSignal = npads[2];
2159 range = ((npads[2]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2160 nPadsDiff = (npads[2]-npads[1])/2;
2161 }
2162
2163 range-=0.5; // dead zone close to the edges
2164
2165 TVector *tv;
2166 TMatrix &signal = *m;
2167
1283eee5 2168 Int_t CentralPad = (nPadsSignal-1)/2;
8c555625 2169 Float_t PadSignal[7]; // signal from a single electron
2170 // Loop over tracks
2171 for(Int_t nt=0;nt<ntracks;nt++){
2172 tv=(TVector*)p->At(nt); // pointer to a track
2173 TVector &v = *tv;
2174 Int_t nElectrons = (tv->GetNrows()-1)/4;
2175 // Loop over electrons
2176 for(Int_t nel=0; nel<nElectrons; nel++){
2177 Int_t idx=nel*4;
2178 xwire=v(idx+1);
2179
2180 if (iFlag==0) xwire+=fTPCParam->GetPadPitchLength();
2181 if (iFlag==2) xwire-=fTPCParam->GetPadPitchLength();
2182 if (iFlag==3) xwire-=fTPCParam->GetPadPitchLength();
2183 if (iFlag==4) xwire+=fTPCParam->GetPadPitchLength();
2184
2185 // electron acceptance for the cross-talk (only the closest wire)
2186
2187 Float_t dxMax = fTPCParam->GetPadPitchLength()*0.5+fTPCParam->GetWWPitch();
2188 if(TMath::Abs(xwire)>dxMax) continue;
2189
2190 Float_t y = v(idx+2);
2191 Float_t z = v(idx+3);
2192 Float_t absy=TMath::Abs(y);
2193
2194 if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
2195 PadNumber=CentralPad;
2196 }
2197 else if (absy < range){
8022212f 2198 PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.);
8c555625 2199 PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
2200 }
2201 else continue; // electron out of sense wire range, lost at the sector edge
2202
2203 Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
2204
2205 Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
2206
a03540ae 2207 for (Int_t i=0;i<7;i++){
b6738b36 2208 PadSignal[i]=fPRF2D->GetPRF(-dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
8c555625 2209
a03540ae 2210 PadSignal[i] *= fTPCParam->GetPadCoupling();
2211 }
8c555625 2212 // real pad range
2213
1283eee5 2214 Int_t LeftPad = TMath::Max(0,PadNumber-3);
2215 Int_t RightPad = TMath::Min(nPadsSignal-1,PadNumber+3);
8c555625 2216
2217 Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
2218 Int_t pmax=RightPad-PadNumber+3; // upper index
2219
2220
1283eee5 2221 Float_t z_drift = z*zwidthm1;
8c555625 2222 Float_t z_offset = z_drift-(Int_t)z_drift;
3c0f9266 2223
2224 Int_t FirstBucket = (Int_t)z_drift; // numbering starts from "0"
2225
8c555625 2226 for (Int_t i2=0;i2<4;i2++){
2227 Int_t TrueTime = FirstBucket+i2; // current time bucket
3c0f9266 2228 Float_t dz = (Float_t(i2)+1.- z_offset)*zwidth;
8c555625 2229 Float_t ampl = fRF->GetRF(dz);
3c0f9266 2230 if((TrueTime>MAXTBKT-1)) break; // beyond the time range
8c555625 2231
2232
2233 // loop over pads, from pmin to pmax
2234
2235 for(Int_t i3=pmin;i3<pmax+1;i3++){
2236 Int_t TruePad = LeftPad+i3-pmin;
2237
1283eee5 2238 if(TruePad<nPadsDiff || TruePad > nPadsSignal-nPadsDiff-1) continue;
8c555625 2239
2240 TruePad -= nPadsDiff;
2241 signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!
2242
2243 } // end of loop over pads
2244 } // end of loop over time bins
2245
2246 } // end of loop over electrons
2247
2248 } // end of loop over tracks
2249
2250} // end of GetCrossTalk
2251
2252
fe4da5cc 2253
2254//_____________________________________________________________________________
2255void AliTPC::Init()
2256{
2257 //
2258 // Initialise TPC detector after definition of geometry
2259 //
2260 Int_t i;
2261 //
2262 printf("\n");
2263 for(i=0;i<35;i++) printf("*");
2264 printf(" TPC_INIT ");
2265 for(i=0;i<35;i++) printf("*");
2266 printf("\n");
2267 //
2268 for(i=0;i<80;i++) printf("*");
2269 printf("\n");
2270}
2271
2272//_____________________________________________________________________________
2273void AliTPC::MakeBranch(Option_t* option)
2274{
2275 //
2276 // Create Tree branches for the TPC.
2277 //
2278 Int_t buffersize = 4000;
2279 char branchname[10];
2280 sprintf(branchname,"%s",GetName());
2281
2282 AliDetector::MakeBranch(option);
2283
2284 char *D = strstr(option,"D");
2285
2286 if (fDigits && gAlice->TreeD() && D) {
2287 gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
2288 printf("Making Branch %s for digits\n",branchname);
2289 }
2290
2291 char *R = strstr(option,"R");
2292
2293 if (fClusters && gAlice->TreeR() && R) {
2294 gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
2295 printf("Making Branch %s for Clusters\n",branchname);
2296 }
2297}
2298
2299//_____________________________________________________________________________
2300void AliTPC::ResetDigits()
2301{
2302 //
2303 // Reset number of digits and the digits array for this detector
2304 // reset clusters
2305 //
2306 fNdigits = 0;
8c555625 2307 // if (fDigits) fDigits->Clear();
2308 if (fDigParam->GetArray()!=0) fDigParam->GetArray()->Clear();
fe4da5cc 2309 fNclusters = 0;
2310 if (fClusters) fClusters->Clear();
2311}
2312
2313//_____________________________________________________________________________
2314void AliTPC::SetSecAL(Int_t sec)
2315{
8c555625 2316 //---------------------------------------------------
fe4da5cc 2317 // Activate/deactivate selection for lower sectors
8c555625 2318 //---------------------------------------------------
2319
2320 //-----------------------------------------------------------------
2321 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2322 //-----------------------------------------------------------------
2323
fe4da5cc 2324 fSecAL = sec;
2325}
2326
2327//_____________________________________________________________________________
2328void AliTPC::SetSecAU(Int_t sec)
2329{
8c555625 2330 //----------------------------------------------------
fe4da5cc 2331 // Activate/deactivate selection for upper sectors
8c555625 2332 //---------------------------------------------------
2333
2334 //-----------------------------------------------------------------
2335 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2336 //-----------------------------------------------------------------
2337
fe4da5cc 2338 fSecAU = sec;
2339}
2340
2341//_____________________________________________________________________________
2342void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2343{
8c555625 2344 //----------------------------------------
fe4da5cc 2345 // Select active lower sectors
8c555625 2346 //----------------------------------------
2347
2348 //-----------------------------------------------------------------
2349 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2350 //-----------------------------------------------------------------
2351
fe4da5cc 2352 fSecLows[0] = s1;
2353 fSecLows[1] = s2;
2354 fSecLows[2] = s3;
2355 fSecLows[3] = s4;
2356 fSecLows[4] = s5;
2357 fSecLows[5] = s6;
2358}
2359
2360//_____________________________________________________________________________
2361void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2362 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2363 Int_t s11 , Int_t s12)
2364{
8c555625 2365 //--------------------------------
fe4da5cc 2366 // Select active upper sectors
8c555625 2367 //--------------------------------
2368
2369 //-----------------------------------------------------------------
2370 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2371 //-----------------------------------------------------------------
2372
fe4da5cc 2373 fSecUps[0] = s1;
2374 fSecUps[1] = s2;
2375 fSecUps[2] = s3;
2376 fSecUps[3] = s4;
2377 fSecUps[4] = s5;
2378 fSecUps[5] = s6;
2379 fSecUps[6] = s7;
2380 fSecUps[7] = s8;
2381 fSecUps[8] = s9;
2382 fSecUps[9] = s10;
2383 fSecUps[10] = s11;
2384 fSecUps[11] = s12;
2385}
2386
2387//_____________________________________________________________________________
2388void AliTPC::SetSens(Int_t sens)
2389{
8c555625 2390
2391 //-------------------------------------------------------------
2392 // Activates/deactivates the sensitive strips at the center of
2393 // the pad row -- this is for the space-point resolution calculations
2394 //-------------------------------------------------------------
2395
2396 //-----------------------------------------------------------------
2397 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2398 //-----------------------------------------------------------------
2399
fe4da5cc 2400 fSens = sens;
2401}
4b0fdcad 2402
2403void AliTPC::SetSide(Float_t side)
2404{
2405 fSide = side;
2406
2407}
1283eee5 2408//____________________________________________________________________________
2409void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2410 Float_t p2,Float_t p3)
2411{
fe4da5cc 2412
1283eee5 2413 fNoComp = nc;
2414
2415 fMixtComp[0]=c1;
2416 fMixtComp[1]=c2;
2417 fMixtComp[2]=c3;
2418
2419 fMixtProp[0]=p1;
2420 fMixtProp[1]=p2;
2421 fMixtProp[2]=p3;
2422
2423
2424}
fe4da5cc 2425//_____________________________________________________________________________
2426void AliTPC::Streamer(TBuffer &R__b)
2427{
2428 //
2429 // Stream an object of class AliTPC.
2430 //
2431 if (R__b.IsReading()) {
2432 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2433 AliDetector::Streamer(R__b);
2434 if (R__v < 2) return;
2435 R__b >> fNsectors;
2436 R__b >> fNclusters;
2437 R__b >> fNtracks;
fe4da5cc 2438 fDigitsIndex = new Int_t[fNsectors+1];
2439 } else {
2440 R__b.WriteVersion(AliTPC::IsA());
2441 AliDetector::Streamer(R__b);
2442 R__b << fNsectors;
2443 R__b << fNclusters;
2444 R__b << fNtracks;
2445 }
2446}
2447
2448ClassImp(AliTPCcluster)
2449
2450//_____________________________________________________________________________
2451AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
2452{
2453 //
2454 // Creates a simulated cluster for the TPC
2455 //
2456 fTracks[0] = lab[0];
2457 fTracks[1] = lab[1];
2458 fTracks[2] = lab[2];
2459 fSector = lab[3];
2460 fPadRow = lab[4];
2461 fY = hits[0];
2462 fZ = hits[1];
2463 fQ = hits[2];
2464 fSigmaY2 = hits[3];
2465 fSigmaZ2 = hits[4];
2466}
2467
2468//_____________________________________________________________________________
8c555625 2469void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const
fe4da5cc 2470{
2471 //
8c555625 2472 // Transformation from local to global coordinate system
fe4da5cc 2473 //
3c0f9266 2474 x[0]=par->GetPadRowRadii(fSector,fPadRow);
8c555625 2475 x[1]=fY;
2476 x[2]=fZ;
3c0f9266 2477 par->CRXYZtoXYZ(x,fSector,fPadRow,1);
2478 x[2]=fZ;
fe4da5cc 2479}
2480
8c555625 2481//_____________________________________________________________________________
2482Int_t AliTPCcluster::Compare(TObject * o)
2483{
2484 //
2485 // compare two clusters according y coordinata
3c0f9266 2486 //
8c555625 2487 AliTPCcluster *cl= (AliTPCcluster *)o;
2488 if (fY<cl->fY) return -1;
2489 if (fY==cl->fY) return 0;
2490 return 1;
2491}
2492
2493Bool_t AliTPCcluster::IsSortable() const
2494{
2495 //
2496 //make AliTPCcluster sortabale
3c0f9266 2497 //
8c555625 2498 return kTRUE;
2499}
2500
2501
2502
fe4da5cc 2503ClassImp(AliTPCdigit)
2504
2505//_____________________________________________________________________________
2506AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2507 AliDigit(tracks)
2508{
2509 //
2510 // Creates a TPC digit object
2511 //
2512 fSector = digits[0];
2513 fPadRow = digits[1];
2514 fPad = digits[2];
2515 fTime = digits[3];
2516 fSignal = digits[4];
2517}
2518
2519
2520ClassImp(AliTPChit)
2521
2522//_____________________________________________________________________________
2523AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2524AliHit(shunt,track)
2525{
2526 //
2527 // Creates a TPC hit object
2528 //
2529 fSector = vol[0];
2530 fPadRow = vol[1];
2531 fX = hits[0];
2532 fY = hits[1];
2533 fZ = hits[2];
2534 fQ = hits[3];
2535}
2536
2537
2538ClassImp(AliTPCtrack)
2539
2540//_____________________________________________________________________________
2541AliTPCtrack::AliTPCtrack(Float_t *hits)
2542{
2543 //
2544 // Default creator for a TPC reconstructed track object
2545 //
3c0f9266 2546 fX=hits[0]; // This is dummy code !
fe4da5cc 2547}
2548
3c0f9266 2549AliTPCtrack::AliTPCtrack(const AliTPCcluster *c,const TVector& xx,
2550 const TMatrix& CC, Double_t xref, Double_t alpha):
2551 x(xx),C(CC),fClusters(200)
fe4da5cc 2552{
3c0f9266 2553 //-----------------------------------------------------------------
2554 // This is the main track constructor.
fe4da5cc 2555 //
3c0f9266 2556 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2557 //-----------------------------------------------------------------
2558 fX=xref;
2559 fAlpha=alpha;
2560 fChi2=0.;
2561 fClusters.AddLast((AliTPCcluster*)(c));
fe4da5cc 2562}
2563
2564//_____________________________________________________________________________
2565AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
3c0f9266 2566 fClusters(t.fClusters.GetEntriesFast())
fe4da5cc 2567{
3c0f9266 2568 //-----------------------------------------------------------------
2569 // This is a track copy constructor.
fe4da5cc 2570 //
3c0f9266 2571 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2572 //-----------------------------------------------------------------
2573 fX=t.fX;
2574 fChi2=t.fChi2;
fe4da5cc 2575 fAlpha=t.fAlpha;
3c0f9266 2576 int n=t.fClusters.GetEntriesFast();
2577 for (int i=0; i<n; i++) fClusters.AddLast(t.fClusters.UncheckedAt(i));
fe4da5cc 2578}
2579
2580//_____________________________________________________________________________
3c0f9266 2581Int_t AliTPCtrack::Compare(TObject *o) {
2582 //-----------------------------------------------------------------
2583 // This function compares tracks according to their curvature.
fe4da5cc 2584 //
3c0f9266 2585 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2586 //-----------------------------------------------------------------
2587 AliTPCtrack *t=(AliTPCtrack*)o;
2588 Double_t co=t->GetSigmaY2();
2589 Double_t c =GetSigmaY2();
2590 if (c>co) return 1;
2591 else if (c<co) return -1;
2592 return 0;
fe4da5cc 2593}
2594
2595//_____________________________________________________________________________
2596int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
2597{
3c0f9266 2598 //-----------------------------------------------------------------
2599 // This function propagates a track to a reference plane x=xk.
fe4da5cc 2600 //
3c0f9266 2601 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2602 //-----------------------------------------------------------------
fe4da5cc 2603 if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
3c0f9266 2604 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
fe4da5cc 2605 return 0;
2606 }
8c555625 2607
3c0f9266 2608 Double_t x1=fX, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
fe4da5cc 2609 Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
2610 Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
2611
2612 x(0) += dx*(c1+c2)/(r1+r2);
2613 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
8c555625 2614
fe4da5cc 2615 TMatrix F(5,5); F.UnitMatrix();
2616 Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
2617 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2618 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2619 Double_t cr=c1*r2+c2*r1;
2620 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2621 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2622 F(1,4)= dx*cc/cr;
2623 TMatrix tmp(F,TMatrix::kMult,C);
2624 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2625
3c0f9266 2626 fX=x2;
fe4da5cc 2627
2628 //Multiple scattering******************
3c0f9266 2629 Double_t ey=x(2)*fX - x(3);
fe4da5cc 2630 Double_t ex=sqrt(1-ey*ey);
2631 Double_t ez=x(4);
2632 TMatrix Q(5,5); Q=0.;
2633 Q(2,2)=ez*ez+ey*ey; Q(2,3)=-ex*ey; Q(2,4)=-ex*ez;
2634 Q(3,2)=Q(2,3); Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
2635 Q(4,2)=Q(2,4); Q(4,3)= Q(3,4); Q(4,4)=1.;
2636
2637 F=0;
2638 F(2,2)=-x(2)*ex; F(2,3)=-x(2)*ey;
3c0f9266 2639 F(3,2)=-ex*(x(2)*fX-ey); F(3,3)=-(1.+ x(2)*fX*ey - ey*ey);
fe4da5cc 2640 F(4,2)=-ez*ex; F(4,3)=-ez*ey; F(4,4)=1.;
2641
2642 tmp.Mult(F,Q);
2643 Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2644
2645 Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
2646 Double_t beta2=p2/(p2 + pm*pm);
3c0f9266 2647 Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
fe4da5cc 2648 d*=2.;
2649 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
2650 Q*=theta2;
2651 C+=Q;
2652
2653 //Energy losses************************
2654 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
2655 if (x1 < x2) dE=-dE;
2656 x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
8c555625 2657 //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE);
fe4da5cc 2658
3c0f9266 2659 x1=fX; x2=xk; y1=x(0); z1=x(1);
fe4da5cc 2660 c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1);
2661 c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2);
2662
2663 x(0) += dx*(c1+c2)/(r1+r2);
2664 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2665
2666 F.UnitMatrix();
2667 rr=r1+r2; cc=c1+c2; xx=x1+x2;
2668 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2669 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2670 cr=c1*r2+c2*r1;
2671 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2672 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2673 F(1,4)= dx*cc/cr;
2674 tmp.Mult(F,C);
2675 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2676
3c0f9266 2677 fX=x2;
fe4da5cc 2678
2679 return 1;
2680}
2681
2682//_____________________________________________________________________________
2683void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
2684{
3c0f9266 2685 //-----------------------------------------------------------------
2686 // This function propagates tracks to the "vertex".
fe4da5cc 2687 //
3c0f9266 2688 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2689 //-----------------------------------------------------------------
2690 Double_t c=x(2)*fX - x(3);
fe4da5cc 2691 Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
2692 Double_t snf=tgf/sqrt(1.+ tgf*tgf);
2693 Double_t xv=(x(3)+snf)/x(2);
2694 PropagateTo(xv,x0,rho,pm);
2695}
2696
2697//_____________________________________________________________________________
2698void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
2699{
3c0f9266 2700 //-----------------------------------------------------------------
2701 // This function associates a clusters with this track.
fe4da5cc 2702 //
3c0f9266 2703 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2704 //-----------------------------------------------------------------
fe4da5cc 2705 TMatrix H(2,5); H.UnitMatrix();
2706 TMatrix Ht(TMatrix::kTransposed,H);
2707 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
2708 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2709
2710 TMatrix tmp(H,TMatrix::kMult,C);
2711 TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
2712
2713 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2714 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2715 R(1,0)*=-1; R(0,1)=R(1,0);
2716 R*=1./det;
2717
2718 //R.Invert();
2719
2720 TMatrix K(C,TMatrix::kMult,Ht); K*=R;
2721
2722 TVector savex=x;
2723 x*=H; x-=m; x*=-1; x*=K; x+=savex;
3c0f9266 2724 if (TMath::Abs(x(2)*fX-x(3)) >= 0.999) {
2725 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
fe4da5cc 2726 x=savex;
2727 return;
2728 }
2729
2730 TMatrix saveC=C;
2731 C.Mult(K,tmp); C-=saveC; C*=-1;
2732
3c0f9266 2733 fClusters.AddLast((AliTPCcluster*)c);
2734 fChi2 += chisq;
fe4da5cc 2735}
2736
2737//_____________________________________________________________________________
2738int AliTPCtrack::Rotate(Double_t alpha)
2739{
3c0f9266 2740 //-----------------------------------------------------------------
2741 // This function rotates this track.
fe4da5cc 2742 //
3c0f9266 2743 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2744 //-----------------------------------------------------------------
fe4da5cc 2745 fAlpha += alpha;
2746
3c0f9266 2747 Double_t x1=fX, y1=x(0);
fe4da5cc 2748 Double_t ca=cos(alpha), sa=sin(alpha);
3c0f9266 2749 Double_t r1=x(2)*fX - x(3);
fe4da5cc 2750
3c0f9266 2751 fX = x1*ca + y1*sa;
fe4da5cc 2752 x(0)=-x1*sa + y1*ca;
2753 x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
2754
3c0f9266 2755 Double_t r2=x(2)*fX - x(3);
fe4da5cc 2756 if (TMath::Abs(r2) >= 0.999) {
3c0f9266 2757 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
fe4da5cc 2758 return 0;
2759 }
2760
2761 Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
2762 if ((x(0)-y0)*x(2) >= 0.) {
3c0f9266 2763 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
fe4da5cc 2764 return 0;
2765 }
2766
2767 TMatrix F(5,5); F.UnitMatrix();
2768 F(0,0)=ca;
2769 F(3,0)=x(2)*sa;
2770 F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa;
2771 F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
2772 TMatrix tmp(F,TMatrix::kMult,C);
2773 // Double_t dy2=C(0,0);
2774 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2775 // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
2776 // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
2777
2778 return 1;
2779}
2780
2781//_____________________________________________________________________________
2782void AliTPCtrack::UseClusters() const
2783{
3c0f9266 2784 //-----------------------------------------------------------------
2785 // This function marks clusters associated with this track.
fe4da5cc 2786 //
3c0f9266 2787 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2788 //-----------------------------------------------------------------
2789 int num_of_clusters=fClusters.GetEntriesFast();
fe4da5cc 2790 for (int i=0; i<num_of_clusters; i++) {
2791 //if (i<=14) continue;
3c0f9266 2792 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
fe4da5cc 2793 c->Use();
2794 }
2795}
2796
2797//_____________________________________________________________________________
2798Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const
2799{
3c0f9266 2800 //-----------------------------------------------------------------
2801 // This function calculates a predicted chi2 increment.
fe4da5cc 2802 //
3c0f9266 2803 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2804 //-----------------------------------------------------------------
fe4da5cc 2805 TMatrix H(2,5); H.UnitMatrix();
2806 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
2807 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2808 TVector res=x; res*=H; res-=m; //res*=-1;
2809 TMatrix tmp(H,TMatrix::kMult,C);
2810 TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
2811
2812 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2813 if (TMath::Abs(det) < 1.e-10) {
3c0f9266 2814 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
fe4da5cc 2815 return 1e10;
2816 }
2817 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2818 R(1,0)*=-1; R(0,1)=R(1,0);
2819 R*=1./det;
2820
2821 //R.Invert();
2822
2823 TVector r=res;
2824 res*=R;
2825 return r*res;
2826}
2827
2828//_____________________________________________________________________________
3c0f9266 2829struct S { int lab; int max; };
2830int AliTPCtrack::GetLabel(int nrows) const
fe4da5cc 2831{
3c0f9266 2832 //-----------------------------------------------------------------
2833 // This function returns the track label. If label<0, this track is fake.
fe4da5cc 2834 //
3c0f9266 2835 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2836 //-----------------------------------------------------------------
2837 int num_of_clusters=fClusters.GetEntriesFast();
2838 S *s=new S[num_of_clusters];
fe4da5cc 2839 int i;
3c0f9266 2840 for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
2841
2842 int lab=123456789;
fe4da5cc 2843 for (i=0; i<num_of_clusters; i++) {
3c0f9266 2844 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
8c555625 2845 lab=TMath::Abs(c->fTracks[0]);
fe4da5cc 2846 int j;
3c0f9266 2847 for (j=0; j<num_of_clusters; j++)
fe4da5cc 2848 if (s[j].lab==lab || s[j].max==0) break;
2849 s[j].lab=lab;
2850 s[j].max++;
2851 }
2852
2853 int max=0;
2854 for (i=0; i<num_of_clusters; i++)
2855 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
3c0f9266 2856
2857 delete[] s;
2858
8c555625 2859 for (i=0; i<num_of_clusters; i++) {
3c0f9266 2860 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
8c555625 2861 if (TMath::Abs(c->fTracks[1]) == lab ||
2862 TMath::Abs(c->fTracks[2]) == lab ) max++;
2863 }
fe4da5cc 2864
2865 if (1.-float(max)/num_of_clusters > 0.10) return -lab;
2866
3c0f9266 2867 int tail=int(0.08*nrows);
2868 if (num_of_clusters < tail) return lab;
fe4da5cc 2869
2870 max=0;
3c0f9266 2871 for (i=1; i<=tail; i++) {
2872 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(num_of_clusters-i);
8c555625 2873 if (lab == TMath::Abs(c->fTracks[0]) ||
2874 lab == TMath::Abs(c->fTracks[1]) ||
2875 lab == TMath::Abs(c->fTracks[2])) max++;
fe4da5cc 2876 }
3c0f9266 2877 if (max < int(0.5*tail)) return -lab;
fe4da5cc 2878
2879 return lab;
2880}
2881
2882//_____________________________________________________________________________
2883void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
2884{
3c0f9266 2885 //-----------------------------------------------------------------
2886 // This function returns reconstructed track momentum in the global system.
fe4da5cc 2887 //
3c0f9266 2888 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2889 //-----------------------------------------------------------------
2890 Double_t pt=TMath::Abs(GetPt()); // GeV/c
2891 Double_t r=x(2)*fX-x(3);
fe4da5cc 2892 Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2893 px=-pt*(x(0)-y0)*x(2); //cos(phi);
3c0f9266 2894 py=-pt*(x(3)-fX*x(2)); //sin(phi);
fe4da5cc 2895 pz=pt*x(4);
2896 Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2897 py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2898 px=tmp;
2899}
2900
2901//_____________________________________________________________________________
3c0f9266 2902Double_t AliTPCtrack::GetdEdX(Double_t low, Double_t up) const {
2903 //-----------------------------------------------------------------
2904 // This funtion calculates dE/dX within the "low" and "up" cuts.
fe4da5cc 2905 //
3c0f9266 2906 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2907 //-----------------------------------------------------------------
2908 int ncl=fClusters.GetEntriesFast();
2909 int n=0;
2910 Double_t *q=new Double_t[ncl];
2911 int i;
2912 for (i=0; i<ncl; i++) {
2913 AliTPCcluster *cl=(AliTPCcluster*)(fClusters.UncheckedAt(i));
2914 // if (cl->fdEdX > 3000) continue;
2915 if (cl->fdEdX <= 0) continue;
2916 q[n++]=cl->fdEdX;
2917 }
2918
2919 //stupid sorting
2920 int swap;
2921 do {
2922 swap=0;
2923 for (i=0; i<n-1; i++) {
2924 if (q[i]<=q[i+1]) continue;
2925 Double_t tmp=q[i]; q[i]=q[i+1]; q[i+1]=tmp;
2926 swap++;
2927 }
2928 } while (swap);
2929
2930 int nl=int(low*n), nu=int(up *n);
2931 Double_t dedx=0.;
2932 for (i=nl; i<=nu; i++) dedx += q[i];
2933 dedx /= (nu-nl+1);
2934 return dedx;
2935}
2936
2937//_________________________________________________________________________
2938//
2939// Classes for internal tracking use
2940//_________________________________________________________________________
2941void AliTPCRow::InsertCluster(const AliTPCcluster* c) {
2942 //-----------------------------------------------------------------------
2943 // Insert a cluster into this pad row in accordence with its y-coordinate
2944 //
2945 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2946 //-----------------------------------------------------------------------
fe4da5cc 2947 if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2948 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2949 }
2950 if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2951 int i=Find(c->fY);
2952 memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2953 clusters[i]=c; num_of_clusters++;
2954}
2955
3c0f9266 2956int AliTPCRow::Find(Double_t y) const {
2957 //-----------------------------------------------------------------------
2958 // Return the index of the nearest cluster
fe4da5cc 2959 //
3c0f9266 2960 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2961 //-----------------------------------------------------------------------
fe4da5cc 2962 if (y <= clusters[0]->fY) return 0;
2963 if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2964 int b=0, e=num_of_clusters-1, m=(b+e)/2;
2965 for (; b<e; m=(b+e)/2) {
2966 if (y > clusters[m]->fY) b=m+1;
2967 else e=m;
2968 }
2969 return m;
2970}
8c555625 2971