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