1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 Revision 1.14 1999/09/29 09:24:33 fca
19 Introduction of the Copyright and cvs Log
23 ///////////////////////////////////////////////////////////////////////////////
25 // Time Projection Chamber //
26 // This class contains the basic functions for the Time Projection Chamber //
27 // detector. Functions specific to one particular geometry are //
28 // contained in the derived classes //
32 <img src="picts/AliTPCClass.gif">
37 ///////////////////////////////////////////////////////////////////////////////
43 #include <TGeometry.h>
46 #include <TObjectTable.h>
47 #include "TParticle.h"
55 #include "AliTPCParam.h"
57 #include "AliTPCPRF2D.h"
58 #include "AliTPCRF1D.h"
64 //_____________________________________________________________________________
68 // Default constructor
79 fDigParam= new AliTPCD();
80 fDigits = fDigParam->GetArray();
83 //_____________________________________________________________________________
84 AliTPC::AliTPC(const char *name, const char *title)
85 : AliDetector(name,title)
88 // Standard constructor
92 // Initialise arrays of hits and digits
93 fHits = new TClonesArray("AliTPChit", 176);
94 // fDigits = new TClonesArray("AliTPCdigit",10000);
96 fDigParam= new AliTPCD;
97 fDigits = fDigParam->GetArray();
99 AliTPCParam *fTPCParam = &(fDigParam->GetParam());
102 // Initialise counters
106 fNsectors = fTPCParam->GetNSector();
109 fDigitsIndex = new Int_t[fNsectors+1];
113 // Initialise color attributes
114 SetMarkerColor(kYellow);
117 //_____________________________________________________________________________
129 if (fDigitsIndex) delete [] fDigitsIndex;
132 //_____________________________________________________________________________
133 void AliTPC::AddCluster(Float_t *hits, Int_t *tracks)
136 // Add a simulated cluster to the list
138 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
139 TClonesArray &lclusters = *fClusters;
140 new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks);
143 //_____________________________________________________________________________
144 void AliTPC::AddCluster(const AliTPCcluster &c)
147 // Add a simulated cluster copy to the list
149 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
150 TClonesArray &lclusters = *fClusters;
151 new(lclusters[fNclusters++]) AliTPCcluster(c);
154 //_____________________________________________________________________________
155 void AliTPC::AddDigit(Int_t *tracks, Int_t *digits)
158 // Add a TPC digit to the list
160 // TClonesArray &ldigits = *fDigits;
162 TClonesArray &ldigits = *fDigParam->GetArray();
163 new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits);
166 //_____________________________________________________________________________
167 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
170 // Add a hit to the list
172 TClonesArray &lhits = *fHits;
173 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
176 //_____________________________________________________________________________
177 void AliTPC::AddTrack(Float_t *hits)
180 // Add a track to the list of tracks
182 TClonesArray <racks = *fTracks;
183 new(ltracks[fNtracks++]) AliTPCtrack(hits);
186 //_____________________________________________________________________________
187 void AliTPC::AddTrack(const AliTPCtrack& t)
190 // Add a track copy to the list of tracks
192 if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000);
193 TClonesArray <racks = *fTracks;
194 new(ltracks[fNtracks++]) AliTPCtrack(t);
197 //_____________________________________________________________________________
198 void AliTPC::BuildGeometry()
201 // Build TPC ROOT TNode geometry for the event display
206 const int kColorTPC=19;
207 char name[5], title[25];
208 const Double_t kDegrad=TMath::Pi()/180;
209 const Double_t kRaddeg=180./TMath::Pi();
211 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
213 Float_t InnerOpenAngle = fTPCParam->GetInnerAngle();
214 Float_t OuterOpenAngle = fTPCParam->GetOuterAngle();
216 Float_t InnerAngleShift = fTPCParam->GetInnerAngleShift();
217 Float_t OuterAngleShift = fTPCParam->GetOuterAngleShift();
219 Int_t nLo = fTPCParam->GetNInnerSector()/2;
220 Int_t nHi = fTPCParam->GetNOuterSector()/2;
222 const Double_t loAng = (Double_t)TMath::Nint(InnerOpenAngle*kRaddeg);
223 const Double_t hiAng = (Double_t)TMath::Nint(OuterOpenAngle*kRaddeg);
224 const Double_t loAngSh = (Double_t)TMath::Nint(InnerAngleShift*kRaddeg);
225 const Double_t hiAngSh = (Double_t)TMath::Nint(OuterAngleShift*kRaddeg);
228 const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad);
229 const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad);
235 // Get ALICE top node
238 Top=gAlice->GetGeometry()->GetNode("alice");
242 rl = fTPCParam->GetInSecLowEdge();
243 ru = fTPCParam->GetInSecUpEdge();
247 sprintf(name,"LS%2.2d",i);
249 sprintf(title,"TPC low sector %3d",i);
252 tubs = new TTUBS(name,title,"void",rl*loCorr,ru*loCorr,250.,
253 loAng*(i-0.5)+loAngSh,loAng*(i+0.5)+loAngSh);
254 tubs->SetNumberOfDivisions(1);
256 Node = new TNode(name,title,name,0,0,0,"");
257 Node->SetLineColor(kColorTPC);
263 rl = fTPCParam->GetOuSecLowEdge();
264 ru = fTPCParam->GetOuSecUpEdge();
267 sprintf(name,"US%2.2d",i);
269 sprintf(title,"TPC upper sector %d",i);
271 tubs = new TTUBS(name,title,"void",rl*hiCorr,ru*hiCorr,250,
272 hiAng*(i-0.5)+hiAngSh,hiAng*(i+0.5)+hiAngSh);
273 tubs->SetNumberOfDivisions(1);
275 Node = new TNode(name,title,name,0,0,0,"");
276 Node->SetLineColor(kColorTPC);
283 //_____________________________________________________________________________
284 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
287 // Calculate distance from TPC to mouse on the display
293 //_____________________________________________________________________________
294 static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
297 // Parametrised error of the cluster reconstruction (pad direction)
299 pt=TMath::Abs(pt)*1000.;
302 Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x;
303 if (s<0.4e-3) s=0.4e-3;
304 s*=1.3; //Iouri Belikov
308 //_____________________________________________________________________________
309 static Double_t SigmaZ2(Double_t r, Double_t tgl)
312 // Parametrised error of the cluster reconstruction (drift direction)
315 Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl;
316 if (s<0.4e-3) s=0.4e-3;
317 s*=1.3; //Iouri Belikov
321 //_____________________________________________________________________________
322 inline Double_t f1(Double_t x1,Double_t y1,
323 Double_t x2,Double_t y2,
324 Double_t x3,Double_t y3)
326 //-----------------------------------------------------------------
327 // Initial approximation of the track curvature
329 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
330 //-----------------------------------------------------------------
331 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
332 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
333 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
334 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
335 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
337 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
339 return -xr*yr/sqrt(xr*xr+yr*yr);
343 //_____________________________________________________________________________
344 inline Double_t f2(Double_t x1,Double_t y1,
345 Double_t x2,Double_t y2,
346 Double_t x3,Double_t y3)
348 //-----------------------------------------------------------------
349 // Initial approximation of the track curvature times center of curvature
351 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
352 //-----------------------------------------------------------------
353 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
354 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
355 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
356 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
357 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
359 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
361 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
364 //_____________________________________________________________________________
365 inline Double_t f3(Double_t x1,Double_t y1,
366 Double_t x2,Double_t y2,
367 Double_t z1,Double_t z2)
369 //-----------------------------------------------------------------
370 // Initial approximation of the tangent of the track dip angle
372 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
373 //-----------------------------------------------------------------
374 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
377 //_____________________________________________________________________________
378 static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec,
381 //-----------------------------------------------------------------
382 // This function tries to find a track prolongation.
384 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
385 //-----------------------------------------------------------------
386 const int ROWS_TO_SKIP=int(0.5*sec->GetNRows());
387 const Float_t MAX_CHI2=12.;
388 int try_again=ROWS_TO_SKIP;
389 Double_t alpha=sec->GetAlpha();
390 int ns=int(2*TMath::Pi()/alpha+0.5);
392 for (int nr=sec->GetRowNumber(t.GetX())-1; nr>=rf; nr--) {
393 Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
394 if (!t.PropagateTo(x)) return 0;
397 Double_t max_chi2=MAX_CHI2;
398 const AliTPCRow& row=sec[s][nr];
399 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
400 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
401 Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
404 if (t>4) cerr<<t<<" FindProlongation warning: Too broad road !\n";
409 for (int i=row.Find(y-road); i<row; i++) {
410 AliTPCcluster* c=(AliTPCcluster*)(row[i]);
411 if (c->fY > y+road) break;
412 if (c->IsUsed()) continue;
413 if ((c->fZ - z)*(c->fZ - z) > 25.*(t.GetSigmaZ2() + sz2)) continue;
414 Double_t chi2=t.GetPredictedChi2(c);
415 if (chi2 > max_chi2) continue;
421 t.Update(cl,max_chi2);
422 Double_t ll=TMath::Sqrt((1+t.GetTgl()*t.GetTgl())/
423 (1-(t.GetC()*x-t.GetEta())*(t.GetC()*x-t.GetEta())));
424 cl->fdEdX = cl->fQ/ll;
425 try_again=ROWS_TO_SKIP;
427 if (try_again==0) break;
430 if (!t.Rotate(alpha)) return 0;
431 } else if (y <-ymax) {
433 if (!t.Rotate(-alpha)) return 0;
443 //_____________________________________________________________________________
444 static void MakeSeeds(TObjArray& seeds,const AliTPCSector *sec, int max_sec,
447 //-----------------------------------------------------------------
448 // This function creates track seeds.
450 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
451 //-----------------------------------------------------------------
452 TMatrix C(5,5); TVector x(5);
453 double alpha=sec->GetAlpha(), shift=sec->GetAlphaShift();
454 double cs=cos(alpha), sn=sin(alpha);
455 for (int ns=0; ns<max_sec; ns++) {
456 int nl=sec[(ns-1+max_sec)%max_sec][i2];
458 int nu=sec[(ns+1)%max_sec][i2];
459 const AliTPCRow& r1=sec[ns][i1];
460 for (int is=0; is < r1; is++) {
461 double x1=sec->GetX(i1), y1=r1[is]->fY, z1=r1[is]->fZ;
462 for (int js=0; js < nl+nm+nu; js++) {
463 const AliTPCcluster *cl;
465 double x2=sec->GetX(i2), y2, z2, tmp;
468 ks=(ns-1+max_sec)%max_sec;
469 const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2];
471 y2=cl->fY; z2=cl->fZ;
473 y2 =-x2*sn+y2*cs; x2=tmp;
477 const AliTPCRow& r2=sec[ns][i2];
479 y2=cl->fY; z2=cl->fZ;
482 const AliTPCRow& r2=sec[(ns+1)%max_sec][i2];
484 y2=cl->fY; z2=cl->fZ;
486 y2 =x2*sn+y2*cs; x2=tmp;
489 double d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
490 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
494 x(2)=f1(x1,y1,x2,y2,0.,0.);
495 x(3)=f2(x1,y1,x2,y2,0.,0.);
496 x(4)=f3(x1,y1,x2,y2,z1,z2);
498 if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue;
500 if (TMath::Abs(x(4)) > 1.2) continue;
502 Double_t a=asin(x(3));
503 Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3)));
504 if (TMath::Abs(zv)>33.) continue;
506 TMatrix X(6,6); X=0.;
507 X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2;
508 X(2,2)=cl->fSigmaY2; X(3,3)=cl->fSigmaZ2;
509 X(4,4)=3./12.; X(5,5)=3./12.;
510 TMatrix F(5,6); F.UnitMatrix();
511 Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1));
512 F(2,0)=(f1(x1,y1+sy,x2,y2,0.,0.)-x(2))/sy;
513 F(2,2)=(f1(x1,y1,x2,y2+sy,0.,0.)-x(2))/sy;
514 F(2,4)=(f1(x1,y1,x2,y2,0.,0.+sy)-x(2))/sy;
515 F(3,0)=(f2(x1,y1+sy,x2,y2,0.,0.)-x(3))/sy;
516 F(3,2)=(f2(x1,y1,x2,y2+sy,0.,0.)-x(3))/sy;
517 F(3,4)=(f2(x1,y1,x2,y2,0.,0.+sy)-x(3))/sy;
518 F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy;
519 F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz;
520 F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy;
521 F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz;
525 TMatrix t(F,TMatrix::kMult,X);
526 C.Mult(t,TMatrix(TMatrix::kTransposed,F));
528 AliTPCtrack *track=new AliTPCtrack(r1[is], x, C, x1, ns*alpha+shift);
529 int rc=FindProlongation(*track,sec,ns,i2);
530 if (rc<0 || *track<(i1-i2)/2) delete track;
531 else seeds.AddLast(track);
537 //_____________________________________________________________________________
538 AliTPCParam *AliTPCSector::param;
539 void AliTPC::Clusters2Tracks()
541 //-----------------------------------------------------------------
542 // This is a track finder.
544 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
545 //-----------------------------------------------------------------
546 if (!fClusters) return;
548 AliTPCParam *p=&fDigParam->GetParam();
549 AliTPCSector::SetParam(p);
551 const int nis=p->GetNInnerSector()/2;
552 AliTPCSSector *ssec=new AliTPCSSector[nis];
553 int nrow_low=ssec->GetNRows();
555 const int nos=p->GetNOuterSector()/2;
556 AliTPCLSector *lsec=new AliTPCLSector[nos];
557 int nrow_up=lsec->GetNRows();
559 int ncl=fClusters->GetEntriesFast();
561 AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl);
562 Int_t sec=c->fSector, row=c->fPadRow;
564 ssec[sec%nis][row].InsertCluster(c);
567 lsec[sec%nos][row].InsertCluster(c);
571 TObjArray seeds(20000);
573 int nrows=nrow_low+nrow_up;
574 int gap=int(0.125*nrows), shift=int(0.5*gap);
575 MakeSeeds(seeds, lsec, nos, nrow_up-1, nrow_up-1-gap);
576 MakeSeeds(seeds, lsec, nos, nrow_up-1-shift, nrow_up-1-shift-gap);
581 int nseed=seeds.GetEntriesFast();
583 for (int s=0; s<nseed; s++) {
584 AliTPCtrack& t=*((AliTPCtrack*)seeds.UncheckedAt(s));
585 double alpha=t.GetAlpha();
586 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
587 if (alpha < 0. ) alpha += 2.*TMath::Pi();
588 int ns=int(alpha/lsec->GetAlpha())%nos;
590 if (!FindProlongation(t,lsec,ns)) continue;
592 alpha=t.GetAlpha() + 0.5*ssec->GetAlpha() - ssec->GetAlphaShift();
593 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
594 if (alpha < 0. ) alpha += 2.*TMath::Pi();
595 ns=int(alpha/ssec->GetAlpha())%nis; //index of the inner sector needed
597 alpha=ns*ssec->GetAlpha() - t.GetAlpha();
598 if (!t.Rotate(alpha)) continue;
600 if (!FindProlongation(t,ssec,ns)) continue;
602 if (t < int(0.4*nrows)) continue;
613 //_____________________________________________________________________________
614 void AliTPC::CreateMaterials()
616 //-----------------------------------------------
617 // Create Materials for for TPC
618 //-----------------------------------------------
620 //-----------------------------------------------------------------
621 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
622 //-----------------------------------------------------------------
624 Int_t ISXFLD=gAlice->Field()->Integ();
625 Float_t SXMGMX=gAlice->Field()->Max();
627 Float_t amat[5]; // atomic numbers
628 Float_t zmat[5]; // z
629 Float_t wmat[5]; // proportions
633 // ********************* Gases *******************
635 //--------------------------------------------------------------
637 //--------------------------------------------------------------
642 Float_t a_ne = 20.18;
647 AliMaterial(20,"Ne",a_ne,z_ne,density,999.,999.);
651 Float_t a_ar = 39.948;
656 AliMaterial(21,"Ar",a_ar,z_ar,density,999.,999.);
664 //--------------------------------------------------------------
666 //--------------------------------------------------------------
683 amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
685 AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
700 amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
702 AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
717 amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
719 AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
721 //----------------------------------------------------------------
722 // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
723 //----------------------------------------------------------------
731 Float_t a,z,rho,absl,X0,buf[1];
734 for(nc = 0;nc<fNoComp;nc++)
737 // retrive material constants
739 gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
744 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
746 am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc]);
747 density += fMixtProp[nc]*rho; // density of the mixture
751 // mixture proportions by weight!
753 for(nc = 0;nc<fNoComp;nc++)
756 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
758 wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc])/am;
762 AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
763 AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
764 AliMixture(33,"Drift gas 3",amat,zmat,density,fNoComp,wmat);
766 AliMedium(2, "Drift gas 1", 31, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
767 AliMedium(3, "Drift gas 2", 32, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
768 AliMedium(4, "Drift gas 3", 33, 1, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
772 AliMaterial(24, "Air", 14.61, 7.3, .001205, 30420., 67500.);
774 AliMedium(24, "Air", 24, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
776 //----------------------------------------------------------------------
778 //----------------------------------------------------------------------
782 AliMaterial(30, "Al", 26.98, 13., 2.7, 8.9, 37.2);
784 AliMedium(0, "Al",30, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
788 AliMaterial(31, "Si", 28.086, 14.,2.33, 9.36, 999.);
790 AliMedium(7, "Al",31, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
809 AliMixture(32, "Mylar",amat,zmat,density,-3,wmat);
811 AliMedium(5, "Mylar",32, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
818 AliMaterial(33,"C normal",12.011,6.,2.265,18.8,999.);
820 AliMedium(6,"C normal",33,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
822 // G10 for inner and outr field cage
823 // G10 is 60% SiO2 + 40% epoxy, right now I use A and Z for SiO2
839 AliMixture(34,"G10 aux.",amat,zmat,density,-2,wmat);
842 gMC->Gfmate((*fIdmate)[34],namate,a,z,rho,X0,absl,buf,nbuf);
844 Float_t thickX0 = 0.0052; // field cage in X0 units
846 Float_t thick = 2.; // in cm
850 rhoFactor = X0*thickX0/thick;
851 density = rho*rhoFactor;
853 AliMaterial(35,"G10-fc",a,z,density,999.,999.);
855 AliMedium(8,"G10-fc",35,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
857 thickX0 = 0.0027; // inner vessel (eta <0.9)
859 rhoFactor = X0*thickX0/thick;
860 density = rho*rhoFactor;
862 AliMaterial(36,"G10-iv",a,z,density,999.,999.);
864 AliMedium(9,"G10-iv",36,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
868 gMC->Gfmate((*fIdmate)[33],namate,a,z,rho,X0,absl,buf,nbuf);
870 thickX0 = 0.0133; // outer vessel
872 rhoFactor = X0*thickX0/thick;
873 density = rho*rhoFactor;
876 AliMaterial(37,"C-ov",a,z,density,999.,999.);
878 AliMedium(10,"C-ov",37,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
880 thickX0=0.015; // inner vessel (cone, eta > 0.9)
882 rhoFactor = X0*thickX0/thick;
883 density = rho*rhoFactor;
885 AliMaterial(38,"C-ivc",a,z,density,999.,999.);
887 AliMedium(11,"C-ivc",38,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
891 AliMedium(12,"CO2",10,0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
897 //_____________________________________________________________________________
899 const AliTPCdigit *dig;
901 Bin() {dig=0; idx=-1;}
904 struct PreCluster : public AliTPCcluster {
905 const AliTPCdigit* summit; //pointer to the maximum digit of this precluster
906 int idx; //index in AliTPC::fClusters
907 int npeaks; //number of peaks in this precluster
908 int ndigits; //number of digits in this precluster
911 PreCluster::PreCluster() : AliTPCcluster() {npeaks=ndigits=0;}
913 //_____________________________________________________________________________
914 static void FindPreCluster(int i,int j,int maxj,Bin *bins,PreCluster &c)
916 //-----------------------------------------------------------------
917 // This function looks for "preclusters".
919 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
920 //-----------------------------------------------------------------
921 Bin& b=bins[i*maxj+j];
922 double q=(double)TMath::Abs(b.dig->fSignal);
924 if (b.idx >= 0 && b.idx != c.idx) {
929 if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
938 b.dig = 0; b.idx = c.idx;
940 if (bins[(i-1)*maxj+j].dig) FindPreCluster(i-1,j,maxj,bins,c);
941 if (bins[i*maxj+(j-1)].dig) FindPreCluster(i,j-1,maxj,bins,c);
942 if (bins[(i+1)*maxj+j].dig) FindPreCluster(i+1,j,maxj,bins,c);
943 if (bins[i*maxj+(j+1)].dig) FindPreCluster(i,j+1,maxj,bins,c);
947 //_____________________________________________________________________________
948 void AliTPC::Digits2Clusters()
950 //-----------------------------------------------------------------
951 // This is a simple cluster finder.
953 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
954 //-----------------------------------------------------------------
955 AliTPCParam *par = &(fDigParam->GetParam());
957 int inp=par->GetNPads(0, par->GetNRowLow()-1);
958 int onp=par->GetNPads(par->GetNSector()-1,par->GetNRowUp() -1);
959 const int MAXY=(inp>onp) ? inp+2 : onp+2;
960 const int MAXTBKT=int((z_end+6*par->GetZSigma())/par->GetZWidth())+1;
961 const int MAXZ=MAXTBKT+2;
962 const int THRESHOLD=20;
964 TTree *t=(TTree*)gDirectory->Get("TreeD0_Param1");
965 t->GetBranch("Digits")->SetAddress(&fDigits);
966 Int_t sectors_by_rows=(Int_t)t->GetEntries();
970 for (Int_t n=0; n<sectors_by_rows; n++) {
971 if (!t->GetEvent(n)) continue;
972 Bin *bins=new Bin[MAXY*MAXZ];
973 AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
974 int sec=dig->fSector, row=dig->fPadRow;
975 int ndigits=fDigits->GetEntriesFast();
979 int nis=par->GetNInnerSector(), nos=par->GetNOuterSector();
981 npads = par->GetNPadsLow(row);
982 sign = (sec < nis/2) ? 1 : -1;
984 npads = par->GetNPadsUp(row);
985 sign = ((sec-nis) < nos/2) ? 1 : -1;
990 for (ndig=0; ndig<ndigits; ndig++) {
991 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
992 int i=dig->fPad+1, j=dig->fTime+1;
994 cerr<<"AliTPC::Digits2Clusters error: pad number is out of range ! ";
995 cerr<<i<<' '<<npads<<endl;
999 cerr<<"AliTPC::Digits2Clusters error: time bucket is out of range ! ";
1000 cerr<<j<<' '<<MAXTBKT<<endl;
1003 if (dig->fSignal >= THRESHOLD) bins[i*MAXZ+j].dig=dig;
1004 if (i==1 || i==npads || j==1 || j==MAXTBKT) dig->fSignal*=-1;
1010 for (i=1; i<MAXY-1; i++) {
1011 for (j=1; j<MAXZ-1; j++) {
1012 if (bins[i*MAXZ+j].dig == 0) continue;
1013 PreCluster c; c.summit=bins[i*MAXZ+j].dig; c.idx=ncls;
1014 FindPreCluster(i, j, MAXZ, bins, c);
1018 double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
1019 c.fSigmaY2 = s2 + 1./12.;
1020 c.fSigmaY2 *= par->GetPadPitchWidth()*par->GetPadPitchWidth();
1021 if (s2 != 0.) c.fSigmaY2 *= 0.17;
1023 s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
1024 c.fSigmaZ2 = s2 + 1./12.;
1025 c.fSigmaZ2 *= par->GetZWidth()*par->GetZWidth();
1026 if (s2 != 0.) c.fSigmaZ2 *= 0.41;
1028 c.fY = (c.fY - 0.5 - 0.5*npads)*par->GetPadPitchWidth();
1029 c.fZ = par->GetZWidth()*c.fZ;
1030 c.fZ -= 3.*par->GetZSigma(); // PASA delay
1031 c.fZ = sign*(z_end - c.fZ);
1035 c.fTracks[0]=c.summit->fTracks[0];
1036 c.fTracks[1]=c.summit->fTracks[1];
1037 c.fTracks[2]=c.summit->fTracks[2];
1039 if (c.summit->fSignal<0) {
1044 AddCluster(c); ncls++; ncl++;
1048 for (ndig=0; ndig<ndigits; ndig++) {
1049 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
1050 int i=dig->fPad+1, j=dig->fTime+1;
1052 cerr<<"AliTPC::Digits2Clusters error: pad number is out of range ! ";
1053 cerr<<i<<' '<<npads<<endl;
1057 cerr<<"AliTPC::Digits2Clusters error: time bucket is out of range ! ";
1058 cerr<<j<<' '<<MAXTBKT<<endl;
1061 if (TMath::Abs(dig->fSignal)>=par->GetZeroSup()) bins[i*MAXZ+j].dig=dig;
1064 for (i=1; i<MAXY-1; i++) {
1065 for (j=1; j<MAXZ-1; j++) {
1066 if (bins[i*MAXZ+j].dig == 0) continue;
1067 PreCluster c; c.summit=bins[i*MAXZ+j].dig; c.idx=ncls;
1068 FindPreCluster(i, j, MAXZ, bins, c);
1069 if (c.ndigits < 2) continue; //noise cluster
1070 if (c.npeaks>1) continue; //overlapped cluster
1074 double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
1075 c.fSigmaY2 = s2 + 1./12.;
1076 c.fSigmaY2 *= par->GetPadPitchWidth()*par->GetPadPitchWidth();
1077 if (s2 != 0.) c.fSigmaY2 *= 0.04;
1079 s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
1080 c.fSigmaZ2 = s2 + 1./12.;
1081 c.fSigmaZ2 *= par->GetZWidth()*par->GetZWidth();
1082 if (s2 != 0.) c.fSigmaZ2 *= 0.10;
1084 c.fY = (c.fY - 0.5 - 0.5*npads)*par->GetPadPitchWidth();
1085 c.fZ = par->GetZWidth()*c.fZ;
1086 c.fZ -= 3.*par->GetZSigma(); // PASA delay
1087 c.fZ = sign*(z_end - c.fZ);
1091 c.fTracks[0]=c.summit->fTracks[0];
1092 c.fTracks[1]=c.summit->fTracks[1];
1093 c.fTracks[2]=c.summit->fTracks[2];
1095 if (c.summit->fSignal<0) {
1100 if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;}
1102 new ((*fClusters)[c.idx]) AliTPCcluster(c);
1107 cerr<<"sector, row, digits, clusters: "
1108 <<sec<<' '<<row<<' '<<ndigits<<' '<<ncl<<" \r";
1116 //_____________________________________________________________________________
1117 void AliTPC::ElDiff(Float_t *xyz)
1119 //--------------------------------------------------
1120 // calculates the diffusion of a single electron
1121 //--------------------------------------------------
1123 //-----------------------------------------------------------------
1124 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1125 //-----------------------------------------------------------------
1126 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1131 driftl=z_end-TMath::Abs(xyz[2]);
1133 if(driftl<0.01) driftl=0.01;
1135 // check the attachment
1137 driftl=TMath::Sqrt(driftl);
1139 // Float_t sig_t = driftl*diff_t;
1140 //Float_t sig_l = driftl*diff_l;
1141 Float_t sig_t = driftl*fTPCParam->GetDiffT();
1142 Float_t sig_l = driftl*fTPCParam->GetDiffL();
1143 xyz[0]=gRandom->Gaus(xyz[0],sig_t);
1144 xyz[1]=gRandom->Gaus(xyz[1],sig_t);
1145 xyz[2]=gRandom->Gaus(xyz[2],sig_l);
1147 if (TMath::Abs(xyz[2])>z_end){
1148 xyz[2]=TMath::Sign(z_end,z0);
1151 Float_t eps = 0.0001;
1152 xyz[2]=TMath::Sign(eps,z0);
1156 //_____________________________________________________________________________
1157 void AliTPC::Hits2Clusters()
1159 //--------------------------------------------------------
1160 // TPC simple cluster generator from hits
1161 // obtained from the TPC Fast Simulator
1162 // The point errors are taken from the parametrization
1163 //--------------------------------------------------------
1165 //-----------------------------------------------------------------
1166 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1167 //-----------------------------------------------------------------
1169 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1170 Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
1172 TParticle *particle; // pointer to a given particle
1173 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1174 TClonesArray *Particles; //pointer to the particle list
1178 Float_t pl,pt,tanth,rpad,ratio;
1181 //---------------------------------------------------------------
1182 // Get the access to the tracks
1183 //---------------------------------------------------------------
1185 TTree *TH = gAlice->TreeH();
1186 Stat_t ntracks = TH->GetEntries();
1187 Particles=gAlice->Particles();
1189 //------------------------------------------------------------
1190 // Loop over all sectors (72 sectors)
1191 // Sectors 0-35 are lower sectors, 0-17 z>0, 18-35 z<0
1192 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
1194 // First cluster for sector 0 starts at "0"
1195 //------------------------------------------------------------
1199 int Nsectors=fDigParam->GetParam().GetNSector();
1200 for(Int_t isec=0; isec<Nsectors; isec++){
1202 fTPCParam->AdjustAngles(isec,cph,sph);
1204 //------------------------------------------------------------
1206 //------------------------------------------------------------
1208 for(Int_t track=0;track<ntracks;track++){
1210 TH->GetEvent(track);
1212 // Get number of the TPC hits and a pointer
1215 nhits=fHits->GetEntriesFast();
1219 for(Int_t hit=0;hit<nhits;hit++){
1220 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1221 sector=tpcHit->fSector; // sector number
1222 if(sector != isec) continue; //terminate iteration
1223 ipart=tpcHit->fTrack;
1224 particle=(TParticle*)Particles->UncheckedAt(ipart);
1227 if(pt < 1.e-9) pt=1.e-9;
1229 tanth = TMath::Abs(tanth);
1230 rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
1231 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1233 // space-point resolutions
1235 sigma_rphi=SigmaY2(rpad,tanth,pt);
1236 sigma_z =SigmaZ2(rpad,tanth );
1240 cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
1241 cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
1243 // temporary protection
1245 if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
1246 if(sigma_z < 0.) sigma_z=0.4e-3;
1247 if(cl_rphi < 0.) cl_rphi=2.5e-3;
1248 if(cl_z < 0.) cl_z=2.5e-5;
1251 // smearing --> rotate sectors firstly,
1252 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1254 Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
1255 Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
1256 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi)); // y
1257 Double_t alpha=(sector < fTPCParam->GetNInnerSector()) ?
1258 fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
1259 if (TMath::Abs(xyz[0]/xprim) > TMath::Tan(0.5*alpha)) xyz[0]=yprim;
1260 xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z
1261 if (TMath::Abs(xyz[1]) > 250) xyz[1]=tpcHit->fZ;
1262 xyz[2]=tpcHit->fQ+1;// q; let it be not equal to zero (Y.Belikov)
1263 xyz[3]=sigma_rphi; // fSigmaY2
1264 xyz[4]=sigma_z; // fSigmaZ2
1266 // and finally add the cluster
1267 Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, tpcHit->fPadRow};
1268 AddCluster(xyz,tracks);
1270 } // end of loop over hits
1271 } // end of loop over tracks
1274 } // end of loop over sectors
1280 void AliTPC::Hits2Digits()
1283 //----------------------------------------------------
1284 // Loop over all sectors (72 sectors)
1285 // Sectors 0-35 are lower sectors, 0-17 z>0, 18-35 z<0
1286 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
1288 int Nsectors=fDigParam->GetParam().GetNSector();
1289 for(Int_t isec=0;isec<Nsectors;isec++) Hits2DigitsSector(isec);
1293 //_____________________________________________________________________________
1294 void AliTPC::Hits2DigitsSector(Int_t isec)
1296 //-------------------------------------------------------------------
1297 // TPC conversion from hits to digits.
1298 //-------------------------------------------------------------------
1300 //-----------------------------------------------------------------
1301 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1302 //-----------------------------------------------------------------
1304 //-------------------------------------------------------
1305 // Get the access to the track hits
1306 //-------------------------------------------------------
1308 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1309 TTree *TH = gAlice->TreeH(); // pointer to the hits tree
1310 Stat_t ntracks = TH->GetEntries();
1314 //-------------------------------------------
1315 // Only if there are any tracks...
1316 //-------------------------------------------
1319 // TObjArrays for three neighouring pad-rows
1321 TObjArray **rowTriplet = new TObjArray* [3];
1323 // TObjArray-s for each pad-row
1327 for (Int_t trip=0;trip<3;trip++){
1328 rowTriplet[trip]=new TObjArray;
1333 printf("*** Processing sector number %d ***\n",isec);
1335 Int_t nrows =fTPCParam->GetNRow(isec);
1337 row= new TObjArray* [nrows];
1339 MakeSector(isec,nrows,TH,ntracks,row);
1341 //--------------------------------------------------------
1342 // Digitize this sector, row by row
1343 // row[i] is the pointer to the TObjArray of TVectors,
1344 // each one containing electrons accepted on this
1345 // row, assigned into tracks
1346 //--------------------------------------------------------
1350 for (i=0;i<nrows;i++){
1352 // Triplets for i = 0 and i=1 are identical!
1353 // The same for i = nrows-1 and nrows!
1355 if(i != 1 && i != nrows-1){
1356 MakeTriplet(i,rowTriplet,row);
1359 DigitizeRow(i,isec,rowTriplet);
1363 Int_t ndig=fDigParam->GetArray()->GetEntriesFast();
1365 printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig);
1367 ResetDigits(); // reset digits for this row after storing them
1369 } // end of the sector digitization
1371 // delete the last triplet
1373 for (i=0;i<3;i++) rowTriplet[i]->Delete();
1375 delete [] row; // delete the array of pointers to TObjArray-s
1378 } // end of Hits2Digits
1379 //_____________________________________________________________________________
1380 void AliTPC::MakeTriplet(Int_t row,
1381 TObjArray **rowTriplet, TObjArray **prow)
1383 //------------------------------------------------------------------
1384 // Makes the "triplet" of the neighbouring pad-row for the
1385 // digitization including the cross-talk between the pad-rows
1386 //------------------------------------------------------------------
1388 //-----------------------------------------------------------------
1389 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1390 //-----------------------------------------------------------------
1392 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1393 Float_t gasgain = fTPCParam->GetGasGain();
1396 Int_t nElements,nElectrons;
1401 //-------------------------------------------------------------------
1402 // pv is an "old" track, i.e. label + triplets of (x,y,z)
1403 // for each electron
1405 //-------------------------------------------------------------------
1411 if(row == 0 || row == 1){
1413 // create entire triplet for the first AND the second row
1415 nTracks[0] = prow[0]->GetEntries();
1416 nTracks[1] = prow[1]->GetEntries();
1417 nTracks[2] = prow[2]->GetEntries();
1419 for(i1=0;i1<3;i1++){
1420 nt = nTracks[i1]; // number of tracks for this row
1422 for(i2=0;i2<nt;i2++){
1423 pv = (TVector*)prow[i1]->At(i2);
1425 nElements = pv->GetNrows();
1426 nElectrons = (nElements-1)/3;
1428 tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1430 v2(0)=v1(0); //track label
1432 for(nel=0;nel<nElectrons;nel++){
1435 // Avalanche, including fluctuations
1436 Int_t aval = (Int_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1437 v2(idx2+1)= v1(idx1+1);
1438 v2(idx2+2)= v1(idx1+2);
1439 v2(idx2+3)= v1(idx1+3);
1440 v2(idx2+4)= (Float_t)aval; // in number of electrons!
1441 } // end of loop over electrons
1443 // Add this track to a row
1446 rowTriplet[i1]->Add(tv);
1449 } // end of loop over tracks for this row
1451 prow[i1]->Delete(); // remove "old" tracks
1452 delete prow[i1]; // delete TObjArray for this row
1453 prow[i1]=0; // set pointer to NULL
1455 } // end of loop over row triplets
1461 rowTriplet[0]->Delete(); // remove old lower row
1463 nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row
1464 nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row
1465 nTracks[2]=prow[row+1]->GetEntries(); // next row
1468 //-------------------------------------------
1469 // shift new tracks downwards
1470 //-------------------------------------------
1472 for(i1=0;i1<nTracks[0];i1++){
1473 pv=(TVector*)rowTriplet[1]->At(i1);
1474 rowTriplet[0]->Add(pv);
1477 rowTriplet[1]->Clear(); // leave tracks on the heap!!!
1479 for(i1=0;i1<nTracks[1];i1++){
1480 pv=(TVector*)rowTriplet[2]->At(i1);
1481 rowTriplet[1]->Add(pv);
1484 rowTriplet[2]->Clear(); // leave tracks on the heap!!!
1486 //---------------------------------------------
1487 // Create new upper row
1488 //---------------------------------------------
1492 for(i1=0;i1<nTracks[2];i1++){
1493 pv = (TVector*)prow[row+1]->At(i1);
1495 nElements = pv->GetNrows();
1496 nElectrons = (nElements-1)/3;
1498 tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1500 v2(0)=v1(0); //track label
1502 for(nel=0;nel<nElectrons;nel++){
1506 // Avalanche, including fluctuations
1507 Int_t aval = (Int_t)
1508 (-gasgain*TMath::Log(gRandom->Rndm()));
1510 v2(idx2+1)= v1(idx1+1);
1511 v2(idx2+2)= v1(idx1+2);
1512 v2(idx2+3)= v1(idx1+3);
1513 v2(idx2+4)= (Float_t)aval; // in number of electrons!
1514 } // end of loop over electrons
1516 rowTriplet[2]->Add(tv);
1518 } // end of loop over tracks
1520 prow[row+1]->Delete(); // delete tracks for this row
1521 delete prow[row+1]; // delete TObjArray for this row
1522 prow[row+1]=0; // set a pointer to NULL
1526 } // end of MakeTriplet
1527 //_____________________________________________________________________________
1528 void AliTPC::ExB(Float_t *xyz)
1530 //------------------------------------------------
1531 // ExB at the wires and wire number calulation
1532 //------------------------------------------------
1534 //-----------------------------------------------------------------
1535 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1536 //-----------------------------------------------------------------
1537 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1540 fTPCParam->GetWire(x1); //calculate nearest wire position
1541 Float_t dx=xyz[0]-x1;
1542 xyz[1]+=dx*fTPCParam->GetOmegaTau();
1545 //_____________________________________________________________________________
1546 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet)
1548 //-----------------------------------------------------------
1549 // Single row digitization, coupling from the neighbouring
1550 // rows taken into account
1551 //-----------------------------------------------------------
1553 //-----------------------------------------------------------------
1554 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1555 //-----------------------------------------------------------------
1558 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1559 Float_t chipgain= fTPCParam->GetChipGain();
1560 Float_t zerosup = fTPCParam->GetZeroSup();
1561 Int_t nrows =fTPCParam->GetNRow(isec);
1563 int((z_end+6*fTPCParam->GetZSigma())/fTPCParam->GetZWidth())+1;
1567 Int_t IndexRange[4];
1572 // iFlag = 0 -> inner row, iFlag = 1 -> the middle one, iFlag = 2 -> the outer one
1575 nTracks[0]=rowTriplet[0]->GetEntries(); // lower row
1576 nTracks[1]=rowTriplet[1]->GetEntries(); // middle row
1577 nTracks[2]=rowTriplet[2]->GetEntries(); // upper row
1582 n_of_pads[0]=fTPCParam->GetNPads(isec,0);
1583 n_of_pads[1]=fTPCParam->GetNPads(isec,1);
1585 else if(irow == nrows-1){
1587 n_of_pads[1]=fTPCParam->GetNPads(isec,irow-1);
1588 n_of_pads[2]=fTPCParam->GetNPads(isec,irow);
1592 for(i1=0;i1<3;i1++){
1593 n_of_pads[i1]=fTPCParam->GetNPads(isec,irow-1+i1);
1598 // Integrated signal for this row
1599 // and a single track signal
1602 TMatrix *m1 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTBKT-1); // integrated
1603 TMatrix *m2 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTBKT-1); // single
1607 TMatrix &Total = *m1;
1609 // Array of pointers to the label-signal list
1611 Int_t NofDigits = n_of_pads[iFlag]*MAXTBKT; // number of digits for this row
1613 Float_t **pList = new Float_t* [NofDigits];
1617 for(lp=0;lp<NofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1620 // Straight signal and cross-talk, cross-talk is integrated over all
1621 // tracks and added to the signal at the very end
1625 for (i1=0;i1<nTracks[iFlag];i1++){
1627 m2->Zero(); // clear single track signal matrix
1629 Float_t TrackLabel =
1630 GetSignal(rowTriplet[iFlag],i1,n_of_pads[iFlag],m2,m1,IndexRange);
1632 GetList(TrackLabel,n_of_pads[iFlag],m2,IndexRange,pList);
1637 // Cross talk from the neighbouring pad-rows
1640 TMatrix *m3 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTBKT-1); // cross-talk
1642 TMatrix &Cross = *m3;
1646 // cross-talk from the outer row only (first pad row)
1648 GetCrossTalk(0,rowTriplet[1],nTracks[1],n_of_pads,m3);
1651 else if(iFlag == 2){
1653 // cross-talk from the inner row only (last pad row)
1655 GetCrossTalk(2,rowTriplet[1],nTracks[1],n_of_pads,m3);
1660 // cross-talk from both inner and outer rows
1662 GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner
1663 GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer
1666 Total += Cross; // add the cross-talk
1669 // Convert analog signal to ADC counts
1676 for(Int_t ip=0;ip<n_of_pads[iFlag];ip++){
1677 for(Int_t it=0;it<MAXTBKT;it++){
1679 Float_t q = Total(ip,it);
1681 Int_t gi =it*n_of_pads[iFlag]+ip; // global index
1683 q = gRandom->Gaus(q,fTPCParam->GetNoise()); // apply noise
1684 q *= (q_el*1.e15); // convert to fC
1685 q *= chipgain; // convert to mV
1686 q *= (adc_sat/dyn_range); // convert to ADC counts
1688 if(q <zerosup) continue; // do not fill zeros
1689 if(q > adc_sat) q = adc_sat; // saturation
1692 // "real" signal or electronic noise (list = -1)?
1695 for(Int_t j1=0;j1<3;j1++){
1696 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1703 digits[4]= (Int_t)q;
1707 AddDigit(tracks,digits);
1709 } // end of loop over time buckets
1710 } // end of lop over pads
1713 // This row has been digitized, delete nonused stuff
1716 for(lp=0;lp<NofDigits;lp++){
1717 if(pList[lp]) delete [] pList[lp];
1726 } // end of DigitizeRow
1727 //_____________________________________________________________________________
1728 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, Int_t np, TMatrix *m1, TMatrix *m2,
1732 //---------------------------------------------------------------
1733 // Calculates 2-D signal (pad,time) for a single track,
1734 // returns a pointer to the signal matrix and the track label
1735 // No digitization is performed at this level!!!
1736 //---------------------------------------------------------------
1738 //-----------------------------------------------------------------
1739 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1740 //-----------------------------------------------------------------
1743 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1744 AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
1745 AliTPCRF1D * fRF = &(fDigParam->GetRF());
1747 int((z_end+6*fTPCParam->GetZSigma())/fTPCParam->GetZWidth())+1;
1749 //to make the code faster we put parameters to the stack
1751 Float_t zwidth = fTPCParam->GetZWidth();
1752 Float_t zwidthm1 =1./zwidth;
1754 tv = (TVector*)p1->At(ntr); // pointer to a track
1757 Float_t label = v(0);
1759 Int_t CentralPad = (np-1)/2;
1761 Int_t nElectrons = (tv->GetNrows()-1)/4;
1762 Float_t range=((np-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); // pad range
1764 range -= 0.5; // dead zone, 5mm from the edge, according to H.G. Fischer
1766 Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1769 Float_t PadSignal[7]; // signal from a single electron
1771 TMatrix &signal = *m1;
1772 TMatrix &total = *m2;
1775 IndexRange[0]=9999; // min pad
1776 IndexRange[1]=-1; // max pad
1777 IndexRange[2]=9999; //min time
1778 IndexRange[3]=-1; // max time
1781 // Loop over all electrons
1784 for(Int_t nel=0; nel<nElectrons; nel++){
1786 Float_t xwire = v(idx+1);
1787 Float_t y = v(idx+2);
1788 Float_t z = v(idx+3);
1791 Float_t absy=TMath::Abs(y);
1793 if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
1794 PadNumber=CentralPad;
1796 else if (absy < range){
1797 PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth()+1.);
1798 PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
1800 else continue; // electron out of pad-range , lost at the sector edge
1802 Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
1805 Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
1806 for (Int_t i=0;i<7;i++){
1807 PadSignal[i]=fPRF2D->GetPRF(-dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
1808 PadSignal[i] *= fTPCParam->GetPadCoupling();
1811 Int_t LeftPad = TMath::Max(0,PadNumber-3);
1812 Int_t RightPad = TMath::Min(np-1,PadNumber+3);
1814 Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
1815 Int_t pmax=RightPad-PadNumber+3; // upper index
1817 Float_t z_drift = z*zwidthm1;
1818 Float_t z_offset = z_drift-(Int_t)z_drift;
1820 Int_t FirstBucket = (Int_t)z_drift; // numbering starts from "0"
1823 // loop over time bins (4 bins is enough - 3 sigma truncated Gaussian)
1825 for (Int_t i2=0;i2<4;i2++){
1826 Int_t TrueTime = FirstBucket+i2; // current time bucket
1827 Float_t dz = (Float_t(i2)+1.-z_offset)*zwidth;
1828 Float_t ampl = fRF->GetRF(dz);
1829 if( (TrueTime>MAXTBKT-1) ) break; // beyond the time range
1831 IndexRange[2]=TMath::Min(IndexRange[2],TrueTime); // min time
1832 IndexRange[3]=TMath::Max(IndexRange[3],TrueTime); // max time
1834 // loop over pads, from pmin to pmax
1835 for(Int_t i3=pmin;i3<=pmax;i3++){
1836 Int_t TruePad = LeftPad+i3-pmin;
1837 IndexRange[0]=TMath::Min(IndexRange[0],TruePad); // min pad
1838 IndexRange[1]=TMath::Max(IndexRange[1],TruePad); // max pad
1839 signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!!
1840 total(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!!
1841 } // end of pads loop
1842 } // end of time loop
1843 } // end of loop over electrons
1845 return label; // returns track label when finished
1848 //_____________________________________________________________________________
1849 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange,
1852 //----------------------------------------------------------------------
1853 // Updates the list of tracks contributing to digits for a given row
1854 //----------------------------------------------------------------------
1856 //-----------------------------------------------------------------
1857 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1858 //-----------------------------------------------------------------
1860 TMatrix &signal = *m;
1862 // lop over nonzero digits
1864 for(Int_t it=IndexRange[2];it<IndexRange[3]+1;it++){
1865 for(Int_t ip=IndexRange[0];ip<IndexRange[1]+1;ip++){
1868 Int_t GlobalIndex = it*np+ip; // GlobalIndex starts from 0!
1870 if(!pList[GlobalIndex]){
1873 // Create new list (6 elements - 3 signals and 3 labels),
1874 // but only if the signal is > 0.
1877 if(signal(ip,it)>0.){
1879 pList[GlobalIndex] = new Float_t [6];
1883 *pList[GlobalIndex] = -1.;
1884 *(pList[GlobalIndex]+1) = -1.;
1885 *(pList[GlobalIndex]+2) = -1.;
1886 *(pList[GlobalIndex]+3) = -1.;
1887 *(pList[GlobalIndex]+4) = -1.;
1888 *(pList[GlobalIndex]+5) = -1.;
1891 *pList[GlobalIndex] = label;
1892 *(pList[GlobalIndex]+3) = signal(ip,it);}
1896 // check the signal magnitude
1898 Float_t highest = *(pList[GlobalIndex]+3);
1899 Float_t middle = *(pList[GlobalIndex]+4);
1900 Float_t lowest = *(pList[GlobalIndex]+5);
1903 // compare the new signal with already existing list
1906 if(signal(ip,it)<lowest) continue; // neglect this track
1910 if (signal(ip,it)>highest){
1911 *(pList[GlobalIndex]+5) = middle;
1912 *(pList[GlobalIndex]+4) = highest;
1913 *(pList[GlobalIndex]+3) = signal(ip,it);
1915 *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1916 *(pList[GlobalIndex]+1) = *pList[GlobalIndex];
1917 *pList[GlobalIndex] = label;
1919 else if (signal(ip,it)>middle){
1920 *(pList[GlobalIndex]+5) = middle;
1921 *(pList[GlobalIndex]+4) = signal(ip,it);
1923 *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1924 *(pList[GlobalIndex]+1) = label;
1927 *(pList[GlobalIndex]+5) = signal(ip,it);
1928 *(pList[GlobalIndex]+2) = label;
1932 } // end of loop over pads
1933 } // end of loop over time bins
1939 //___________________________________________________________________
1940 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1941 Stat_t ntracks,TObjArray **row)
1944 //-----------------------------------------------------------------
1945 // Prepares the sector digitization, creates the vectors of
1946 // tracks for each row of this sector. The track vector
1947 // contains the track label and the position of electrons.
1948 //-----------------------------------------------------------------
1950 //-----------------------------------------------------------------
1951 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1952 //-----------------------------------------------------------------
1954 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1958 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1960 //----------------------------------------------
1961 // Create TObjArray-s, one for each row,
1962 // each TObjArray will store the TVectors
1963 // of electrons, one TVector per each track.
1964 //----------------------------------------------
1966 for(i=0; i<nrows; i++){
1967 row[i] = new TObjArray;
1969 Int_t *n_of_electrons = new Int_t [nrows]; // electron counter for each row
1970 TVector **tr = new TVector* [nrows]; //pointers to the track vectors
1972 //--------------------------------------------------------------------
1973 // Loop over tracks, the "track" contains the full history
1974 //--------------------------------------------------------------------
1976 Int_t previousTrack,currentTrack;
1977 previousTrack = -1; // nothing to store so far!
1979 for(Int_t track=0;track<ntracks;track++){
1983 TH->GetEvent(track); // get next track
1984 Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1986 if(nhits == 0) continue; // no hits in the TPC for this track
1988 //--------------------------------------------------------------
1990 //--------------------------------------------------------------
1992 for(Int_t hit=0;hit<nhits;hit++){
1994 tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1996 Int_t sector=tpcHit->fSector; // sector number
1997 if(sector != isec) continue; //terminate iteration
1999 currentTrack = tpcHit->fTrack; // track number
2000 if(currentTrack != previousTrack){
2002 // store already filled fTrack
2004 for(i=0;i<nrows;i++){
2005 if(previousTrack != -1){
2006 if(n_of_electrons[i]>0){
2007 TVector &v = *tr[i];
2008 v(0) = previousTrack;
2009 tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary
2013 delete tr[i]; // delete empty TVector
2018 n_of_electrons[i]=0;
2019 tr[i] = new TVector(361); // TVectors for the next fTrack
2021 } // end of loop over rows
2023 previousTrack=currentTrack; // update track label
2026 Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2028 //---------------------------------------------------
2029 // Calculate the electron attachment probability
2030 //---------------------------------------------------
2032 Float_t time = 1.e6*(z_end-TMath::Abs(tpcHit->fZ))/fTPCParam->GetDriftV();
2034 Float_t AttProb = fTPCParam->GetAttCoef()*
2035 fTPCParam->GetOxyCont()*time; // fraction!
2037 //-----------------------------------------------
2038 // Loop over electrons
2039 //-----------------------------------------------
2041 for(Int_t nel=0;nel<QI;nel++){
2042 // skip if electron lost due to the attachment
2043 if((gRandom->Rndm(0)) < AttProb) continue; // electron lost!
2047 ElDiff(xyz); // Appply the diffusion
2049 fTPCParam->XYZtoCRXYZ(xyz,isec,row_number,3);
2051 //transform position to local coordinates
2052 //option 3 means that we calculate x position relative to
2055 if ((row_number<0)||row_number>=fTPCParam->GetNRow(isec)) continue;
2056 ExB(xyz); // ExB effect at the sense wires
2057 n_of_electrons[row_number]++;
2058 //----------------------------------
2059 // Expand vector if necessary
2060 //----------------------------------
2061 if(n_of_electrons[row_number]>120){
2062 Int_t range = tr[row_number]->GetNrows();
2063 if(n_of_electrons[row_number] > (range-1)/3){
2064 tr[row_number]->ResizeTo(range+150); // Add 50 electrons
2068 TVector &v = *tr[row_number];
2069 Int_t idx = 3*n_of_electrons[row_number]-2;
2071 v(idx)= xyz[0]; // X
2072 v(idx+1)=xyz[1]; // Y (along the pad-row)
2073 v(idx+2)=xyz[2]; // Z
2075 } // end of loop over electrons
2077 } // end of loop over hits
2078 } // end of loop over tracks
2081 // store remaining track (the last one) if not empty
2084 for(i=0;i<nrows;i++){
2085 if(n_of_electrons[i]>0){
2086 TVector &v = *tr[i];
2087 v(0) = previousTrack;
2088 tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary
2098 delete [] n_of_electrons;
2100 } // end of MakeSector
2101 //_____________________________________________________________________________
2102 void AliTPC::GetCrossTalk (Int_t iFlag,TObjArray *p,Int_t ntracks,Int_t *npads,
2106 //-------------------------------------------------------------
2107 // Calculates the cross-talk from one row (inner or outer)
2108 //-------------------------------------------------------------
2110 //-----------------------------------------------------------------
2111 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2112 //-----------------------------------------------------------------
2115 // iFlag=2 & 3 --> cross-talk from the inner row
2116 // iFlag=0 & 4 --> cross-talk from the outer row
2118 AliTPCParam * fTPCParam = &(fDigParam->GetParam());
2119 AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
2120 AliTPCRF1D * fRF = &(fDigParam->GetRF());
2122 int((z_end+6*fTPCParam->GetZSigma())/fTPCParam->GetZWidth())+1;
2124 //to make code faster
2126 Float_t zwidth = fTPCParam->GetZWidth();
2127 Float_t zwidthm1 =1/fTPCParam->GetZWidth();
2132 Int_t nPadsSignal; // for this pads the signal is calculated
2133 Float_t range; // sense wire range
2136 Float_t IneffFactor=0.5; // gain inefficiency close to the edges
2140 nPadsSignal = npads[1];
2141 range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2142 nPadsDiff = (npads[1]-npads[0])/2;
2144 else if (iFlag == 2){
2146 nPadsSignal = npads[2];
2147 range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2150 else if (iFlag == 3){
2152 nPadsSignal = npads[1];
2153 range = ((npads[0]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2158 nPadsSignal = npads[2];
2159 range = ((npads[2]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth();
2160 nPadsDiff = (npads[2]-npads[1])/2;
2163 range-=0.5; // dead zone close to the edges
2166 TMatrix &signal = *m;
2168 Int_t CentralPad = (nPadsSignal-1)/2;
2169 Float_t PadSignal[7]; // signal from a single electron
2171 for(Int_t nt=0;nt<ntracks;nt++){
2172 tv=(TVector*)p->At(nt); // pointer to a track
2174 Int_t nElectrons = (tv->GetNrows()-1)/4;
2175 // Loop over electrons
2176 for(Int_t nel=0; nel<nElectrons; nel++){
2180 if (iFlag==0) xwire+=fTPCParam->GetPadPitchLength();
2181 if (iFlag==2) xwire-=fTPCParam->GetPadPitchLength();
2182 if (iFlag==3) xwire-=fTPCParam->GetPadPitchLength();
2183 if (iFlag==4) xwire+=fTPCParam->GetPadPitchLength();
2185 // electron acceptance for the cross-talk (only the closest wire)
2187 Float_t dxMax = fTPCParam->GetPadPitchLength()*0.5+fTPCParam->GetWWPitch();
2188 if(TMath::Abs(xwire)>dxMax) continue;
2190 Float_t y = v(idx+2);
2191 Float_t z = v(idx+3);
2192 Float_t absy=TMath::Abs(y);
2194 if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
2195 PadNumber=CentralPad;
2197 else if (absy < range){
2198 PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.);
2199 PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
2201 else continue; // electron out of sense wire range, lost at the sector edge
2203 Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
2205 Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
2207 for (Int_t i=0;i<7;i++){
2208 PadSignal[i]=fPRF2D->GetPRF(-dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
2210 PadSignal[i] *= fTPCParam->GetPadCoupling();
2214 Int_t LeftPad = TMath::Max(0,PadNumber-3);
2215 Int_t RightPad = TMath::Min(nPadsSignal-1,PadNumber+3);
2217 Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
2218 Int_t pmax=RightPad-PadNumber+3; // upper index
2221 Float_t z_drift = z*zwidthm1;
2222 Float_t z_offset = z_drift-(Int_t)z_drift;
2224 Int_t FirstBucket = (Int_t)z_drift; // numbering starts from "0"
2226 for (Int_t i2=0;i2<4;i2++){
2227 Int_t TrueTime = FirstBucket+i2; // current time bucket
2228 Float_t dz = (Float_t(i2)+1.- z_offset)*zwidth;
2229 Float_t ampl = fRF->GetRF(dz);
2230 if((TrueTime>MAXTBKT-1)) break; // beyond the time range
2233 // loop over pads, from pmin to pmax
2235 for(Int_t i3=pmin;i3<pmax+1;i3++){
2236 Int_t TruePad = LeftPad+i3-pmin;
2238 if(TruePad<nPadsDiff || TruePad > nPadsSignal-nPadsDiff-1) continue;
2240 TruePad -= nPadsDiff;
2241 signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!
2243 } // end of loop over pads
2244 } // end of loop over time bins
2246 } // end of loop over electrons
2248 } // end of loop over tracks
2250 } // end of GetCrossTalk
2254 //_____________________________________________________________________________
2258 // Initialise TPC detector after definition of geometry
2263 for(i=0;i<35;i++) printf("*");
2264 printf(" TPC_INIT ");
2265 for(i=0;i<35;i++) printf("*");
2268 for(i=0;i<80;i++) printf("*");
2272 //_____________________________________________________________________________
2273 void AliTPC::MakeBranch(Option_t* option)
2276 // Create Tree branches for the TPC.
2278 Int_t buffersize = 4000;
2279 char branchname[10];
2280 sprintf(branchname,"%s",GetName());
2282 AliDetector::MakeBranch(option);
2284 char *D = strstr(option,"D");
2286 if (fDigits && gAlice->TreeD() && D) {
2287 gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
2288 printf("Making Branch %s for digits\n",branchname);
2291 char *R = strstr(option,"R");
2293 if (fClusters && gAlice->TreeR() && R) {
2294 gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
2295 printf("Making Branch %s for Clusters\n",branchname);
2299 //_____________________________________________________________________________
2300 void AliTPC::ResetDigits()
2303 // Reset number of digits and the digits array for this detector
2307 // if (fDigits) fDigits->Clear();
2308 if (fDigParam->GetArray()!=0) fDigParam->GetArray()->Clear();
2310 if (fClusters) fClusters->Clear();
2313 //_____________________________________________________________________________
2314 void AliTPC::SetSecAL(Int_t sec)
2316 //---------------------------------------------------
2317 // Activate/deactivate selection for lower sectors
2318 //---------------------------------------------------
2320 //-----------------------------------------------------------------
2321 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2322 //-----------------------------------------------------------------
2327 //_____________________________________________________________________________
2328 void AliTPC::SetSecAU(Int_t sec)
2330 //----------------------------------------------------
2331 // Activate/deactivate selection for upper sectors
2332 //---------------------------------------------------
2334 //-----------------------------------------------------------------
2335 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2336 //-----------------------------------------------------------------
2341 //_____________________________________________________________________________
2342 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2344 //----------------------------------------
2345 // Select active lower sectors
2346 //----------------------------------------
2348 //-----------------------------------------------------------------
2349 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2350 //-----------------------------------------------------------------
2360 //_____________________________________________________________________________
2361 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2362 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2363 Int_t s11 , Int_t s12)
2365 //--------------------------------
2366 // Select active upper sectors
2367 //--------------------------------
2369 //-----------------------------------------------------------------
2370 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2371 //-----------------------------------------------------------------
2387 //_____________________________________________________________________________
2388 void AliTPC::SetSens(Int_t sens)
2391 //-------------------------------------------------------------
2392 // Activates/deactivates the sensitive strips at the center of
2393 // the pad row -- this is for the space-point resolution calculations
2394 //-------------------------------------------------------------
2396 //-----------------------------------------------------------------
2397 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2398 //-----------------------------------------------------------------
2403 void AliTPC::SetSide(Float_t side)
2408 //____________________________________________________________________________
2409 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2410 Float_t p2,Float_t p3)
2425 //_____________________________________________________________________________
2426 void AliTPC::Streamer(TBuffer &R__b)
2429 // Stream an object of class AliTPC.
2431 if (R__b.IsReading()) {
2432 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2433 AliDetector::Streamer(R__b);
2434 if (R__v < 2) return;
2438 fDigitsIndex = new Int_t[fNsectors+1];
2440 R__b.WriteVersion(AliTPC::IsA());
2441 AliDetector::Streamer(R__b);
2448 ClassImp(AliTPCcluster)
2450 //_____________________________________________________________________________
2451 AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
2454 // Creates a simulated cluster for the TPC
2456 fTracks[0] = lab[0];
2457 fTracks[1] = lab[1];
2458 fTracks[2] = lab[2];
2468 //_____________________________________________________________________________
2469 void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const
2472 // Transformation from local to global coordinate system
2474 x[0]=par->GetPadRowRadii(fSector,fPadRow);
2477 par->CRXYZtoXYZ(x,fSector,fPadRow,1);
2481 //_____________________________________________________________________________
2482 Int_t AliTPCcluster::Compare(TObject * o)
2485 // compare two clusters according y coordinata
2487 AliTPCcluster *cl= (AliTPCcluster *)o;
2488 if (fY<cl->fY) return -1;
2489 if (fY==cl->fY) return 0;
2493 Bool_t AliTPCcluster::IsSortable() const
2496 //make AliTPCcluster sortabale
2503 ClassImp(AliTPCdigit)
2505 //_____________________________________________________________________________
2506 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2510 // Creates a TPC digit object
2512 fSector = digits[0];
2513 fPadRow = digits[1];
2516 fSignal = digits[4];
2522 //_____________________________________________________________________________
2523 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2527 // Creates a TPC hit object
2538 ClassImp(AliTPCtrack)
2540 //_____________________________________________________________________________
2541 AliTPCtrack::AliTPCtrack(Float_t *hits)
2544 // Default creator for a TPC reconstructed track object
2546 fX=hits[0]; // This is dummy code !
2549 AliTPCtrack::AliTPCtrack(const AliTPCcluster *c,const TVector& xx,
2550 const TMatrix& CC, Double_t xref, Double_t alpha):
2551 x(xx),C(CC),fClusters(200)
2553 //-----------------------------------------------------------------
2554 // This is the main track constructor.
2556 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2557 //-----------------------------------------------------------------
2561 fClusters.AddLast((AliTPCcluster*)(c));
2564 //_____________________________________________________________________________
2565 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
2566 fClusters(t.fClusters.GetEntriesFast())
2568 //-----------------------------------------------------------------
2569 // This is a track copy constructor.
2571 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2572 //-----------------------------------------------------------------
2576 int n=t.fClusters.GetEntriesFast();
2577 for (int i=0; i<n; i++) fClusters.AddLast(t.fClusters.UncheckedAt(i));
2580 //_____________________________________________________________________________
2581 Int_t AliTPCtrack::Compare(TObject *o) {
2582 //-----------------------------------------------------------------
2583 // This function compares tracks according to their curvature.
2585 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2586 //-----------------------------------------------------------------
2587 AliTPCtrack *t=(AliTPCtrack*)o;
2588 Double_t co=t->GetSigmaY2();
2589 Double_t c =GetSigmaY2();
2591 else if (c<co) return -1;
2595 //_____________________________________________________________________________
2596 int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
2598 //-----------------------------------------------------------------
2599 // This function propagates a track to a reference plane x=xk.
2601 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2602 //-----------------------------------------------------------------
2603 if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
2604 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
2608 Double_t x1=fX, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
2609 Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
2610 Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
2612 x(0) += dx*(c1+c2)/(r1+r2);
2613 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2615 TMatrix F(5,5); F.UnitMatrix();
2616 Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
2617 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2618 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2619 Double_t cr=c1*r2+c2*r1;
2620 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2621 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2623 TMatrix tmp(F,TMatrix::kMult,C);
2624 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2628 //Multiple scattering******************
2629 Double_t ey=x(2)*fX - x(3);
2630 Double_t ex=sqrt(1-ey*ey);
2632 TMatrix Q(5,5); Q=0.;
2633 Q(2,2)=ez*ez+ey*ey; Q(2,3)=-ex*ey; Q(2,4)=-ex*ez;
2634 Q(3,2)=Q(2,3); Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
2635 Q(4,2)=Q(2,4); Q(4,3)= Q(3,4); Q(4,4)=1.;
2638 F(2,2)=-x(2)*ex; F(2,3)=-x(2)*ey;
2639 F(3,2)=-ex*(x(2)*fX-ey); F(3,3)=-(1.+ x(2)*fX*ey - ey*ey);
2640 F(4,2)=-ez*ex; F(4,3)=-ez*ey; F(4,4)=1.;
2643 Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2645 Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
2646 Double_t beta2=p2/(p2 + pm*pm);
2647 Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
2649 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
2653 //Energy losses************************
2654 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
2655 if (x1 < x2) dE=-dE;
2656 x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2657 //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2659 x1=fX; x2=xk; y1=x(0); z1=x(1);
2660 c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1);
2661 c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2);
2663 x(0) += dx*(c1+c2)/(r1+r2);
2664 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2667 rr=r1+r2; cc=c1+c2; xx=x1+x2;
2668 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2669 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2671 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2672 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2675 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2682 //_____________________________________________________________________________
2683 void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
2685 //-----------------------------------------------------------------
2686 // This function propagates tracks to the "vertex".
2688 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2689 //-----------------------------------------------------------------
2690 Double_t c=x(2)*fX - x(3);
2691 Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
2692 Double_t snf=tgf/sqrt(1.+ tgf*tgf);
2693 Double_t xv=(x(3)+snf)/x(2);
2694 PropagateTo(xv,x0,rho,pm);
2697 //_____________________________________________________________________________
2698 void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
2700 //-----------------------------------------------------------------
2701 // This function associates a clusters with this track.
2703 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2704 //-----------------------------------------------------------------
2705 TMatrix H(2,5); H.UnitMatrix();
2706 TMatrix Ht(TMatrix::kTransposed,H);
2707 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
2708 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2710 TMatrix tmp(H,TMatrix::kMult,C);
2711 TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
2713 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2714 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2715 R(1,0)*=-1; R(0,1)=R(1,0);
2720 TMatrix K(C,TMatrix::kMult,Ht); K*=R;
2723 x*=H; x-=m; x*=-1; x*=K; x+=savex;
2724 if (TMath::Abs(x(2)*fX-x(3)) >= 0.999) {
2725 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
2731 C.Mult(K,tmp); C-=saveC; C*=-1;
2733 fClusters.AddLast((AliTPCcluster*)c);
2737 //_____________________________________________________________________________
2738 int AliTPCtrack::Rotate(Double_t alpha)
2740 //-----------------------------------------------------------------
2741 // This function rotates this track.
2743 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2744 //-----------------------------------------------------------------
2747 Double_t x1=fX, y1=x(0);
2748 Double_t ca=cos(alpha), sa=sin(alpha);
2749 Double_t r1=x(2)*fX - x(3);
2752 x(0)=-x1*sa + y1*ca;
2753 x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
2755 Double_t r2=x(2)*fX - x(3);
2756 if (TMath::Abs(r2) >= 0.999) {
2757 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
2761 Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
2762 if ((x(0)-y0)*x(2) >= 0.) {
2763 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
2767 TMatrix F(5,5); F.UnitMatrix();
2770 F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa;
2771 F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
2772 TMatrix tmp(F,TMatrix::kMult,C);
2773 // Double_t dy2=C(0,0);
2774 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2775 // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
2776 // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
2781 //_____________________________________________________________________________
2782 void AliTPCtrack::UseClusters() const
2784 //-----------------------------------------------------------------
2785 // This function marks clusters associated with this track.
2787 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2788 //-----------------------------------------------------------------
2789 int num_of_clusters=fClusters.GetEntriesFast();
2790 for (int i=0; i<num_of_clusters; i++) {
2791 //if (i<=14) continue;
2792 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
2797 //_____________________________________________________________________________
2798 Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const
2800 //-----------------------------------------------------------------
2801 // This function calculates a predicted chi2 increment.
2803 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2804 //-----------------------------------------------------------------
2805 TMatrix H(2,5); H.UnitMatrix();
2806 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
2807 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2808 TVector res=x; res*=H; res-=m; //res*=-1;
2809 TMatrix tmp(H,TMatrix::kMult,C);
2810 TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
2812 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2813 if (TMath::Abs(det) < 1.e-10) {
2814 if (*this>4) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
2817 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2818 R(1,0)*=-1; R(0,1)=R(1,0);
2828 //_____________________________________________________________________________
2829 struct S { int lab; int max; };
2830 int AliTPCtrack::GetLabel(int nrows) const
2832 //-----------------------------------------------------------------
2833 // This function returns the track label. If label<0, this track is fake.
2835 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2836 //-----------------------------------------------------------------
2837 int num_of_clusters=fClusters.GetEntriesFast();
2838 S *s=new S[num_of_clusters];
2840 for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
2843 for (i=0; i<num_of_clusters; i++) {
2844 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
2845 lab=TMath::Abs(c->fTracks[0]);
2847 for (j=0; j<num_of_clusters; j++)
2848 if (s[j].lab==lab || s[j].max==0) break;
2854 for (i=0; i<num_of_clusters; i++)
2855 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
2859 for (i=0; i<num_of_clusters; i++) {
2860 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
2861 if (TMath::Abs(c->fTracks[1]) == lab ||
2862 TMath::Abs(c->fTracks[2]) == lab ) max++;
2865 if (1.-float(max)/num_of_clusters > 0.10) return -lab;
2867 int tail=int(0.08*nrows);
2868 if (num_of_clusters < tail) return lab;
2871 for (i=1; i<=tail; i++) {
2872 AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(num_of_clusters-i);
2873 if (lab == TMath::Abs(c->fTracks[0]) ||
2874 lab == TMath::Abs(c->fTracks[1]) ||
2875 lab == TMath::Abs(c->fTracks[2])) max++;
2877 if (max < int(0.5*tail)) return -lab;
2882 //_____________________________________________________________________________
2883 void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
2885 //-----------------------------------------------------------------
2886 // This function returns reconstructed track momentum in the global system.
2888 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2889 //-----------------------------------------------------------------
2890 Double_t pt=TMath::Abs(GetPt()); // GeV/c
2891 Double_t r=x(2)*fX-x(3);
2892 Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2893 px=-pt*(x(0)-y0)*x(2); //cos(phi);
2894 py=-pt*(x(3)-fX*x(2)); //sin(phi);
2896 Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2897 py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2901 //_____________________________________________________________________________
2902 Double_t AliTPCtrack::GetdEdX(Double_t low, Double_t up) const {
2903 //-----------------------------------------------------------------
2904 // This funtion calculates dE/dX within the "low" and "up" cuts.
2906 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2907 //-----------------------------------------------------------------
2908 int ncl=fClusters.GetEntriesFast();
2910 Double_t *q=new Double_t[ncl];
2912 for (i=0; i<ncl; i++) {
2913 AliTPCcluster *cl=(AliTPCcluster*)(fClusters.UncheckedAt(i));
2914 // if (cl->fdEdX > 3000) continue;
2915 if (cl->fdEdX <= 0) continue;
2923 for (i=0; i<n-1; i++) {
2924 if (q[i]<=q[i+1]) continue;
2925 Double_t tmp=q[i]; q[i]=q[i+1]; q[i+1]=tmp;
2930 int nl=int(low*n), nu=int(up *n);
2932 for (i=nl; i<=nu; i++) dedx += q[i];
2937 //_________________________________________________________________________
2939 // Classes for internal tracking use
2940 //_________________________________________________________________________
2941 void AliTPCRow::InsertCluster(const AliTPCcluster* c) {
2942 //-----------------------------------------------------------------------
2943 // Insert a cluster into this pad row in accordence with its y-coordinate
2945 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2946 //-----------------------------------------------------------------------
2947 if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2948 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2950 if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2952 memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2953 clusters[i]=c; num_of_clusters++;
2956 int AliTPCRow::Find(Double_t y) const {
2957 //-----------------------------------------------------------------------
2958 // Return the index of the nearest cluster
2960 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2961 //-----------------------------------------------------------------------
2962 if (y <= clusters[0]->fY) return 0;
2963 if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2964 int b=0, e=num_of_clusters-1, m=(b+e)/2;
2965 for (; b<e; m=(b+e)/2) {
2966 if (y > clusters[m]->fY) b=m+1;