1 ///////////////////////////////////////////////////////////////////////////////
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 //
10 <img src="gif/AliTPCClass.gif">
15 ///////////////////////////////////////////////////////////////////////////////
21 #include <TGeometry.h>
24 #include <TObjectTable.h>
25 #include "GParticle.h"
33 #include "AliTPCParam.h"
35 #include "AliTPCPRF2D.h"
36 #include "AliTPCRF1D.h"
41 //_____________________________________________________________________________
45 // Default constructor
56 fDigParam= new AliTPCD();
57 fDigits = fDigParam->GetArray();
60 //_____________________________________________________________________________
61 AliTPC::AliTPC(const char *name, const char *title)
62 : AliDetector(name,title)
65 // Standard constructor
69 // Initialise arrays of hits and digits
70 fHits = new TClonesArray("AliTPChit", 176);
71 // fDigits = new TClonesArray("AliTPCdigit",10000);
73 fDigParam= new AliTPCD;
74 fDigits = fDigParam->GetArray();
76 // Initialise counters
82 fDigitsIndex = new Int_t[fNsectors+1];
83 fClustersIndex = new Int_t[fNsectors+1];
87 // Initialise color attributes
88 SetMarkerColor(kYellow);
91 //_____________________________________________________________________________
103 if (fDigitsIndex) delete [] fDigitsIndex;
104 if (fClustersIndex) delete [] fClustersIndex;
107 //_____________________________________________________________________________
108 void AliTPC::AddCluster(Float_t *hits, Int_t *tracks)
111 // Add a simulated cluster to the list
113 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
114 TClonesArray &lclusters = *fClusters;
115 new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks);
118 //_____________________________________________________________________________
119 void AliTPC::AddCluster(const AliTPCcluster &c)
122 // Add a simulated cluster copy to the list
124 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
125 TClonesArray &lclusters = *fClusters;
126 new(lclusters[fNclusters++]) AliTPCcluster(c);
129 //_____________________________________________________________________________
130 void AliTPC::AddDigit(Int_t *tracks, Int_t *digits)
133 // Add a TPC digit to the list
135 // TClonesArray &ldigits = *fDigits;
137 TClonesArray &ldigits = *fDigParam->GetArray();
138 new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits);
141 //_____________________________________________________________________________
142 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
145 // Add a hit to the list
147 TClonesArray &lhits = *fHits;
148 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
151 //_____________________________________________________________________________
152 void AliTPC::AddTrack(Float_t *hits)
155 // Add a track to the list of tracks
157 TClonesArray <racks = *fTracks;
158 new(ltracks[fNtracks++]) AliTPCtrack(hits);
161 //_____________________________________________________________________________
162 void AliTPC::AddTrack(const AliTPCtrack& t)
165 // Add a track copy to the list of tracks
167 if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000);
168 TClonesArray <racks = *fTracks;
169 new(ltracks[fNtracks++]) AliTPCtrack(t);
172 //_____________________________________________________________________________
173 void AliTPC::BuildGeometry()
176 // Build TPC ROOT TNode geometry for the event display
181 const int kColorTPC=19;
182 char name[5], title[20];
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);
191 // Get ALICE top node
192 Top=gAlice->GetGeometry()->GetNode("alice");
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);
201 Node = new TNode(name,title,name,0,0,0,"");
202 Node->SetLineColor(kColorTPC);
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);
212 Node = new TNode(name,title,name,0,0,0,"");
213 Node->SetLineColor(kColorTPC);
217 //_____________________________________________________________________________
218 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
221 // Calculate distance from TPC to mouse on the display
227 //_____________________________________________________________________________
228 //const int MAX_CLUSTER=nrow_low+nrow_up;
229 const int MAX_CLUSTER=200;
230 const int S_MAXSEC=24;
231 const int L_MAXSEC=48;
232 const int ROWS_TO_SKIP=21;
233 const Float_t MAX_CHI2=12.;
236 //_____________________________________________________________________________
237 static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
240 // Calculate spread in Y
242 pt=TMath::Abs(pt)*1000.;
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;
250 //_____________________________________________________________________________
251 static Double_t SigmaZ2(Double_t r, Double_t tgl)
254 // Calculate spread in Z
257 Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl;
258 if (s<0.4e-3) s=0.4e-3;
262 //_____________________________________________________________________________
263 inline Double_t f1(Double_t x1,Double_t y1, //C
264 Double_t x2,Double_t y2,
265 Double_t x3,Double_t y3)
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));
276 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
278 return -xr*yr/sqrt(xr*xr+yr*yr);
282 //_____________________________________________________________________________
283 inline 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)
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));
296 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
298 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
301 //_____________________________________________________________________________
302 inline Double_t f3(Double_t x1,Double_t y1, //tgl
303 Double_t x2,Double_t y2,
304 Double_t z1,Double_t z2)
309 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
312 //_____________________________________________________________________________
313 static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec,
314 int s, int ri, int rf=0)
319 int try_again=ROWS_TO_SKIP;
320 Double_t alpha=sec->GetAlpha();
321 int ns=int(2*TMath::Pi()/alpha)+1;
323 for (int nr=ri; nr>=rf; nr--) {
324 Double_t x=sec[s].GetX(nr), ymax=sec[s].GetMaxY(nr);
325 if (!t.PropagateTo(x)) return -1;
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());
332 Double_t road=3.*sqrt(t.GetSigmaY2() + 4*sy2), y=t.GetY(), z=t.GetZ();
335 if (t>3) cerr<<t<<" AliTPCtrack warning: Too broad road !\n";
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;
344 if ((c->fZ - z)*(c->fZ - z) > 9.*(t.GetSigmaZ2() + 4*sz2)) continue;
345 Double_t chi2=t.GetPredictedChi2(c);
346 if (chi2 > max_chi2) continue;
352 t.Update(cl,max_chi2);
353 try_again=ROWS_TO_SKIP;
355 if (try_again==0) break;
358 if (!t.Rotate(alpha)) return -1;
359 } else if (y <-ymax) {
361 if (!t.Rotate(-alpha)) return -1;
371 //_____________________________________________________________________________
372 static void MakeSeeds(TObjArray& seeds,const AliTPCSector* sec,int i1,int i2,
373 const AliTPCParam *p)
376 // Find seed for tracking
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];
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;
394 ks=(ns-1+max_sec)%max_sec;
395 const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2];
397 cs=cos(alpha); sn=sin(alpha);
401 const AliTPCRow& r2=sec[ns][i2];
406 const AliTPCRow& r2=sec[(ns+1)%max_sec][i2];
408 cs=cos(alpha); sn=-sin(alpha);
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;
414 Double_t tmp= x2*cs+y2*sn;
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);
424 if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue;
426 if (TMath::Abs(x(4)) > 1.2) continue;
428 Double_t a=asin(x(3));
429 Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3)));
430 if (TMath::Abs(zv)>33.) continue;
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;
451 TMatrix t(F,TMatrix::kMult,X);
452 C.Mult(t,TMatrix(TMatrix::kTransposed,F));
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);
463 //_____________________________________________________________________________
464 void AliTPC::Clusters2Tracks()
467 // TPC Track finder from clusters.
469 if (!fClusters) return;
471 AliTPCParam *p=&fDigParam->GetParam();
472 Int_t nrow_low=p->GetNRowLow();
473 Int_t nrow_up=p->GetNRowUp();
475 AliTPCSSector ssec[S_MAXSEC/2];
476 for (int i=0; i<S_MAXSEC/2; i++) ssec[i].SetUp(p);
478 AliTPCLSector lsec[L_MAXSEC/2];
479 for (int j=0; j<L_MAXSEC/2; j++) lsec[j].SetUp(p);
481 int ncl=fClusters->GetEntriesFast();
483 AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl);
485 int sec=int(c->fSector)-1, row=int(c->fPadRow)-1;
488 if (row<0 || row>nrow_low) {cerr<<"low !!!"<<row<<endl; continue;}
489 ssec[sec%12][row].InsertCluster(c);
491 if (row<0 || row>nrow_up ) {cerr<<"up !!!"<<row<<endl; continue;}
493 lsec[sec%24][row].InsertCluster(c);
498 TObjArray seeds(20000);
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);
505 int nseed=seeds.GetEntriesFast();
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);
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;
518 else {cerr<<x<<" =x !!!\n"; continue;}
520 int ls=FindProlongation(t,lsec,ns,nr-1);
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); //
529 if (FindProlongation(t,ssec,ss,nrow_low-1)<0) continue;
530 if (t < 30) continue;
538 //_____________________________________________________________________________
539 void AliTPC::CreateMaterials()
541 //-----------------------------------------------
542 // Create Materials for for TPC
543 //-----------------------------------------------
545 //-----------------------------------------------------------------
546 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
547 //-----------------------------------------------------------------
549 AliMC* pMC = AliMC::GetMC();
551 Int_t ISXFLD=gAlice->Field()->Integ();
552 Float_t SXMGMX=gAlice->Field()->Max();
554 Float_t absl, radl, a, d, z;
560 // --- Methane (CH4) ---
561 Float_t am[2] = { 12.,1. };
562 Float_t zm[2] = { 6.,1. };
563 Float_t wm[2] = { 1.,4. };
564 Float_t dm = 7.17e-4;
565 // --- The Neon CO2 90/10 mixture ---
566 Float_t ag[2] = { 20.18 };
567 Float_t zg[2] = { 10. };
568 Float_t wg[2] = { .8,.2 };
569 Float_t dne = 9e-4; // --- Neon density in g/cm3 ---
571 // --- Mylar (C5H4O2) ---
572 Float_t amy[3] = { 12.,1.,16. };
573 Float_t zmy[3] = { 6.,1.,8. };
574 Float_t wmy[3] = { 5.,4.,2. };
577 Float_t ac[2] = { 12.,16. };
578 Float_t zc[2] = { 6.,8. };
579 Float_t wc[2] = { 1.,2. };
580 Float_t dc = .001977;
581 // --- Carbon density and radiation length ---
582 Float_t densc = 2.265;
583 Float_t radlc = 18.8;
588 Float_t radsi = 9.36;
590 // --- Define the various materials for GEANT ---
591 AliMaterial(0, "Al $", 26.98, 13., 2.7, 8.9, 37.2);
593 AliMaterial(1, "Ne $", 20.18, 10., dne, x0ne, 999.);
595 // -- Methane, defined by the proportions of atoms
597 AliMixture(2, "Methane$", am, zm, dm, -2, wm);
599 // --- CO2, defined by the proportion of atoms
601 AliMixture(7, "CO2$", ac, zc, dc, -2, wc);
603 // -- Get A,Z etc. for CO2
606 pMC->Gfmate((*fIdmate)[7], namate, a, z, d, radl, absl, buf, nbuf);
609 dg = dne * .9 + dc * .1;
611 // -- Create Ne/CO2 90/10 mixture
613 AliMixture(3, "Gas-mixt $", ag, zg, dg, 2, wg);
614 AliMixture(4, "Gas-mixt $", ag, zg, dg, 2, wg);
616 AliMaterial(5, "G10$", 20., 10., 1.7, 19.4, 999.);
617 AliMixture(6, "Mylar$", amy, zmy, dmy, -3, wmy);
621 AliMaterial(8, "Carbon", a, z, densc, radlc, 999.);
623 AliMaterial(9, "Silicon", asi, zsi, desi, radsi, 999.);
624 AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500.);
626 AliMedium(400, "Al wall$", 0, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
627 AliMedium(402, "Gas mix1$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
628 AliMedium(403, "Gas mix2$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
629 AliMedium(404, "Gas mix3$", 4, 1, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
630 AliMedium(405, "G10 pln$", 5, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
631 AliMedium(406, "Mylar $", 6, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
632 AliMedium(407, "CO2 $", 7, 0, ISXFLD, SXMGMX, 10., .01,.1, .01, .01);
633 AliMedium(408, "Carbon $", 8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
634 AliMedium(409, "Silicon$", 9, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
635 AliMedium(499, "Air gap$", 99, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
638 //_____________________________________________________________________________
640 const AliTPCdigit *dig;
642 Bin() {dig=0; idx=-1;}
645 struct PreCluster : public AliTPCcluster {
646 const AliTPCdigit* summit;
652 PreCluster::PreCluster() : AliTPCcluster() {cut=npeaks=0;}
655 //_____________________________________________________________________________
656 static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c)
662 double q=double(b.dig->fSignal);
664 if (q<0) { // digit is at the edge of the pad row
668 if (b.idx >= 0 && b.idx != c.idx) {
673 if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
681 b.dig = 0; b.idx = c.idx;
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 if (bins[i+1][j].dig) FindCluster(i+1,j,bins,c);
686 if (bins[i][j+1].dig) FindCluster(i,j+1,bins,c);
690 //_____________________________________________________________________________
691 void AliTPC::Digits2Clusters()
694 // simple TPC cluster finder from digits.
697 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
699 const Int_t MAX_PAD=200+2, MAX_BUCKET=MAXTPCTBK+2;
700 const Int_t Q_min=60;
701 const Int_t THRESHOLD=20;
703 TTree *t=(TTree*)gDirectory->Get("TreeD0_Param1");
704 t->GetBranch("Digits")->SetAddress(&fDigits);
705 Int_t sectors_by_rows=(Int_t)t->GetEntries();
709 for (Int_t n=0; n<sectors_by_rows; n++) {
710 if (!t->GetEvent(n)) continue;
711 Bin bins[MAX_PAD][MAX_BUCKET];
712 AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
713 Int_t nsec=dig->fSector, nrow=dig->fPadRow;
714 Int_t ndigits=fDigits->GetEntriesFast();
716 int npads; int sign_z;
718 sign_z=(nsec<13) ? 1 : -1;
719 npads=fTPCParam->GetNPadsLow(nrow-1);
721 sign_z=(nsec<49) ? 1 : -1;
722 npads=fTPCParam->GetNPadsUp(nrow-1);
726 for (ndig=0; ndig<ndigits; ndig++) {
727 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
728 int i=dig->fPad, j=dig->fTime;
729 if (dig->fSignal >= THRESHOLD) bins[i][j].dig=dig;
730 if (i==1 || i==npads || j==1 || j==MAXTPCTBK) dig->fSignal*=-1;
736 for (i=1; i<MAX_PAD-1; i++) {
737 for (j=1; j<MAX_BUCKET-1; j++) {
738 if (bins[i][j].dig == 0) continue;
739 PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
740 FindCluster(i, j, bins, c);
744 double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
745 c.fSigmaY2 = s2 + 1./12.;
746 c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()*
747 fTPCParam->GetPadPitchWidth();
748 if (s2 != 0.) c.fSigmaY2 *= 0.022*8*4;
750 s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
751 c.fSigmaZ2 = s2 + 1./12.;
752 c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth();
753 if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*4;
755 c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth();
756 c.fZ = fTPCParam->GetZWidth()*(c.fZ+1);
757 c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay
758 c.fZ = sign_z*(z_end - c.fZ);
762 c.fTracks[0]=c.summit->fTracks[0];
763 c.fTracks[1]=c.summit->fTracks[1];
764 c.fTracks[2]=c.summit->fTracks[2];
771 AddCluster(c); ncls++; ncl++;
775 for (ndig=0; ndig<ndigits; ndig++) {
776 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
777 if (TMath::Abs(dig->fSignal) >= 0)
778 bins[dig->fPad][dig->fTime].dig=dig;
781 for (i=1; i<MAX_PAD-1; i++) {
782 for (j=1; j<MAX_BUCKET-1; j++) {
783 if (bins[i][j].dig == 0) continue;
784 PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
785 FindCluster(i, j, bins, c);
786 if (c.fQ <= Q_min) continue; //noise cluster
787 if (c.npeaks>1) continue; //overlapped cluster
791 double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
792 c.fSigmaY2 = s2 + 1./12.;
793 c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()*
794 fTPCParam->GetPadPitchWidth();
795 if (s2 != 0.) c.fSigmaY2 *= 0.022*4*0.6*4;
797 s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
798 c.fSigmaZ2 = s2 + 1./12.;
799 c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth();
800 if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*0.4;
802 c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth();
803 c.fZ = fTPCParam->GetZWidth()*(c.fZ+1);
804 c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay
805 c.fZ = sign_z*(z_end - c.fZ);
809 c.fTracks[0]=c.summit->fTracks[0];
810 c.fTracks[1]=c.summit->fTracks[1];
811 c.fTracks[2]=c.summit->fTracks[2];
818 if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;}
820 new ((*fClusters)[c.idx]) AliTPCcluster(c);
825 cerr<<"sector, row, digits, clusters: "
826 <<nsec<<' '<<nrow<<' '<<ndigits<<' '<<ncl<<" \r";
833 //_____________________________________________________________________________
834 void AliTPC::ElDiff(Float_t *xyz)
836 //--------------------------------------------------
837 // calculates the diffusion of a single electron
838 //--------------------------------------------------
840 //-----------------------------------------------------------------
841 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
842 //-----------------------------------------------------------------
843 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
848 driftl=z_end-TMath::Abs(xyz[2]);
850 if(driftl<0.01) driftl=0.01;
852 // check the attachment
854 driftl=TMath::Sqrt(driftl);
856 // Float_t sig_t = driftl*diff_t;
857 //Float_t sig_l = driftl*diff_l;
858 Float_t sig_t = driftl*fTPCParam->GetDiffT();
859 Float_t sig_l = driftl*fTPCParam->GetDiffL();
860 xyz[0]=gRandom->Gaus(xyz[0],sig_t);
861 xyz[1]=gRandom->Gaus(xyz[1],sig_t);
862 xyz[2]=gRandom->Gaus(xyz[2],sig_l);
864 if (TMath::Abs(xyz[2])>z_end){
865 xyz[2]=TMath::Sign(z_end,z0);
868 Float_t eps = 0.0001;
869 xyz[2]=TMath::Sign(eps,z0);
873 //_____________________________________________________________________________
874 void AliTPC::Hits2Clusters()
876 //--------------------------------------------------------
877 // TPC simple cluster generator from hits
878 // obtained from the TPC Fast Simulator
879 // The point errors are taken from the parametrization
880 //--------------------------------------------------------
882 //-----------------------------------------------------------------
883 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
884 //-----------------------------------------------------------------
885 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
886 Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
888 GParticle *particle; // pointer to a given particle
889 AliTPChit *tpcHit; // pointer to a sigle TPC hit
890 TClonesArray *Particles; //pointer to the particle list
894 Float_t pl,pt,tanth,rpad,ratio;
897 //---------------------------------------------------------------
898 // Get the access to the tracks
899 //---------------------------------------------------------------
901 TTree *TH = gAlice->TreeH();
902 Stat_t ntracks = TH->GetEntries();
904 //------------------------------------------------------------
905 // Loop over all sectors (72 sectors)
906 // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
907 // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
909 // First cluster for sector 1 starts at "0"
910 //------------------------------------------------------------
913 fClustersIndex[0] = 0;
916 for(Int_t isec=1;isec<fNsectors+1;isec++){
918 fTPCParam->AdjustAngles(isec,cph,sph);
920 //------------------------------------------------------------
922 //------------------------------------------------------------
924 for(Int_t track=0;track<ntracks;track++){
928 // Get number of the TPC hits and a pointer
931 nhits=fHits->GetEntriesFast();
932 Particles=gAlice->Particles();
936 for(Int_t hit=0;hit<nhits;hit++){
937 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
938 sector=tpcHit->fSector; // sector number
939 if(sector != isec) continue; //terminate iteration
940 ipart=tpcHit->fTrack;
941 particle=(GParticle*)Particles->UncheckedAt(ipart);
942 pl=particle->GetPz();
943 pt=particle->GetPT();
944 if(pt < 1.e-9) pt=1.e-9;
946 tanth = TMath::Abs(tanth);
947 rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
948 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
950 // space-point resolutions
952 sigma_rphi=SigmaY2(rpad,tanth,pt);
953 sigma_z =SigmaZ2(rpad,tanth );
957 cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
958 cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
960 // temporary protection
962 if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
963 if(sigma_z < 0.) sigma_z=0.4e-3;
964 if(cl_rphi < 0.) cl_rphi=2.5e-3;
965 if(cl_z < 0.) cl_z=2.5e-5;
970 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
971 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
973 Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
974 Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
975 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi)); // y
976 Double_t alpha=(sector<25) ? alpha_low : alpha_up;
977 if (TMath::Abs(xyz[0]/xprim) > TMath::Tan(0.5*alpha)) xyz[0]=yprim;
978 xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z
979 if (TMath::Abs(xyz[1]) > 250) xyz[1]=tpcHit->fZ;
980 xyz[2]=tpcHit->fQ; // q
981 xyz[3]=sigma_rphi; // fSigmaY2
982 xyz[4]=sigma_z; // fSigmaZ2
986 Int_t row = fTPCParam->GetPadRow(sector,xprim) ;
987 // and finally add the cluster
988 Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, row+1};
989 AddCluster(xyz,tracks);
991 } // end of loop over hits
992 } // end of loop over tracks
994 fClustersIndex[isec] = fNclusters; // update clusters index
996 } // end of loop over sectors
998 fClustersIndex[fNsectors]--; // set end of the clusters buffer
1003 void AliTPC::Hits2Digits()
1006 //----------------------------------------------------
1007 // Loop over all sectors (72 sectors)
1008 // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
1009 // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
1011 for(Int_t isec=1;isec<fNsectors+1;isec++) Hits2DigitsSector(isec);
1015 //_____________________________________________________________________________
1016 void AliTPC::Hits2DigitsSector(Int_t isec)
1018 //-------------------------------------------------------------------
1019 // TPC conversion from hits to digits.
1020 //-------------------------------------------------------------------
1022 //-----------------------------------------------------------------
1023 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1024 //-----------------------------------------------------------------
1026 //-------------------------------------------------------
1027 // Get the access to the track hits
1028 //-------------------------------------------------------
1030 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1031 TTree *TH = gAlice->TreeH(); // pointer to the hits tree
1032 Stat_t ntracks = TH->GetEntries();
1036 //-------------------------------------------
1037 // Only if there are any tracks...
1038 //-------------------------------------------
1041 // TObjArrays for three neighouring pad-rows
1043 TObjArray **rowTriplet = new TObjArray* [3];
1045 // TObjArray-s for each pad-row
1049 for (Int_t trip=0;trip<3;trip++){
1050 rowTriplet[trip]=new TObjArray;
1055 printf("*** Processing sector number %d ***\n",isec);
1057 Int_t nrows =fTPCParam->GetNRow(isec);
1059 row= new TObjArray* [nrows];
1061 MakeSector(isec,nrows,TH,ntracks,row);
1063 //--------------------------------------------------------
1064 // Digitize this sector, row by row
1065 // row[i] is the pointer to the TObjArray of TVectors,
1066 // each one containing electrons accepted on this
1067 // row, assigned into tracks
1068 //--------------------------------------------------------
1072 for (i=0;i<nrows;i++){
1074 // Triplets for i = 0 and i=1 are identical!
1075 // The same for i = nrows-1 and nrows!
1077 if(i != 1 && i != nrows-1){
1078 MakeTriplet(i,rowTriplet,row);
1081 DigitizeRow(i,isec,rowTriplet);
1085 Int_t ndig=fDigParam->GetArray()->GetEntriesFast();
1087 printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig);
1089 ResetDigits(); // reset digits for this row after storing them
1091 } // end of the sector digitization
1093 // delete the last triplet
1095 for (i=0;i<3;i++) rowTriplet[i]->Delete();
1097 delete [] row; // delete the array of pointers to TObjArray-s
1100 } // end of Hits2Digits
1101 //_____________________________________________________________________________
1102 void AliTPC::MakeTriplet(Int_t row,
1103 TObjArray **rowTriplet, TObjArray **prow)
1105 //------------------------------------------------------------------
1106 // Makes the "triplet" of the neighbouring pad-row for the
1107 // digitization including the cross-talk between the pad-rows
1108 //------------------------------------------------------------------
1110 //-----------------------------------------------------------------
1111 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1112 //-----------------------------------------------------------------
1114 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1115 Float_t gasgain = fTPCParam->GetGasGain();
1118 Int_t nElements,nElectrons;
1123 //-------------------------------------------------------------------
1124 // pv is an "old" track, i.e. label + triplets of (x,y,z)
1125 // for each electron
1127 //-------------------------------------------------------------------
1133 if(row == 0 || row == 1){
1135 // create entire triplet for the first AND the second row
1137 nTracks[0] = prow[0]->GetEntries();
1138 nTracks[1] = prow[1]->GetEntries();
1139 nTracks[2] = prow[2]->GetEntries();
1141 for(i1=0;i1<3;i1++){
1142 nt = nTracks[i1]; // number of tracks for this row
1144 for(i2=0;i2<nt;i2++){
1145 pv = (TVector*)prow[i1]->At(i2);
1147 nElements = pv->GetNrows();
1148 nElectrons = (nElements-1)/3;
1150 tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1152 v2(0)=v1(0); //track label
1154 for(nel=0;nel<nElectrons;nel++){
1157 // Avalanche, including fluctuations
1158 Int_t aval = (Int_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1159 v2(idx2+1)= v1(idx1+1);
1160 v2(idx2+2)= v1(idx1+2);
1161 v2(idx2+3)= v1(idx1+3);
1162 v2(idx2+4)= (Float_t)aval; // in number of electrons!
1163 } // end of loop over electrons
1165 // Add this track to a row
1168 rowTriplet[i1]->Add(tv);
1171 } // end of loop over tracks for this row
1173 prow[i1]->Delete(); // remove "old" tracks
1174 delete prow[i1]; // delete TObjArray for this row
1175 prow[i1]=0; // set pointer to NULL
1177 } // end of loop over row triplets
1183 rowTriplet[0]->Delete(); // remove old lower row
1185 nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row
1186 nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row
1187 nTracks[2]=prow[row+1]->GetEntries(); // next row
1190 //-------------------------------------------
1191 // shift new tracks downwards
1192 //-------------------------------------------
1194 for(i1=0;i1<nTracks[0];i1++){
1195 pv=(TVector*)rowTriplet[1]->At(i1);
1196 rowTriplet[0]->Add(pv);
1199 rowTriplet[1]->Clear(); // leave tracks on the heap!!!
1201 for(i1=0;i1<nTracks[1];i1++){
1202 pv=(TVector*)rowTriplet[2]->At(i1);
1203 rowTriplet[1]->Add(pv);
1206 rowTriplet[2]->Clear(); // leave tracks on the heap!!!
1208 //---------------------------------------------
1209 // Create new upper row
1210 //---------------------------------------------
1214 for(i1=0;i1<nTracks[2];i1++){
1215 pv = (TVector*)prow[row+1]->At(i1);
1217 nElements = pv->GetNrows();
1218 nElectrons = (nElements-1)/3;
1220 tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1222 v2(0)=v1(0); //track label
1224 for(nel=0;nel<nElectrons;nel++){
1228 // Avalanche, including fluctuations
1229 Int_t aval = (Int_t)
1230 (-gasgain*TMath::Log(gRandom->Rndm()));
1232 v2(idx2+1)= v1(idx1+1);
1233 v2(idx2+2)= v1(idx1+2);
1234 v2(idx2+3)= v1(idx1+3);
1235 v2(idx2+4)= (Float_t)aval; // in number of electrons!
1236 } // end of loop over electrons
1238 rowTriplet[2]->Add(tv);
1240 } // end of loop over tracks
1242 prow[row+1]->Delete(); // delete tracks for this row
1243 delete prow[row+1]; // delete TObjArray for this row
1244 prow[row+1]=0; // set a pointer to NULL
1248 } // end of MakeTriplet
1249 //_____________________________________________________________________________
1250 void AliTPC::ExB(Float_t *xyz)
1252 //------------------------------------------------
1253 // ExB at the wires and wire number calulation
1254 //------------------------------------------------
1256 //-----------------------------------------------------------------
1257 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1258 //-----------------------------------------------------------------
1259 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1262 fTPCParam->GetWire(x1); //calculate nearest wire position
1263 Float_t dx=xyz[0]-x1;
1264 xyz[1]+=dx*fTPCParam->GetOmegaTau();
1267 //_____________________________________________________________________________
1268 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet)
1270 //-----------------------------------------------------------
1271 // Single row digitization, coupling from the neighbouring
1272 // rows taken into account
1273 //-----------------------------------------------------------
1275 //-----------------------------------------------------------------
1276 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1277 //-----------------------------------------------------------------
1280 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1281 Float_t chipgain= fTPCParam->GetChipGain();
1282 Float_t zerosup = fTPCParam->GetZeroSup();
1283 Int_t nrows =fTPCParam->GetNRow(isec);
1287 Int_t IndexRange[4];
1292 // iFlag = 0 -> inner row, iFlag = 1 -> the middle one, iFlag = 2 -> the outer one
1295 nTracks[0]=rowTriplet[0]->GetEntries(); // lower row
1296 nTracks[1]=rowTriplet[1]->GetEntries(); // middle row
1297 nTracks[2]=rowTriplet[2]->GetEntries(); // upper row
1302 n_of_pads[0]=fTPCParam->GetNPads(isec,0);
1303 n_of_pads[1]=fTPCParam->GetNPads(isec,1);
1305 else if(irow == nrows-1){
1307 n_of_pads[1]=fTPCParam->GetNPads(isec,irow-1);
1308 n_of_pads[2]=fTPCParam->GetNPads(isec,irow);
1312 for(i1=0;i1<3;i1++){
1313 n_of_pads[i1]=fTPCParam->GetNPads(isec,irow-1+i1);
1318 // Integrated signal for this row
1319 // and a single track signal
1322 TMatrix *m1 = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // integrated
1323 TMatrix *m2 = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // single
1327 TMatrix &Total = *m1;
1329 // Array of pointers to the label-signal list
1331 Int_t NofDigits = n_of_pads[iFlag]*MAXTPCTBK; // number of digits for this row
1333 Float_t **pList = new Float_t* [NofDigits];
1337 for(lp=0;lp<NofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1340 // Straight signal and cross-talk, cross-talk is integrated over all
1341 // tracks and added to the signal at the very end
1345 for (i1=0;i1<nTracks[iFlag];i1++){
1347 m2->Zero(); // clear single track signal matrix
1349 Float_t TrackLabel =
1350 GetSignal(rowTriplet[iFlag],i1,n_of_pads[iFlag],m2,m1,IndexRange);
1352 GetList(TrackLabel,n_of_pads[iFlag],m2,IndexRange,pList);
1357 // Cross talk from the neighbouring pad-rows
1360 TMatrix *m3 = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // cross-talk
1362 TMatrix &Cross = *m3;
1366 // cross-talk from the outer row only (first pad row)
1368 GetCrossTalk(0,rowTriplet[1],nTracks[1],n_of_pads,m3);
1371 else if(iFlag == 2){
1373 // cross-talk from the inner row only (last pad row)
1375 GetCrossTalk(2,rowTriplet[1],nTracks[1],n_of_pads,m3);
1380 // cross-talk from both inner and outer rows
1382 GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner
1383 GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer
1386 Total += Cross; // add the cross-talk
1389 // Convert analog signal to ADC counts
1396 for(Int_t ip=1;ip<n_of_pads[iFlag]+1;ip++){
1397 for(Int_t it=1;it<MAXTPCTBK+1;it++){
1399 Float_t q = Total(ip,it);
1401 Int_t gi =(it-1)*n_of_pads[iFlag]+ip-1; // global index
1403 q = gRandom->Gaus(q,fTPCParam->GetNoise()); // apply noise
1404 q *= (q_el*1.e15); // convert to fC
1405 q *= chipgain; // convert to mV
1406 q *= (adc_sat/dyn_range); // convert to ADC counts
1408 if(q <zerosup) continue; // do not fill zeros
1409 if(q > adc_sat) q = adc_sat; // saturation
1412 // "real" signal or electronic noise (list = -1)?
1415 for(Int_t j1=0;j1<3;j1++){
1416 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1423 digits[4]= (Int_t)q;
1427 AddDigit(tracks,digits);
1429 } // end of loop over time buckets
1430 } // end of lop over pads
1433 // This row has been digitized, delete nonused stuff
1436 for(lp=0;lp<NofDigits;lp++){
1437 if(pList[lp]) delete [] pList[lp];
1446 } // end of DigitizeRow
1447 //_____________________________________________________________________________
1448 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, Int_t np, TMatrix *m1, TMatrix *m2,
1452 //---------------------------------------------------------------
1453 // Calculates 2-D signal (pad,time) for a single track,
1454 // returns a pointer to the signal matrix and the track label
1455 // No digitization is performed at this level!!!
1456 //---------------------------------------------------------------
1458 //-----------------------------------------------------------------
1459 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1460 //-----------------------------------------------------------------
1463 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1464 AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
1465 AliTPCRF1D * fRF = &(fDigParam->GetRF());
1467 //to make the code faster we put parameters to the stack
1469 Float_t zwidth = fTPCParam->GetZWidth();
1470 Float_t zwidthm1 =1./fTPCParam->GetZWidth();
1472 tv = (TVector*)p1->At(ntr); // pointer to a track
1475 Float_t label = v(0);
1477 Int_t CentralPad = (np+1)/2;
1479 Int_t nElectrons = (tv->GetNrows()-1)/4;
1480 Float_t range=((np-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); // pad range
1481 range -= 0.5; // dead zone, 5mm from the edge, according to H.G. Fischer
1483 Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1486 Float_t PadSignal[7]; // signal from a single electron
1488 TMatrix &signal = *m1;
1489 TMatrix &total = *m2;
1492 IndexRange[0]=9999; // min pad
1493 IndexRange[1]=-1; // max pad
1494 IndexRange[2]=9999; //min time
1495 IndexRange[3]=-1; // max time
1498 // Loop over all electrons
1501 for(Int_t nel=0; nel<nElectrons; nel++){
1503 Float_t xwire = v(idx+1);
1504 Float_t y = v(idx+2);
1505 Float_t z = v(idx+3);
1508 Float_t absy=TMath::Abs(y);
1510 if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
1511 PadNumber=CentralPad;
1513 else if (absy < range){
1514 PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.);
1515 PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
1517 else continue; // electron out of pad-range , lost at the sector edge
1519 Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
1522 Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
1523 for (Int_t i=0;i<7;i++){
1524 PadSignal[i]=fPRF2D->GetPRF(dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
1525 PadSignal[i] *= fTPCParam->GetPadCoupling();
1528 Int_t LeftPad = TMath::Max(1,PadNumber-3);
1529 Int_t RightPad = TMath::Min(np,PadNumber+3);
1531 Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
1532 Int_t pmax=RightPad-PadNumber+3; // upper index
1534 Float_t z_drift = (z_end-z)*zwidthm1;
1535 Float_t z_offset = z_drift-(Int_t)z_drift;
1536 //distance to the centre of nearest time bin (in time bin units)
1537 Int_t FirstBucket = (Int_t)z_drift+1;
1540 // loop over time bins (4 bins is enough - 3 sigma truncated Gaussian)
1541 for (Int_t i2=0;i2<4;i2++){
1542 Int_t TrueTime = FirstBucket+i2; // current time bucket
1543 Float_t dz = (Float_t(i2)+z_offset)*zwidth;
1544 Float_t ampl = fRF->GetRF(dz);
1545 if( (TrueTime>MAXTPCTBK) ) break; // beyond the time range
1547 IndexRange[2]=TMath::Min(IndexRange[2],TrueTime); // min time
1548 IndexRange[3]=TMath::Max(IndexRange[3],TrueTime); // max time
1550 // loop over pads, from pmin to pmax
1551 for(Int_t i3=pmin;i3<=pmax;i3++){
1552 Int_t TruePad = LeftPad+i3-pmin;
1553 IndexRange[0]=TMath::Min(IndexRange[0],TruePad); // min pad
1554 IndexRange[1]=TMath::Max(IndexRange[1],TruePad); // max pad
1555 signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!!
1556 total(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!!
1557 } // end of pads loop
1558 } // end of time loop
1559 } // end of loop over electrons
1561 return label; // returns track label when finished
1564 //_____________________________________________________________________________
1565 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange,
1568 //----------------------------------------------------------------------
1569 // Updates the list of tracks contributing to digits for a given row
1570 //----------------------------------------------------------------------
1572 //-----------------------------------------------------------------
1573 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1574 //-----------------------------------------------------------------
1576 TMatrix &signal = *m;
1578 // lop over nonzero digits
1580 for(Int_t it=IndexRange[2];it<IndexRange[3]+1;it++){
1581 for(Int_t ip=IndexRange[0];ip<IndexRange[1]+1;ip++){
1584 Int_t GlobalIndex = (it-1)*np+ip-1; // GlobalIndex starts from 0!
1586 if(!pList[GlobalIndex]){
1589 // Create new list (6 elements - 3 signals and 3 labels),
1590 // but only if the signal is > 0.
1593 if(signal(ip,it)>0.){
1595 pList[GlobalIndex] = new Float_t [6];
1599 *pList[GlobalIndex] = -1.;
1600 *(pList[GlobalIndex]+1) = -1.;
1601 *(pList[GlobalIndex]+2) = -1.;
1602 *(pList[GlobalIndex]+3) = -1.;
1603 *(pList[GlobalIndex]+4) = -1.;
1604 *(pList[GlobalIndex]+5) = -1.;
1607 *pList[GlobalIndex] = label;
1608 *(pList[GlobalIndex]+3) = signal(ip,it);}
1612 // check the signal magnitude
1614 Float_t highest = *(pList[GlobalIndex]+3);
1615 Float_t middle = *(pList[GlobalIndex]+4);
1616 Float_t lowest = *(pList[GlobalIndex]+5);
1619 // compare the new signal with already existing list
1622 if(signal(ip,it)<lowest) continue; // neglect this track
1626 if (signal(ip,it)>highest){
1627 *(pList[GlobalIndex]+5) = middle;
1628 *(pList[GlobalIndex]+4) = highest;
1629 *(pList[GlobalIndex]+3) = signal(ip,it);
1631 *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1632 *(pList[GlobalIndex]+1) = *pList[GlobalIndex];
1633 *pList[GlobalIndex] = label;
1635 else if (signal(ip,it)>middle){
1636 *(pList[GlobalIndex]+5) = middle;
1637 *(pList[GlobalIndex]+4) = signal(ip,it);
1639 *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1640 *(pList[GlobalIndex]+1) = label;
1643 *(pList[GlobalIndex]+5) = signal(ip,it);
1644 *(pList[GlobalIndex]+2) = label;
1648 } // end of loop over pads
1649 } // end of loop over time bins
1655 //___________________________________________________________________
1656 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1657 Stat_t ntracks,TObjArray **row)
1660 //-----------------------------------------------------------------
1661 // Prepares the sector digitization, creates the vectors of
1662 // tracks for each row of this sector. The track vector
1663 // contains the track label and the position of electrons.
1664 //-----------------------------------------------------------------
1666 //-----------------------------------------------------------------
1667 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1668 //-----------------------------------------------------------------
1670 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1674 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1676 //----------------------------------------------
1677 // Create TObjArray-s, one for each row,
1678 // each TObjArray will store the TVectors
1679 // of electrons, one TVector per each track.
1680 //----------------------------------------------
1682 for(i=0; i<nrows; i++){
1683 row[i] = new TObjArray;
1685 Int_t *n_of_electrons = new Int_t [nrows]; // electron counter for each row
1686 TVector **tr = new TVector* [nrows]; //pointers to the track vectors
1688 //--------------------------------------------------------------------
1689 // Loop over tracks, the "track" contains the full history
1690 //--------------------------------------------------------------------
1692 Int_t previousTrack,currentTrack;
1693 previousTrack = -1; // nothing to store so far!
1695 for(Int_t track=0;track<ntracks;track++){
1699 TH->GetEvent(track); // get next track
1700 Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1702 if(nhits == 0) continue; // no hits in the TPC for this track
1704 //--------------------------------------------------------------
1706 //--------------------------------------------------------------
1708 for(Int_t hit=0;hit<nhits;hit++){
1710 tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1712 Int_t sector=tpcHit->fSector; // sector number
1713 if(sector != isec) continue; //terminate iteration
1715 currentTrack = tpcHit->fTrack; // track number
1716 if(currentTrack != previousTrack){
1718 // store already filled fTrack
1720 for(i=0;i<nrows;i++){
1721 if(previousTrack != -1){
1722 if(n_of_electrons[i]>0){
1723 TVector &v = *tr[i];
1724 v(0) = previousTrack;
1725 tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary
1729 delete tr[i]; // delete empty TVector
1734 n_of_electrons[i]=0;
1735 tr[i] = new TVector(361); // TVectors for the next fTrack
1737 } // end of loop over rows
1739 previousTrack=currentTrack; // update track label
1742 Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1744 //---------------------------------------------------
1745 // Calculate the electron attachment probability
1746 //---------------------------------------------------
1748 Float_t time = 1.e6*(z_end-TMath::Abs(tpcHit->fZ))/fTPCParam->GetDriftV();
1750 Float_t AttProb = fTPCParam->GetAttCoef()*
1751 fTPCParam->GetOxyCont()*time; // fraction!
1753 //-----------------------------------------------
1754 // Loop over electrons
1755 //-----------------------------------------------
1757 for(Int_t nel=0;nel<QI;nel++){
1758 // skip if electron lost due to the attachment
1759 if((gRandom->Rndm(0)) < AttProb) continue; // electron lost!
1763 ElDiff(xyz); // Appply the diffusion
1765 fTPCParam->XYZtoCRXYZ(xyz,isec,row_number,3);
1767 //transform position to local coordinates
1768 //option 3 means that we calculate x position relative to
1771 if ((row_number<0)||row_number>=fTPCParam->GetNRow(isec)) continue;
1772 ExB(xyz); // ExB effect at the sense wires
1773 n_of_electrons[row_number]++;
1774 //----------------------------------
1775 // Expand vector if necessary
1776 //----------------------------------
1777 if(n_of_electrons[row_number]>120){
1778 Int_t range = tr[row_number]->GetNrows();
1779 if(n_of_electrons[row_number] > (range-1)/3){
1780 tr[row_number]->ResizeTo(range+150); // Add 50 electrons
1784 TVector &v = *tr[row_number];
1785 Int_t idx = 3*n_of_electrons[row_number]-2;
1787 v(idx)= xyz[0]; // X
1788 v(idx+1)=xyz[1]; // Y (along the pad-row)
1789 v(idx+2)=xyz[2]; // Z
1791 } // end of loop over electrons
1793 } // end of loop over hits
1794 } // end of loop over tracks
1797 // store remaining track (the last one) if not empty
1800 for(i=0;i<nrows;i++){
1801 if(n_of_electrons[i]>0){
1802 TVector &v = *tr[i];
1803 v(0) = previousTrack;
1804 tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary
1814 delete [] n_of_electrons;
1816 } // end of MakeSector
1817 //_____________________________________________________________________________
1818 void AliTPC::GetCrossTalk (Int_t iFlag,TObjArray *p,Int_t ntracks,Int_t *npads,
1822 //-------------------------------------------------------------
1823 // Calculates the cross-talk from one row (inner or outer)
1824 //-------------------------------------------------------------
1826 //-----------------------------------------------------------------
1827 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1828 //-----------------------------------------------------------------
1831 // iFlag=2 & 3 --> cross-talk from the inner row
1832 // iFlag=0 & 4 --> cross-talk from the outer row
1834 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1835 AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
1836 AliTPCRF1D * fRF = &(fDigParam->GetRF());
1838 //to make code faster
1840 Float_t zwidth = fTPCParam->GetZWidth();
1841 Float_t zwidthm1 =1/fTPCParam->GetZWidth();
1846 Int_t nPadsSignal; // for this pads the signal is calculated
1847 Float_t range; // sense wire range
1850 Float_t IneffFactor=0.5; // gain inefficiency close to the edges
1854 nPadsSignal = npads[1];
1855 range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
1856 nPadsDiff = (npads[1]-npads[0])/2;
1858 else if (iFlag == 2){
1860 nPadsSignal = npads[2];
1861 range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
1864 else if (iFlag == 3){
1866 nPadsSignal = npads[1];
1867 range = ((npads[0]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
1872 nPadsSignal = npads[2];
1873 range = ((npads[2]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
1874 nPadsDiff = (npads[2]-npads[1])/2;
1877 range-=0.5; // dead zone close to the edges
1880 TMatrix &signal = *m;
1882 Int_t CentralPad = (nPadsSignal+1)/2;
1883 Float_t PadSignal[7]; // signal from a single electron
1885 for(Int_t nt=0;nt<ntracks;nt++){
1886 tv=(TVector*)p->At(nt); // pointer to a track
1888 Int_t nElectrons = (tv->GetNrows()-1)/4;
1889 // Loop over electrons
1890 for(Int_t nel=0; nel<nElectrons; nel++){
1894 if (iFlag==0) xwire+=fTPCParam->GetPadPitchLength();
1895 if (iFlag==2) xwire-=fTPCParam->GetPadPitchLength();
1896 if (iFlag==3) xwire-=fTPCParam->GetPadPitchLength();
1897 if (iFlag==4) xwire+=fTPCParam->GetPadPitchLength();
1899 // electron acceptance for the cross-talk (only the closest wire)
1901 Float_t dxMax = fTPCParam->GetPadPitchLength()*0.5+fTPCParam->GetWWPitch();
1902 if(TMath::Abs(xwire)>dxMax) continue;
1904 Float_t y = v(idx+2);
1905 Float_t z = v(idx+3);
1906 Float_t absy=TMath::Abs(y);
1908 if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
1909 PadNumber=CentralPad;
1911 else if (absy < range){
1912 PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.);
1913 PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
1915 else continue; // electron out of sense wire range, lost at the sector edge
1917 Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
1919 Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
1921 for (Int_t i=0;i<7;i++){
1922 PadSignal[i]=fPRF2D->GetPRF(dist+(3-i)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
1924 PadSignal[i] *= fTPCParam->GetPadCoupling();
1928 Int_t LeftPad = TMath::Max(1,PadNumber-3);
1929 Int_t RightPad = TMath::Min(nPadsSignal,PadNumber+3);
1931 Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
1932 Int_t pmax=RightPad-PadNumber+3; // upper index
1935 Float_t z_drift = (z_end-z)*zwidthm1;
1936 Float_t z_offset = z_drift-(Int_t)z_drift;
1937 //distance to the centre of nearest time bin (in time bin units)
1938 Int_t FirstBucket = (Int_t)z_drift+1;
1939 // MI check it --time offset
1940 for (Int_t i2=0;i2<4;i2++){
1941 Int_t TrueTime = FirstBucket+i2; // current time bucket
1942 Float_t dz = (Float_t(i2)+z_offset)*zwidth;
1943 Float_t ampl = fRF->GetRF(dz);
1944 if((TrueTime>MAXTPCTBK)) break; // beyond the time range
1947 // loop over pads, from pmin to pmax
1949 for(Int_t i3=pmin;i3<pmax+1;i3++){
1950 Int_t TruePad = LeftPad+i3-pmin;
1952 if(TruePad<nPadsDiff+1 || TruePad > nPadsSignal-nPadsDiff) continue;
1954 TruePad -= nPadsDiff;
1955 signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!
1957 } // end of loop over pads
1958 } // end of loop over time bins
1960 } // end of loop over electrons
1962 } // end of loop over tracks
1964 } // end of GetCrossTalk
1968 //_____________________________________________________________________________
1972 // Initialise TPC detector after definition of geometry
1977 for(i=0;i<35;i++) printf("*");
1978 printf(" TPC_INIT ");
1979 for(i=0;i<35;i++) printf("*");
1982 for(i=0;i<80;i++) printf("*");
1986 //_____________________________________________________________________________
1987 void AliTPC::MakeBranch(Option_t* option)
1990 // Create Tree branches for the TPC.
1992 Int_t buffersize = 4000;
1993 char branchname[10];
1994 sprintf(branchname,"%s",GetName());
1996 AliDetector::MakeBranch(option);
1998 char *D = strstr(option,"D");
2000 if (fDigits && gAlice->TreeD() && D) {
2001 gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
2002 printf("Making Branch %s for digits\n",branchname);
2005 char *R = strstr(option,"R");
2007 if (fClusters && gAlice->TreeR() && R) {
2008 gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
2009 printf("Making Branch %s for Clusters\n",branchname);
2013 //_____________________________________________________________________________
2014 void AliTPC::ResetDigits()
2017 // Reset number of digits and the digits array for this detector
2021 // if (fDigits) fDigits->Clear();
2022 if (fDigParam->GetArray()!=0) fDigParam->GetArray()->Clear();
2024 if (fClusters) fClusters->Clear();
2027 //_____________________________________________________________________________
2028 void AliTPC::SetSecAL(Int_t sec)
2030 //---------------------------------------------------
2031 // Activate/deactivate selection for lower sectors
2032 //---------------------------------------------------
2034 //-----------------------------------------------------------------
2035 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2036 //-----------------------------------------------------------------
2041 //_____________________________________________________________________________
2042 void AliTPC::SetSecAU(Int_t sec)
2044 //----------------------------------------------------
2045 // Activate/deactivate selection for upper sectors
2046 //---------------------------------------------------
2048 //-----------------------------------------------------------------
2049 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2050 //-----------------------------------------------------------------
2055 //_____________________________________________________________________________
2056 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2058 //----------------------------------------
2059 // Select active lower sectors
2060 //----------------------------------------
2062 //-----------------------------------------------------------------
2063 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2064 //-----------------------------------------------------------------
2074 //_____________________________________________________________________________
2075 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2076 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2077 Int_t s11 , Int_t s12)
2079 //--------------------------------
2080 // Select active upper sectors
2081 //--------------------------------
2083 //-----------------------------------------------------------------
2084 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2085 //-----------------------------------------------------------------
2101 //_____________________________________________________________________________
2102 void AliTPC::SetSens(Int_t sens)
2105 //-------------------------------------------------------------
2106 // Activates/deactivates the sensitive strips at the center of
2107 // the pad row -- this is for the space-point resolution calculations
2108 //-------------------------------------------------------------
2110 //-----------------------------------------------------------------
2111 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2112 //-----------------------------------------------------------------
2117 //_____________________________________________________________________________
2118 void AliTPC::Streamer(TBuffer &R__b)
2121 // Stream an object of class AliTPC.
2123 if (R__b.IsReading()) {
2124 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2125 AliDetector::Streamer(R__b);
2126 if (R__v < 2) return;
2130 fClustersIndex = new Int_t[fNsectors+1];
2131 fDigitsIndex = new Int_t[fNsectors+1];
2133 R__b.WriteVersion(AliTPC::IsA());
2134 AliDetector::Streamer(R__b);
2141 ClassImp(AliTPCcluster)
2143 //_____________________________________________________________________________
2144 AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
2147 // Creates a simulated cluster for the TPC
2149 fTracks[0] = lab[0];
2150 fTracks[1] = lab[1];
2151 fTracks[2] = lab[2];
2161 //_____________________________________________________________________________
2162 void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const
2165 // Transformation from local to global coordinate system
2167 x[0]=par->GetPadRowRadii(fSector,fPadRow-1);
2170 par->CRXYZtoXYZ(x,fSector,fPadRow-1,1);
2173 //_____________________________________________________________________________
2174 Int_t AliTPCcluster::Compare(TObject * o)
2177 // compare two clusters according y coordinata
2178 AliTPCcluster *cl= (AliTPCcluster *)o;
2179 if (fY<cl->fY) return -1;
2180 if (fY==cl->fY) return 0;
2184 Bool_t AliTPCcluster::IsSortable() const
2187 //make AliTPCcluster sortabale
2193 ClassImp(AliTPCdigit)
2195 //_____________________________________________________________________________
2196 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2200 // Creates a TPC digit object
2202 fSector = digits[0];
2203 fPadRow = digits[1];
2206 fSignal = digits[4];
2212 //_____________________________________________________________________________
2213 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2217 // Creates a TPC hit object
2228 ClassImp(AliTPCtrack)
2230 //_____________________________________________________________________________
2231 AliTPCtrack::AliTPCtrack(Float_t *hits)
2234 // Default creator for a TPC reconstructed track object
2236 ref=hits[0]; // This is dummy code !
2239 AliTPCtrack::AliTPCtrack(const AliTPCcluster& c,const TVector& xx,
2240 const TMatrix& CC, const AliTPCParam *p):
2241 x(xx),C(CC),clusters(MAX_CLUSTER)
2244 // Standard creator for a TPC reconstructed track object
2247 int sec=c.fSector-1, row=c.fPadRow-1;
2248 ref=p->GetPadRowRadii(sec+1,row);
2251 fAlpha=(sec%12)*alpha_low;
2253 fAlpha=(sec%24)*alpha_up;
2255 clusters.AddLast((AliTPCcluster*)(&c));
2258 //_____________________________________________________________________________
2259 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
2260 clusters(t.clusters.GetEntriesFast())
2263 // Copy creator for a TPC reconstructed track
2268 int n=t.clusters.GetEntriesFast();
2269 for (int i=0; i<n; i++) clusters.AddLast(t.clusters.UncheckedAt(i));
2272 //_____________________________________________________________________________
2273 Double_t AliTPCtrack::GetY(Double_t xk) const
2278 Double_t c2=x(2)*xk - x(3);
2279 if (TMath::Abs(c2) >= 0.999) {
2280 if (*this>10) cerr<<*this<<" AliTPCtrack warning: No y for this x !\n";
2283 Double_t c1=x(2)*ref - x(3);
2284 Double_t r1=sqrt(1.-c1*c1), r2=sqrt(1.-c2*c2);
2286 return x(0) + dx*(c1+c2)/(r1+r2);
2289 //_____________________________________________________________________________
2290 int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
2293 // Propagate a TPC reconstructed track
2295 if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
2296 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
2300 Double_t x1=ref, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
2301 Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
2302 Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
2304 x(0) += dx*(c1+c2)/(r1+r2);
2305 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2307 TMatrix F(5,5); F.UnitMatrix();
2308 Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
2309 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2310 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2311 Double_t cr=c1*r2+c2*r1;
2312 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2313 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2315 TMatrix tmp(F,TMatrix::kMult,C);
2316 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2320 //Multiple scattering******************
2321 Double_t ey=x(2)*ref - x(3);
2322 Double_t ex=sqrt(1-ey*ey);
2324 TMatrix Q(5,5); Q=0.;
2325 Q(2,2)=ez*ez+ey*ey; Q(2,3)=-ex*ey; Q(2,4)=-ex*ez;
2326 Q(3,2)=Q(2,3); Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
2327 Q(4,2)=Q(2,4); Q(4,3)= Q(3,4); Q(4,4)=1.;
2330 F(2,2)=-x(2)*ex; F(2,3)=-x(2)*ey;
2331 F(3,2)=-ex*(x(2)*ref-ey); F(3,3)=-(1.+ x(2)*ref*ey - ey*ey);
2332 F(4,2)=-ez*ex; F(4,3)=-ez*ey; F(4,4)=1.;
2335 Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2337 Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
2338 Double_t beta2=p2/(p2 + pm*pm);
2339 Double_t d=sqrt((x1-ref)*(x1-ref)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
2341 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
2345 //Energy losses************************
2346 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
2347 if (x1 < x2) dE=-dE;
2348 x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2349 //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2351 x1=ref; x2=xk; y1=x(0); z1=x(1);
2352 c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1);
2353 c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2);
2355 x(0) += dx*(c1+c2)/(r1+r2);
2356 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2359 rr=r1+r2; cc=c1+c2; xx=x1+x2;
2360 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2361 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2363 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2364 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2367 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2374 //_____________________________________________________________________________
2375 void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
2378 // Propagate a reconstructed track from the vertex
2380 Double_t c=x(2)*ref - x(3);
2381 Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
2382 Double_t snf=tgf/sqrt(1.+ tgf*tgf);
2383 Double_t xv=(x(3)+snf)/x(2);
2384 PropagateTo(xv,x0,rho,pm);
2387 //_____________________________________________________________________________
2388 void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
2391 // Update statistics for a reconstructed TPC track
2393 TMatrix H(2,5); H.UnitMatrix();
2394 TMatrix Ht(TMatrix::kTransposed,H);
2395 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
2396 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2398 TMatrix tmp(H,TMatrix::kMult,C);
2399 TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
2401 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2402 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2403 R(1,0)*=-1; R(0,1)=R(1,0);
2408 TMatrix K(C,TMatrix::kMult,Ht); K*=R;
2411 x*=H; x-=m; x*=-1; x*=K; x+=savex;
2412 if (TMath::Abs(x(2)*ref-x(3)) >= 0.999) {
2413 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
2419 C.Mult(K,tmp); C-=saveC; C*=-1;
2421 clusters.AddLast((AliTPCcluster*)c);
2425 //_____________________________________________________________________________
2426 int AliTPCtrack::Rotate(Double_t alpha)
2429 // Rotate a reconstructed TPC track
2433 Double_t x1=ref, y1=x(0);
2434 Double_t ca=cos(alpha), sa=sin(alpha);
2435 Double_t r1=x(2)*ref - x(3);
2437 ref = x1*ca + y1*sa;
2438 x(0)=-x1*sa + y1*ca;
2439 x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
2441 Double_t r2=x(2)*ref - x(3);
2442 if (TMath::Abs(r2) >= 0.999) {
2443 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
2447 Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
2448 if ((x(0)-y0)*x(2) >= 0.) {
2449 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
2453 TMatrix F(5,5); F.UnitMatrix();
2456 F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa;
2457 F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
2458 TMatrix tmp(F,TMatrix::kMult,C);
2459 // Double_t dy2=C(0,0);
2460 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2461 // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
2462 // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
2467 //_____________________________________________________________________________
2468 void AliTPCtrack::UseClusters() const
2473 int num_of_clusters=clusters.GetEntriesFast();
2474 for (int i=0; i<num_of_clusters; i++) {
2475 //if (i<=14) continue;
2476 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
2481 //_____________________________________________________________________________
2482 Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const
2485 // Calculate chi2 for a reconstructed TPC track
2487 TMatrix H(2,5); H.UnitMatrix();
2488 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
2489 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2490 TVector res=x; res*=H; res-=m; //res*=-1;
2491 TMatrix tmp(H,TMatrix::kMult,C);
2492 TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
2494 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2495 if (TMath::Abs(det) < 1.e-10) {
2496 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
2499 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2500 R(1,0)*=-1; R(0,1)=R(1,0);
2510 //_____________________________________________________________________________
2511 int AliTPCtrack::GetLab() const
2520 } s[MAX_CLUSTER]={{0,0}};
2523 int num_of_clusters=clusters.GetEntriesFast();
2524 for (i=0; i<num_of_clusters; i++) {
2525 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
2526 lab=TMath::Abs(c->fTracks[0]);
2528 for (j=0; j<MAX_CLUSTER; j++)
2529 if (s[j].lab==lab || s[j].max==0) break;
2535 for (i=0; i<num_of_clusters; i++)
2536 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
2538 for (i=0; i<num_of_clusters; i++) {
2539 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
2540 if (TMath::Abs(c->fTracks[1]) == lab ||
2541 TMath::Abs(c->fTracks[2]) == lab ) max++;
2544 if (1.-float(max)/num_of_clusters > 0.10) return -lab;
2546 if (num_of_clusters < 6) return lab;
2549 for (i=1; i<=6; i++) {
2550 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(num_of_clusters-i);
2551 if (lab == TMath::Abs(c->fTracks[0]) ||
2552 lab == TMath::Abs(c->fTracks[1]) ||
2553 lab == TMath::Abs(c->fTracks[2])) max++;
2555 if (max<3) return -lab;
2560 //_____________________________________________________________________________
2561 void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
2564 // Get reconstructed TPC track momentum
2566 Double_t pt=0.3*FIELD/TMath::Abs(x(2))/100; // GeV/c
2567 Double_t r=x(2)*ref-x(3);
2568 Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2569 px=-pt*(x(0)-y0)*x(2); //cos(phi);
2570 py=-pt*(x(3)-ref*x(2)); //sin(phi);
2572 Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2573 py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2577 //_____________________________________________________________________________
2579 // Classes for internal tracking use
2582 //_____________________________________________________________________________
2583 void AliTPCRow::InsertCluster(const AliTPCcluster* c)
2586 // Insert a cluster in the list
2588 if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2589 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2591 if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2593 memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2594 clusters[i]=c; num_of_clusters++;
2597 //_____________________________________________________________________________
2598 int AliTPCRow::Find(Double_t y) const
2603 if (y <= clusters[0]->fY) return 0;
2604 if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2605 int b=0, e=num_of_clusters-1, m=(b+e)/2;
2606 for (; b<e; m=(b+e)/2) {
2607 if (y > clusters[m]->fY) b=m+1;