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.19.2.4 2000/06/26 07:39:42 kowal2
19 Changes to obey the coding rules
21 Revision 1.19.2.3 2000/06/25 08:38:41 kowal2
22 Splitted from AliTPCtracking
24 Revision 1.19.2.2 2000/06/16 12:59:28 kowal2
25 Changed parameter settings
27 Revision 1.19.2.1 2000/06/09 07:15:07 kowal2
29 Defaults loaded automatically (hard-wired)
30 Optional parameters can be set via macro called in the constructor
32 Revision 1.19 2000/04/18 19:00:59 fca
33 Small bug fixes to TPC files
35 Revision 1.18 2000/04/17 09:37:33 kowal2
36 removed obsolete AliTPCDigitsDisplay.C
38 Revision 1.17.2.2 2000/04/10 08:15:12 kowal2
40 New, experimental data structure from M. Ivanov
41 New tracking algorithm
42 Different pad geometry for different sectors
43 Digitization rewritten
45 Revision 1.17.2.1 2000/04/10 07:56:53 kowal2
46 Not used anymore - removed
48 Revision 1.17 2000/01/19 17:17:30 fca
49 Introducing a list of lists of hits -- more hits allowed for detector now
51 Revision 1.16 1999/11/05 09:29:23 fca
52 Accept only signals > 0
54 Revision 1.15 1999/10/08 06:26:53 fca
55 Removed ClustersIndex - not used anymore
57 Revision 1.14 1999/09/29 09:24:33 fca
58 Introduction of the Copyright and cvs Log
62 ///////////////////////////////////////////////////////////////////////////////
64 // Time Projection Chamber //
65 // This class contains the basic functions for the Time Projection Chamber //
66 // detector. Functions specific to one particular geometry are //
67 // contained in the derived classes //
71 <img src="picts/AliTPCClass.gif">
76 ///////////////////////////////////////////////////////////////////////////////
84 #include <TGeometry.h>
87 #include <TObjectTable.h>
88 #include "TParticle.h"
97 #include "AliTPCParamSR.h"
98 #include "AliTPCPRF2D.h"
99 #include "AliTPCRF1D.h"
100 #include "AliDigits.h"
101 #include "AliSimDigits.h"
103 #include "AliTPCDigitsArray.h"
104 #include "AliCluster.h"
105 #include "AliClusters.h"
106 #include "AliTPCClustersRow.h"
107 #include "AliTPCClustersArray.h"
109 #include "AliTPCcluster.h"
110 #include "AliTPCclusterer.h"
111 #include "AliTPCtracker.h"
113 #include <TInterpreter.h>
117 //_____________________________________________________________________________
121 // Default constructor
133 //_____________________________________________________________________________
134 AliTPC::AliTPC(const char *name, const char *title)
135 : AliDetector(name,title)
138 // Standard constructor
142 // Initialise arrays of hits and digits
143 fHits = new TClonesArray("AliTPChit", 176);
144 gAlice->AddHitList(fHits);
149 // Initialise counters
155 // Initialise color attributes
156 SetMarkerColor(kYellow);
159 // Set TPC parameters
162 if (!strcmp(title,"Default")) {
163 fTPCParam = new AliTPCParamSR;
165 cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
171 //_____________________________________________________________________________
184 //_____________________________________________________________________________
185 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
188 // Add a hit to the list
190 TClonesArray &lhits = *fHits;
191 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
194 //_____________________________________________________________________________
195 void AliTPC::BuildGeometry()
199 // Build TPC ROOT TNode geometry for the event display
204 const int kColorTPC=19;
205 char name[5], title[25];
206 const Double_t kDegrad=TMath::Pi()/180;
207 const Double_t kRaddeg=180./TMath::Pi();
210 Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
211 Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
213 Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
214 Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
216 Int_t nLo = fTPCParam->GetNInnerSector()/2;
217 Int_t nHi = fTPCParam->GetNOuterSector()/2;
219 const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
220 const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
221 const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
222 const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
225 const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
226 const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
232 // Get ALICE top node
235 nTop=gAlice->GetGeometry()->GetNode("alice");
239 rl = fTPCParam->GetInnerRadiusLow();
240 ru = fTPCParam->GetInnerRadiusUp();
244 sprintf(name,"LS%2.2d",i);
246 sprintf(title,"TPC low sector %3d",i);
249 tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
250 kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
251 tubs->SetNumberOfDivisions(1);
253 nNode = new TNode(name,title,name,0,0,0,"");
254 nNode->SetLineColor(kColorTPC);
260 rl = fTPCParam->GetOuterRadiusLow();
261 ru = fTPCParam->GetOuterRadiusUp();
264 sprintf(name,"US%2.2d",i);
266 sprintf(title,"TPC upper sector %d",i);
268 tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
269 khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
270 tubs->SetNumberOfDivisions(1);
272 nNode = new TNode(name,title,name,0,0,0,"");
273 nNode->SetLineColor(kColorTPC);
279 //_____________________________________________________________________________
280 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
283 // Calculate distance from TPC to mouse on the display
289 void AliTPC::Clusters2Tracks(TFile *of) {
290 //-----------------------------------------------------------------
291 // This is a track finder.
292 //-----------------------------------------------------------------
293 AliTPCtracker::Clusters2Tracks(fTPCParam,of);
296 //_____________________________________________________________________________
297 void AliTPC::CreateMaterials()
299 //-----------------------------------------------
300 // Create Materials for for TPC simulations
301 //-----------------------------------------------
303 //-----------------------------------------------------------------
304 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
305 //-----------------------------------------------------------------
307 Int_t iSXFLD=gAlice->Field()->Integ();
308 Float_t sXMGMX=gAlice->Field()->Max();
310 Float_t amat[5]; // atomic numbers
311 Float_t zmat[5]; // z
312 Float_t wmat[5]; // proportions
318 //***************** Gases *************************
320 //-------------------------------------------------
322 //-------------------------------------------------
333 AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
343 AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
346 //--------------------------------------------------------------
348 //--------------------------------------------------------------
365 amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
367 AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
382 amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
384 AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
400 amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
402 AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
404 //----------------------------------------------------------------
405 // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
406 //----------------------------------------------------------------
412 Float_t rho,absl,X0,buf[1];
416 for(nc = 0;nc<fNoComp;nc++)
419 // retrive material constants
421 gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
426 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
428 am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
429 density += fMixtProp[nc]*rho; // density of the mixture
433 // mixture proportions by weight!
435 for(nc = 0;nc<fNoComp;nc++)
438 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
440 wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
441 apure[nnc] : amol[nnc])/am;
445 // Drift gases 1 - nonsensitive, 2 - sensitive
447 AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
448 AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
457 AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
460 //----------------------------------------------------------------------
462 //----------------------------------------------------------------------
484 AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
506 AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
524 AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
542 AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
544 // G10 60% SiO2 + 40% epoxy, I use A and Z for SiO2
557 AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz
559 gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf);
561 AliMaterial(39,"G10",amat[0],zmat[0],density,999.,999.);
570 AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
579 AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
588 AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
606 AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
625 AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
629 //----------------------------------------------------------
630 // tracking media for gases
631 //----------------------------------------------------------
633 AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
634 AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
635 AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
636 AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
638 //-----------------------------------------------------------
639 // tracking media for solids
640 //-----------------------------------------------------------
642 AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
643 AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
644 AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
645 AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
646 AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
647 AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
648 AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
649 AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
650 AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
651 AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
656 void AliTPC::Digits2Clusters(TFile *of)
658 //-----------------------------------------------------------------
659 // This is a simple cluster finder.
660 //-----------------------------------------------------------------
661 AliTPCclusterer::Digits2Clusters(fTPCParam,of);
664 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
665 extern Double_t SigmaZ2(Double_t, Double_t);
666 //_____________________________________________________________________________
667 void AliTPC::Hits2Clusters(TFile *of)
669 //--------------------------------------------------------
670 // TPC simple cluster generator from hits
671 // obtained from the TPC Fast Simulator
672 // The point errors are taken from the parametrization
673 //--------------------------------------------------------
675 //-----------------------------------------------------------------
676 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
677 //-----------------------------------------------------------------
678 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
679 // Jouri.Belikov@cern.ch
680 //----------------------------------------------------------------
682 /////////////////////////////////////////////////////////////////////////////
684 //---------------------------------------------------------------------
685 // ALICE TPC Cluster Parameters
686 //--------------------------------------------------------------------
690 // Cluster width in rphi
691 const Float_t kACrphi=0.18322;
692 const Float_t kBCrphi=0.59551e-3;
693 const Float_t kCCrphi=0.60952e-1;
694 // Cluster width in z
695 const Float_t kACz=0.19081;
696 const Float_t kBCz=0.55938e-3;
697 const Float_t kCCz=0.30428;
699 TDirectory *savedir=gDirectory;
702 cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
707 printf("AliTPCParam MUST be created firstly\n");
711 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
713 TParticle *particle; // pointer to a given particle
714 AliTPChit *tpcHit; // pointer to a sigle TPC hit
715 TClonesArray *particles; //pointer to the particle list
719 Float_t pl,pt,tanth,rpad,ratio;
722 //---------------------------------------------------------------
723 // Get the access to the tracks
724 //---------------------------------------------------------------
726 TTree *tH = gAlice->TreeH();
727 Stat_t ntracks = tH->GetEntries();
728 particles=gAlice->Particles();
730 //Switch to the output file
733 fTPCParam->Write(fTPCParam->GetTitle());
734 AliTPCClustersArray carray;
735 carray.Setup(fTPCParam);
736 carray.SetClusterType("AliTPCcluster");
739 Int_t nclusters=0; //cluster counter
741 //------------------------------------------------------------
742 // Loop over all sectors (72 sectors for 20 deg
743 // segmentation for both lower and upper sectors)
744 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
745 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
747 // First cluster for sector 0 starts at "0"
748 //------------------------------------------------------------
750 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
752 fTPCParam->AdjustCosSin(isec,cph,sph);
754 //------------------------------------------------------------
756 //------------------------------------------------------------
758 for(Int_t track=0;track<ntracks;track++){
762 // Get number of the TPC hits
764 nhits=fHits->GetEntriesFast();
768 for(Int_t hit=0;hit<nhits;hit++){
769 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
770 if (tpcHit->fQ == 0.) continue; //information about track (I.Belikov)
771 sector=tpcHit->fSector; // sector number
772 if(sector != isec) continue; //terminate iteration
773 ipart=tpcHit->fTrack;
774 particle=(TParticle*)particles->UncheckedAt(ipart);
777 if(pt < 1.e-9) pt=1.e-9;
779 tanth = TMath::Abs(tanth);
780 rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
781 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
783 // space-point resolutions
785 sigmaRphi=SigmaY2(rpad,tanth,pt);
786 sigmaZ =SigmaZ2(rpad,tanth );
790 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
791 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
793 // temporary protection
795 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
796 if(sigmaZ < 0.) sigmaZ=0.4e-3;
797 if(clRphi < 0.) clRphi=2.5e-3;
798 if(clZ < 0.) clZ=2.5e-5;
803 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
804 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
806 Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
807 Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
808 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
809 Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
810 fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
811 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
812 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
813 xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigmaZ)); // z
814 if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->fZ;
815 xyz[2]=tpcHit->fQ; // q
816 xyz[3]=sigmaRphi; // fSigmaY2
817 xyz[4]=sigmaZ; // fSigmaZ2
819 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
820 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
822 Int_t tracks[3]={tpcHit->fTrack, -1, -1};
823 AliTPCcluster cluster(xyz,tracks);
825 clrow->InsertCluster(&cluster); nclusters++;
827 } // end of loop over hits
829 } // end of loop over tracks
831 Int_t nrows=fTPCParam->GetNRow(isec);
832 for (Int_t irow=0; irow<nrows; irow++) {
833 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
834 if (!clrow) continue;
835 carray.StoreRow(isec,irow);
836 carray.ClearRow(isec,irow);
839 } // end of loop over sectors
841 cerr<<"Number of made clusters : "<<nclusters<<" \n";
843 carray.GetTree()->Write();
845 savedir->cd(); //switch back to the input file
849 //_________________________________________________________________
850 void AliTPC::Hits2ExactClustersSector(Int_t isec)
852 //--------------------------------------------------------
853 //calculate exact cross point of track and given pad row
854 //resulting values are expressed in "digit" coordinata
855 //--------------------------------------------------------
857 //-----------------------------------------------------------------
858 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
859 //-----------------------------------------------------------------
861 if (fClustersArray==0){
865 TParticle *particle; // pointer to a given particle
866 AliTPChit *tpcHit; // pointer to a sigle TPC hit
867 TClonesArray *particles; //pointer to the particle list
870 const Int_t kcmaxhits=30000;
871 TVector * xxxx = new TVector(kcmaxhits*4);
872 TVector & xxx = *xxxx;
873 Int_t maxhits = kcmaxhits;
874 //construct array for each padrow
875 for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
876 fClustersArray->CreateRow(isec,i);
878 //---------------------------------------------------------------
879 // Get the access to the tracks
880 //---------------------------------------------------------------
882 TTree *tH = gAlice->TreeH();
883 Stat_t ntracks = tH->GetEntries();
884 particles=gAlice->Particles();
885 Int_t npart = particles->GetEntriesFast();
887 //------------------------------------------------------------
889 //------------------------------------------------------------
891 for(Int_t track=0;track<ntracks;track++){
895 // Get number of the TPC hits and a pointer
898 nhits=fHits->GetEntriesFast();
902 Int_t currentIndex=0;
903 Int_t lastrow=-1; //last writen row
904 for(Int_t hit=0;hit<nhits;hit++){
905 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
906 if (tpcHit==0) continue;
907 sector=tpcHit->fSector; // sector number
908 if(sector != isec) continue;
909 ipart=tpcHit->fTrack;
910 if (ipart<npart) particle=(TParticle*)particles->UncheckedAt(ipart);
914 Float_t x[3]={tpcHit->fX,tpcHit->fY,tpcHit->fZ};
915 Int_t index[3]={1,isec,0};
916 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
917 if (currentrow<0) continue;
918 if (lastrow<0) lastrow=currentrow;
919 if (currentrow==lastrow){
920 if ( currentIndex>=maxhits){
922 xxx.ResizeTo(4*maxhits);
924 xxx(currentIndex*4)=x[0];
925 xxx(currentIndex*4+1)=x[1];
926 xxx(currentIndex*4+2)=x[2];
927 xxx(currentIndex*4+3)=tpcHit->fQ;
943 for (Int_t index=0;index<currentIndex;index++){
945 x=x2=x3=x4=xxx(index*4);
953 sumy+=xxx(index*4+1);
954 sumxy+=xxx(index*4+1)*x;
955 sumx2y+=xxx(index*4+1)*x2;
956 sumz+=xxx(index*4+2);
957 sumxz+=xxx(index*4+2)*x;
958 sumx2z+=xxx(index*4+2)*x2;
959 sumq+=xxx(index*4+3);
961 Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
962 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
963 sumx2*(sumx*sumx3-sumx2*sumx2);
965 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
966 sumx2*(sumxy*sumx3-sumx2y*sumx2);
967 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
968 sumx2*(sumxz*sumx3-sumx2z*sumx2);
970 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
971 sumx2*(sumx*sumx2y-sumx2*sumxy);
972 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
973 sumx2*(sumx*sumx2z-sumx2*sumxz);
975 Float_t y=detay/det+centralPad;
977 Float_t by=detby/det; //y angle
978 Float_t bz=detbz/det; //z angle
979 sumy/=Float_t(currentIndex);
980 sumz/=Float_t(currentIndex);
989 AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
990 if (row!=0) row->InsertCluster(&cl);
993 } //end of calculating cluster for given row
997 } // end of loop over hits
998 } // end of loop over tracks
999 //write padrows to tree
1000 for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1001 fClustersArray->StoreRow(isec,ii);
1002 fClustersArray->ClearRow(isec,ii);
1008 //__________________________________________________________________
1009 void AliTPC::Hits2Digits()
1011 //----------------------------------------------------
1012 // Loop over all sectors
1013 //----------------------------------------------------
1016 printf("AliTPCParam MUST be created firstly\n");
1020 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1025 //_____________________________________________________________________________
1026 void AliTPC::Hits2DigitsSector(Int_t isec)
1028 //-------------------------------------------------------------------
1029 // TPC conversion from hits to digits.
1030 //-------------------------------------------------------------------
1032 //-----------------------------------------------------------------
1033 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1034 //-----------------------------------------------------------------
1036 //-------------------------------------------------------
1037 // Get the access to the track hits
1038 //-------------------------------------------------------
1041 TTree *tH = gAlice->TreeH(); // pointer to the hits tree
1042 Stat_t ntracks = tH->GetEntries();
1046 //-------------------------------------------
1047 // Only if there are any tracks...
1048 //-------------------------------------------
1052 printf("*** Processing sector number %d ***\n",isec);
1054 Int_t nrows =fTPCParam->GetNRow(isec);
1056 row= new TObjArray* [nrows];
1058 MakeSector(isec,nrows,tH,ntracks,row);
1060 //--------------------------------------------------------
1061 // Digitize this sector, row by row
1062 // row[i] is the pointer to the TObjArray of TVectors,
1063 // each one containing electrons accepted on this
1064 // row, assigned into tracks
1065 //--------------------------------------------------------
1069 if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree();
1071 for (i=0;i<nrows;i++){
1073 AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1075 DigitizeRow(i,isec,row);
1077 fDigitsArray->StoreRow(isec,i);
1079 Int_t ndig = dig->GetDigitSize();
1081 printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1083 fDigitsArray->ClearRow(isec,i);
1086 } // end of the sector digitization
1088 for(i=0;i<nrows;i++){
1092 delete [] row; // delete the array of pointers to TObjArray-s
1096 } // end of Hits2DigitsSector
1099 //_____________________________________________________________________________
1100 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1102 //-----------------------------------------------------------
1103 // Single row digitization, coupling from the neighbouring
1104 // rows taken into account
1105 //-----------------------------------------------------------
1107 //-----------------------------------------------------------------
1108 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1109 // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1110 //-----------------------------------------------------------------
1113 Float_t zerosup = fTPCParam->GetZeroSup();
1114 Int_t nrows =fTPCParam->GetNRow(isec);
1115 fCurrentIndex[1]= isec;
1118 Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1119 Int_t nofTbins = fTPCParam->GetMaxTBin();
1120 Int_t indexRange[4];
1122 // Integrated signal for this row
1123 // and a single track signal
1125 TMatrix *m1 = new TMatrix(0,nofPads,0,nofTbins); // integrated
1126 TMatrix *m2 = new TMatrix(0,nofPads,0,nofTbins); // single
1128 TMatrix &total = *m1;
1130 // Array of pointers to the label-signal list
1132 Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1133 Float_t **pList = new Float_t* [nofDigits];
1137 for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1141 Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1142 Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1143 for (Int_t row= row1;row<=row2;row++){
1144 Int_t nTracks= rows[row]->GetEntries();
1145 for (i1=0;i1<nTracks;i1++){
1146 fCurrentIndex[2]= row;
1147 fCurrentIndex[3]=irow;
1149 m2->Zero(); // clear single track signal matrix
1150 Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1151 GetList(trackLabel,nofPads,m2,indexRange,pList);
1153 else GetSignal(rows[row],i1,0,m1,indexRange);
1159 AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1160 for(Int_t ip=0;ip<nofPads;ip++){
1161 for(Int_t it=0;it<nofTbins;it++){
1163 Float_t q = total(ip,it);
1165 Int_t gi =it*nofPads+ip; // global index
1167 q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
1171 if(q <=zerosup) continue; // do not fill zeros
1172 if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation
1175 // "real" signal or electronic noise (list = -1)?
1178 for(Int_t j1=0;j1<3;j1++){
1179 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1184 <A NAME="AliDigits"></A>
1185 using of AliDigits object
1188 dig->SetDigitFast((Short_t)q,it,ip);
1189 if (fDigitsArray->IsSimulated())
1191 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1192 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1193 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1197 } // end of loop over time buckets
1198 } // end of lop over pads
1201 // This row has been digitized, delete nonused stuff
1204 for(lp=0;lp<nofDigits;lp++){
1205 if(pList[lp]) delete [] pList[lp];
1214 } // end of DigitizeRow
1216 //_____________________________________________________________________________
1218 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
1222 //---------------------------------------------------------------
1223 // Calculates 2-D signal (pad,time) for a single track,
1224 // returns a pointer to the signal matrix and the track label
1225 // No digitization is performed at this level!!!
1226 //---------------------------------------------------------------
1228 //-----------------------------------------------------------------
1229 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1230 // Modified: Marian Ivanov
1231 //-----------------------------------------------------------------
1235 tv = (TVector*)p1->At(ntr); // pointer to a track
1238 Float_t label = v(0);
1239 Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
1241 Int_t nElectrons = (tv->GetNrows()-1)/4;
1242 indexRange[0]=9999; // min pad
1243 indexRange[1]=-1; // max pad
1244 indexRange[2]=9999; //min time
1245 indexRange[3]=-1; // max time
1247 // Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1249 TMatrix &signal = *m1;
1250 TMatrix &total = *m2;
1252 // Loop over all electrons
1254 for(Int_t nel=0; nel<nElectrons; nel++){
1256 Float_t aval = v(idx+4);
1257 Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1258 Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1259 Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
1261 if (n>0) for (Int_t i =0; i<n; i++){
1262 Int_t *index = fTPCParam->GetResBin(i);
1263 Int_t pad=index[1]+centralPad; //in digit coordinates central pad has coordinate 0
1264 if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1265 Int_t time=index[2];
1266 Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1267 weight *= eltoadcfac;
1269 if (m1!=0) signal(pad,time)+=weight;
1270 total(pad,time)+=weight;
1271 indexRange[0]=TMath::Min(indexRange[0],pad);
1272 indexRange[1]=TMath::Max(indexRange[1],pad);
1273 indexRange[2]=TMath::Min(indexRange[2],time);
1274 indexRange[3]=TMath::Max(indexRange[3],time);
1277 } // end of loop over electrons
1279 return label; // returns track label when finished
1282 //_____________________________________________________________________________
1283 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
1286 //----------------------------------------------------------------------
1287 // Updates the list of tracks contributing to digits for a given row
1288 //----------------------------------------------------------------------
1290 //-----------------------------------------------------------------
1291 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1292 //-----------------------------------------------------------------
1294 TMatrix &signal = *m;
1296 // lop over nonzero digits
1298 for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1299 for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1302 // accept only the contribution larger than 500 electrons (1/2 s_noise)
1304 if(signal(ip,it)<0.5) continue;
1307 Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1309 if(!pList[globalIndex]){
1312 // Create new list (6 elements - 3 signals and 3 labels),
1315 pList[globalIndex] = new Float_t [6];
1319 *pList[globalIndex] = -1.;
1320 *(pList[globalIndex]+1) = -1.;
1321 *(pList[globalIndex]+2) = -1.;
1322 *(pList[globalIndex]+3) = -1.;
1323 *(pList[globalIndex]+4) = -1.;
1324 *(pList[globalIndex]+5) = -1.;
1327 *pList[globalIndex] = label;
1328 *(pList[globalIndex]+3) = signal(ip,it);
1332 // check the signal magnitude
1334 Float_t highest = *(pList[globalIndex]+3);
1335 Float_t middle = *(pList[globalIndex]+4);
1336 Float_t lowest = *(pList[globalIndex]+5);
1339 // compare the new signal with already existing list
1342 if(signal(ip,it)<lowest) continue; // neglect this track
1346 if (signal(ip,it)>highest){
1347 *(pList[globalIndex]+5) = middle;
1348 *(pList[globalIndex]+4) = highest;
1349 *(pList[globalIndex]+3) = signal(ip,it);
1351 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1352 *(pList[globalIndex]+1) = *pList[globalIndex];
1353 *pList[globalIndex] = label;
1355 else if (signal(ip,it)>middle){
1356 *(pList[globalIndex]+5) = middle;
1357 *(pList[globalIndex]+4) = signal(ip,it);
1359 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1360 *(pList[globalIndex]+1) = label;
1363 *(pList[globalIndex]+5) = signal(ip,it);
1364 *(pList[globalIndex]+2) = label;
1368 } // end of loop over pads
1369 } // end of loop over time bins
1374 //___________________________________________________________________
1375 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1376 Stat_t ntracks,TObjArray **row)
1379 //-----------------------------------------------------------------
1380 // Prepares the sector digitization, creates the vectors of
1381 // tracks for each row of this sector. The track vector
1382 // contains the track label and the position of electrons.
1383 //-----------------------------------------------------------------
1385 //-----------------------------------------------------------------
1386 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1387 //-----------------------------------------------------------------
1389 Float_t gasgain = fTPCParam->GetGasGain();
1393 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1395 //----------------------------------------------
1396 // Create TObjArray-s, one for each row,
1397 // each TObjArray will store the TVectors
1398 // of electrons, one TVector per each track.
1399 //----------------------------------------------
1401 Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
1402 TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
1403 for(i=0; i<nrows; i++){
1404 row[i] = new TObjArray;
1411 //--------------------------------------------------------------------
1412 // Loop over tracks, the "track" contains the full history
1413 //--------------------------------------------------------------------
1415 Int_t previousTrack,currentTrack;
1416 previousTrack = -1; // nothing to store so far!
1418 for(Int_t track=0;track<ntracks;track++){
1422 TH->GetEvent(track); // get next track
1423 Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1425 if(nhits == 0) continue; // no hits in the TPC for this track
1427 //--------------------------------------------------------------
1429 //--------------------------------------------------------------
1431 for(Int_t hit=0;hit<nhits;hit++){
1433 tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1435 Int_t sector=tpcHit->fSector; // sector number
1436 if(sector != isec) continue;
1438 currentTrack = tpcHit->fTrack; // track number
1439 if(currentTrack != previousTrack){
1441 // store already filled fTrack
1443 for(i=0;i<nrows;i++){
1444 if(previousTrack != -1){
1445 if(nofElectrons[i]>0){
1446 TVector &v = *tracks[i];
1447 v(0) = previousTrack;
1448 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1449 row[i]->Add(tracks[i]);
1452 delete tracks[i]; // delete empty TVector
1458 tracks[i] = new TVector(481); // TVectors for the next fTrack
1460 } // end of loop over rows
1462 previousTrack=currentTrack; // update track label
1465 Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1467 //---------------------------------------------------
1468 // Calculate the electron attachment probability
1469 //---------------------------------------------------
1472 Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->fZ))
1473 /fTPCParam->GetDriftV();
1475 Float_t attProb = fTPCParam->GetAttCoef()*
1476 fTPCParam->GetOxyCont()*time; // fraction!
1478 //-----------------------------------------------
1479 // Loop over electrons
1480 //-----------------------------------------------
1483 for(Int_t nel=0;nel<qI;nel++){
1484 // skip if electron lost due to the attachment
1485 if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1489 xyz[3]= (Float_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1492 TransportElectron(xyz,index); //MI change -august
1494 fTPCParam->GetPadRow(xyz,index); //MI change august
1495 rowNumber = index[2];
1496 //transform position to local digit coordinates
1497 //relative to nearest pad row
1498 if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;
1499 nofElectrons[rowNumber]++;
1500 //----------------------------------
1501 // Expand vector if necessary
1502 //----------------------------------
1503 if(nofElectrons[rowNumber]>120){
1504 Int_t range = tracks[rowNumber]->GetNrows();
1505 if((nofElectrons[rowNumber])>(range-1)/4){
1507 tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
1511 TVector &v = *tracks[rowNumber];
1512 Int_t idx = 4*nofElectrons[rowNumber]-3;
1514 v(idx)= xyz[0]; // X - pad row coordinate
1515 v(idx+1)=xyz[1]; // Y - pad coordinate (along the pad-row)
1516 v(idx+2)=xyz[2]; // Z - time bin coordinate
1517 v(idx+3)=xyz[3]; // avalanche size
1518 } // end of loop over electrons
1520 } // end of loop over hits
1521 } // end of loop over tracks
1524 // store remaining track (the last one) if not empty
1527 for(i=0;i<nrows;i++){
1528 if(nofElectrons[i]>0){
1529 TVector &v = *tracks[i];
1530 v(0) = previousTrack;
1531 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1532 row[i]->Add(tracks[i]);
1541 delete [] nofElectrons;
1544 } // end of MakeSector
1547 //_____________________________________________________________________________
1551 // Initialise TPC detector after definition of geometry
1556 for(i=0;i<35;i++) printf("*");
1557 printf(" TPC_INIT ");
1558 for(i=0;i<35;i++) printf("*");
1561 for(i=0;i<80;i++) printf("*");
1565 //_____________________________________________________________________________
1566 void AliTPC::MakeBranch(Option_t* option)
1569 // Create Tree branches for the TPC.
1571 Int_t buffersize = 4000;
1572 char branchname[10];
1573 sprintf(branchname,"%s",GetName());
1575 AliDetector::MakeBranch(option);
1577 char *d = strstr(option,"D");
1579 if (fDigits && gAlice->TreeD() && d) {
1580 gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1581 printf("Making Branch %s for digits\n",branchname);
1585 //_____________________________________________________________________________
1586 void AliTPC::ResetDigits()
1589 // Reset number of digits and the digits array for this detector
1592 if (fDigits) fDigits->Clear();
1595 //_____________________________________________________________________________
1596 void AliTPC::SetSecAL(Int_t sec)
1598 //---------------------------------------------------
1599 // Activate/deactivate selection for lower sectors
1600 //---------------------------------------------------
1602 //-----------------------------------------------------------------
1603 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1604 //-----------------------------------------------------------------
1609 //_____________________________________________________________________________
1610 void AliTPC::SetSecAU(Int_t sec)
1612 //----------------------------------------------------
1613 // Activate/deactivate selection for upper sectors
1614 //---------------------------------------------------
1616 //-----------------------------------------------------------------
1617 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1618 //-----------------------------------------------------------------
1623 //_____________________________________________________________________________
1624 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1626 //----------------------------------------
1627 // Select active lower sectors
1628 //----------------------------------------
1630 //-----------------------------------------------------------------
1631 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1632 //-----------------------------------------------------------------
1642 //_____________________________________________________________________________
1643 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1644 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
1645 Int_t s11 , Int_t s12)
1647 //--------------------------------
1648 // Select active upper sectors
1649 //--------------------------------
1651 //-----------------------------------------------------------------
1652 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1653 //-----------------------------------------------------------------
1669 //_____________________________________________________________________________
1670 void AliTPC::SetSens(Int_t sens)
1673 //-------------------------------------------------------------
1674 // Activates/deactivates the sensitive strips at the center of
1675 // the pad row -- this is for the space-point resolution calculations
1676 //-------------------------------------------------------------
1678 //-----------------------------------------------------------------
1679 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1680 //-----------------------------------------------------------------
1685 void AliTPC::SetSide(Float_t side=0.)
1687 // choice of the TPC side
1692 //____________________________________________________________________________
1693 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
1694 Float_t p2,Float_t p3)
1697 // gax mixture definition
1711 //_____________________________________________________________________________
1713 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
1716 // electron transport taking into account:
1718 // 2.ExB at the wires
1719 // 3. nonisochronity
1721 // xyz and index must be already transformed to system 1
1724 fTPCParam->Transform1to2(xyz,index);
1727 Float_t driftl=xyz[2];
1728 if(driftl<0.01) driftl=0.01;
1729 driftl=TMath::Sqrt(driftl);
1730 Float_t sigT = driftl*(fTPCParam->GetDiffT());
1731 Float_t sigL = driftl*(fTPCParam->GetDiffL());
1732 xyz[0]=gRandom->Gaus(xyz[0],sigT);
1733 xyz[1]=gRandom->Gaus(xyz[1],sigT);
1734 xyz[2]=gRandom->Gaus(xyz[2],sigL);
1738 if (fTPCParam->GetMWPCReadout()==kTRUE){
1740 fTPCParam->Transform2to2NearestWire(xyz,index);
1741 Float_t dx=xyz[0]-x1;
1742 xyz[1]+=dx*(fTPCParam->GetOmegaTau());
1744 //add nonisochronity (not implemented yet)
1747 //_____________________________________________________________________________
1748 void AliTPC::Streamer(TBuffer &R__b)
1751 // Stream an object of class AliTPC.
1753 if (R__b.IsReading()) {
1754 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1755 AliDetector::Streamer(R__b);
1756 if (R__v < 2) return;
1759 R__b.WriteVersion(AliTPC::IsA());
1760 AliDetector::Streamer(R__b);
1765 ClassImp(AliTPCdigit)
1767 //_____________________________________________________________________________
1768 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
1772 // Creates a TPC digit object
1774 fSector = digits[0];
1775 fPadRow = digits[1];
1778 fSignal = digits[4];
1784 //_____________________________________________________________________________
1785 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1789 // Creates a TPC hit object