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