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 ///////////////////////////////////////////////////////////////////////////////
20 // Time Projection Chamber //
21 // This class contains the basic functions for the Time Projection Chamber //
22 // detector. Functions specific to one particular geometry are //
23 // contained in the derived classes //
27 <img src="picts/AliTPCClass.gif">
32 ///////////////////////////////////////////////////////////////////////////////
36 #include <Riostream.h>
40 #include <TGeometry.h>
41 #include <TInterpreter.h>
45 #include <TObjectTable.h>
46 #include <TParticle.h>
53 #include <TVirtualMC.h>
57 #include "AliArrayBranch.h"
58 #include "AliClusters.h"
59 #include "AliComplexCluster.h"
60 #include "AliDigits.h"
62 #include "AliPoints.h"
64 #include "AliRunLoader.h"
65 #include "AliSimDigits.h"
68 #include "AliTPCClustersArray.h"
69 #include "AliTPCClustersRow.h"
70 #include "AliTPCDigitsArray.h"
71 #include "AliTPCLoader.h"
72 #include "AliTPCPRF2D.h"
73 #include "AliTPCParamSR.h"
74 #include "AliTPCRF1D.h"
75 #include "AliTPCTrackHits.h"
76 #include "AliTPCTrackHitsV2.h"
77 #include "AliTPCcluster.h"
78 #include "AliTrackReference.h"
80 #include "AliTPCDigitizer.h"
81 #include "AliTPCclustererMI.h"
82 #include "AliTPCtrackerMI.h"
83 #include "AliTPCpidESD.h"
88 //_____________________________________________________________________________
89 // helper class for fast matrix and vector manipulation - no range checking
90 // origin - Marian Ivanov
92 class AliTPCFastMatrix : public TMatrix {
94 AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb);
95 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
96 Float_t & UncheckedAt(Int_t rown, Int_t coln) const {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces
97 Float_t UncheckedAtFast(Int_t rown, Int_t coln) const {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces
99 Float_t & UncheckedAt(Int_t rown, Int_t coln) const {return (fIndex[coln])[rown];} //fast acces
100 Float_t UncheckedAtFast(Int_t rown, Int_t coln) const {return (fIndex[coln])[rown];} //fast acces
104 AliTPCFastMatrix::AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb):
105 TMatrix(rowlwb, rowupb,collwb,colupb)
108 //_____________________________________________________________________________
109 class AliTPCFastVector : public TVector {
111 AliTPCFastVector(Int_t size);
112 virtual ~AliTPCFastVector(){;}
113 Float_t & UncheckedAt(Int_t index) const {return fElements[index];} //fast acces
117 AliTPCFastVector::AliTPCFastVector(Int_t size):TVector(size){
120 //_____________________________________________________________________________
124 // Default constructor
135 fHitType = 2; //default CONTAINERS - based on ROOT structure
142 //_____________________________________________________________________________
143 AliTPC::AliTPC(const char *name, const char *title)
144 : AliDetector(name,title)
147 // Standard constructor
151 // Initialise arrays of hits and digits
152 fHits = new TClonesArray("AliTPChit", 176);
153 gAlice->GetMCApp()->AddHitList(fHits);
158 fTrackHits = new AliTPCTrackHitsV2;
159 fTrackHits->SetHitPrecision(0.002);
160 fTrackHits->SetStepPrecision(0.003);
161 fTrackHits->SetMaxDistance(100);
163 fTrackHitsOld = new AliTPCTrackHits; //MI - 13.09.2000
164 fTrackHitsOld->SetHitPrecision(0.002);
165 fTrackHitsOld->SetStepPrecision(0.003);
166 fTrackHitsOld->SetMaxDistance(100);
173 // Initialise counters
179 // Initialise color attributes
180 SetMarkerColor(kYellow);
183 // Set TPC parameters
187 if (!strcmp(title,"Default")) {
188 fTPCParam = new AliTPCParamSR;
190 cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
196 //_____________________________________________________________________________
197 AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
199 // dummy copy constructor
212 delete fTrackHits; //MI 15.09.2000
213 delete fTrackHitsOld; //MI 10.12.2001
214 if (fNoiseTable) delete [] fNoiseTable;
218 //_____________________________________________________________________________
219 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
222 // Add a hit to the list
224 // TClonesArray &lhits = *fHits;
225 // new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
227 TClonesArray &lhits = *fHits;
228 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
231 AddHit2(track,vol,hits);
234 //_____________________________________________________________________________
235 void AliTPC::BuildGeometry()
239 // Build TPC ROOT TNode geometry for the event display
244 const int kColorTPC=19;
245 char name[5], title[25];
246 const Double_t kDegrad=TMath::Pi()/180;
247 const Double_t kRaddeg=180./TMath::Pi();
250 Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
251 Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
253 Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
254 Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
256 Int_t nLo = fTPCParam->GetNInnerSector()/2;
257 Int_t nHi = fTPCParam->GetNOuterSector()/2;
259 const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
260 const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
261 const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
262 const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
265 const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
266 const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
272 // Get ALICE top node
275 nTop=gAlice->GetGeometry()->GetNode("alice");
279 rl = fTPCParam->GetInnerRadiusLow();
280 ru = fTPCParam->GetInnerRadiusUp();
284 sprintf(name,"LS%2.2d",i);
286 sprintf(title,"TPC low sector %3d",i);
289 tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
290 kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
291 tubs->SetNumberOfDivisions(1);
293 nNode = new TNode(name,title,name,0,0,0,"");
294 nNode->SetLineColor(kColorTPC);
300 rl = fTPCParam->GetOuterRadiusLow();
301 ru = fTPCParam->GetOuterRadiusUp();
304 sprintf(name,"US%2.2d",i);
306 sprintf(title,"TPC upper sector %d",i);
308 tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
309 khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
310 tubs->SetNumberOfDivisions(1);
312 nNode = new TNode(name,title,name,0,0,0,"");
313 nNode->SetLineColor(kColorTPC);
319 //_____________________________________________________________________________
320 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) const
323 // Calculate distance from TPC to mouse on the display
329 void AliTPC::Clusters2Tracks() const
331 //-----------------------------------------------------------------
332 // This is a track finder.
333 //-----------------------------------------------------------------
334 Error("Clusters2Tracks",
335 "Dummy function ! Call AliTPCtracker::Clusters2Tracks(...) instead !");
339 //_____________________________________________________________________________
340 void AliTPC::Reconstruct() const
342 // reconstruct clusters
344 AliLoader* loader = GetLoader();
345 loader->LoadRecPoints("recreate");
346 loader->LoadDigits("read");
348 AliTPCclustererMI clusterer(fTPCParam);
349 AliRunLoader* runLoader = loader->GetRunLoader();
350 Int_t nEvents = runLoader->GetNumberOfEvents();
352 for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
353 runLoader->GetEvent(iEvent);
355 TTree* treeClusters = loader->TreeR();
357 loader->MakeTree("R");
358 treeClusters = loader->TreeR();
360 TTree* treeDigits = loader->TreeD();
362 Error("Reconstruct", "Can't get digits tree !");
366 // clusterer.Digits2Clusters(treeDigits, treeClusters);
367 clusterer.SetInput(treeDigits);
368 clusterer.SetOutput(treeClusters);
369 clusterer.Digits2Clusters();
371 loader->WriteRecPoints("OVERWRITE");
374 loader->UnloadRecPoints();
375 loader->UnloadDigits();
378 //_____________________________________________________________________________
379 AliTracker* AliTPC::CreateTracker() const
381 // create a TPC tracker
383 return new AliTPCtrackerMI(fTPCParam);
386 //_____________________________________________________________________________
387 void AliTPC::FillESD(AliESD* esd) const
391 Double_t parTPC[] = {47., 0.10, 10.};
392 AliTPCpidESD tpcPID(parTPC);
397 //_____________________________________________________________________________
398 void AliTPC::CreateMaterials()
400 //-----------------------------------------------
401 // Create Materials for for TPC simulations
402 //-----------------------------------------------
404 //-----------------------------------------------------------------
405 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
406 //-----------------------------------------------------------------
408 Int_t iSXFLD=gAlice->Field()->Integ();
409 Float_t sXMGMX=gAlice->Field()->Max();
411 Float_t amat[5]; // atomic numbers
412 Float_t zmat[5]; // z
413 Float_t wmat[5]; // proportions
419 //***************** Gases *************************
421 //-------------------------------------------------
423 //-------------------------------------------------
434 AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
444 AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
447 //--------------------------------------------------------------
449 //--------------------------------------------------------------
466 amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
468 AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
483 amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
485 AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
501 amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
503 AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
505 //----------------------------------------------------------------
506 // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
507 //----------------------------------------------------------------
513 Float_t rho,absl,x0,buf[1];
517 for(nc = 0;nc<fNoComp;nc++)
520 // retrive material constants
522 gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf);
527 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
529 am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
530 density += fMixtProp[nc]*rho; // density of the mixture
534 // mixture proportions by weight!
536 for(nc = 0;nc<fNoComp;nc++)
539 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
541 wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
542 apure[nnc] : amol[nnc])/am;
546 // Drift gases 1 - nonsensitive, 2 - sensitive
548 AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
549 AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
558 AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
561 //----------------------------------------------------------------------
563 //----------------------------------------------------------------------
585 AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
607 AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
625 AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
643 AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
645 // SiO2 - used later for the glass fiber
657 AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2)
666 AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
675 AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
684 AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
702 AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
721 AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
723 // Epoxy - C14 H20 O3
740 AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat);
748 AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.);
752 gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,x0,absl,buf,nbuf);
756 wmat[0]=0.644; // by weight!
759 density=0.5*(1.25+2.265);
761 AliMixture(47,"Cfiber",amat,zmat,density,2,wmat);
765 gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,x0,absl,buf,nbuf);
767 wmat[0]=0.725; // by weight!
772 AliMixture(39,"G10",amat,zmat,density,2,wmat);
777 //----------------------------------------------------------
778 // tracking media for gases
779 //----------------------------------------------------------
781 AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
782 AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
783 AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
784 AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
786 //-----------------------------------------------------------
787 // tracking media for solids
788 //-----------------------------------------------------------
790 AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
791 AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
792 AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
793 AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
794 AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
795 AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
796 AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
797 AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
798 AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
799 AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
800 AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
801 AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
805 void AliTPC::GenerNoise(Int_t tablesize)
808 //Generate table with noise
815 if (fNoiseTable) delete[] fNoiseTable;
816 fNoiseTable = new Float_t[tablesize];
817 fNoiseDepth = tablesize;
818 fCurrentNoise =0; //!index of the noise in the noise table
820 Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
821 for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
824 Float_t AliTPC::GetNoise()
826 // get noise from table
827 // if ((fCurrentNoise%10)==0)
828 // fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
829 if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
830 return fNoiseTable[fCurrentNoise++];
831 //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
835 Bool_t AliTPC::IsSectorActive(Int_t sec) const
838 // check if the sector is active
839 if (!fActiveSectors) return kTRUE;
840 else return fActiveSectors[sec];
843 void AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
845 // activate interesting sectors
846 SetTreeAddress();//just for security
847 if (fActiveSectors) delete [] fActiveSectors;
848 fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
849 for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
850 for (Int_t i=0;i<n;i++)
851 if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector()) fActiveSectors[sectors[i]]=kTRUE;
855 void AliTPC::SetActiveSectors(Int_t flag)
858 // activate sectors which were hitted by tracks
860 SetTreeAddress();//just for security
861 if (fHitType==0) return; // if Clones hit - not short volume ID information
862 if (fActiveSectors) delete [] fActiveSectors;
863 fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
865 for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
868 for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
872 Fatal("SetActiveSectors","Can not find TreeH in folder");
875 if (fHitType>1) branch = TreeH()->GetBranch("TPC2");
876 else branch = TreeH()->GetBranch("TPC");
877 Stat_t ntracks = TreeH()->GetEntries();
878 // loop over all hits
879 cout<<"\nAliTPC::SetActiveSectors(): Got "<<ntracks<<" tracks\n";
881 for(Int_t track=0;track<ntracks;track++)
885 if (fTrackHits && fHitType&4) {
886 TBranch * br1 = TreeH()->GetBranch("fVolumes");
887 TBranch * br2 = TreeH()->GetBranch("fNVolumes");
888 br1->GetEvent(track);
889 br2->GetEvent(track);
890 Int_t *volumes = fTrackHits->GetVolumes();
891 for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
892 fActiveSectors[volumes[j]]=kTRUE;
896 if (fTrackHitsOld && fHitType&2) {
897 TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
899 AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
900 for (UInt_t j=0;j<ar->GetSize();j++){
901 fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
910 void AliTPC::Digits2Clusters(Int_t /*eventnumber*/) const
912 //-----------------------------------------------------------------
913 // This is a simple cluster finder.
914 //-----------------------------------------------------------------
915 Error("Digits2Clusters",
916 "Dummy function ! Call AliTPCclusterer::Digits2Clusters(...) instead !");
919 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
920 extern Double_t SigmaZ2(Double_t, Double_t);
921 //_____________________________________________________________________________
922 void AliTPC::Hits2Clusters(Int_t /*eventn*/)
924 //--------------------------------------------------------
925 // TPC simple cluster generator from hits
926 // obtained from the TPC Fast Simulator
927 // The point errors are taken from the parametrization
928 //--------------------------------------------------------
930 //-----------------------------------------------------------------
931 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
932 //-----------------------------------------------------------------
933 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
934 // Jouri.Belikov@cern.ch
935 //----------------------------------------------------------------
937 /////////////////////////////////////////////////////////////////////////////
939 //---------------------------------------------------------------------
940 // ALICE TPC Cluster Parameters
941 //--------------------------------------------------------------------
945 // Cluster width in rphi
946 const Float_t kACrphi=0.18322;
947 const Float_t kBCrphi=0.59551e-3;
948 const Float_t kCCrphi=0.60952e-1;
949 // Cluster width in z
950 const Float_t kACz=0.19081;
951 const Float_t kBCz=0.55938e-3;
952 const Float_t kCCz=0.30428;
956 cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
960 //if(fDefaults == 0) SetDefaults();
962 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
964 TParticle *particle; // pointer to a given particle
965 AliTPChit *tpcHit; // pointer to a sigle TPC hit
969 Float_t pl,pt,tanth,rpad,ratio;
972 //---------------------------------------------------------------
973 // Get the access to the tracks
974 //---------------------------------------------------------------
979 Fatal("Hits2Clusters","Can not find TreeH in folder");
984 Stat_t ntracks = tH->GetEntries();
986 //Switch to the output file
988 if (fLoader->TreeR() == 0x0) fLoader->MakeTree("R");
990 cout<<"fTPCParam->GetTitle() = "<<fTPCParam->GetTitle()<<endl;
992 AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
994 //fTPCParam->Write(fTPCParam->GetTitle());
996 AliTPCClustersArray carray;
997 carray.Setup(fTPCParam);
998 carray.SetClusterType("AliTPCcluster");
999 carray.MakeTree(fLoader->TreeR());
1001 Int_t nclusters=0; //cluster counter
1003 //------------------------------------------------------------
1004 // Loop over all sectors (72 sectors for 20 deg
1005 // segmentation for both lower and upper sectors)
1006 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
1007 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
1009 // First cluster for sector 0 starts at "0"
1010 //------------------------------------------------------------
1012 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
1014 fTPCParam->AdjustCosSin(isec,cph,sph);
1016 //------------------------------------------------------------
1018 //------------------------------------------------------------
1020 for(Int_t track=0;track<ntracks;track++){
1022 tH->GetEvent(track);
1024 // Get number of the TPC hits
1026 tpcHit = (AliTPChit*)FirstHit(-1);
1032 if (tpcHit->fQ == 0.) {
1033 tpcHit = (AliTPChit*) NextHit();
1034 continue; //information about track (I.Belikov)
1036 sector=tpcHit->fSector; // sector number
1039 tpcHit = (AliTPChit*) NextHit();
1042 ipart=tpcHit->Track();
1043 particle=gAlice->GetMCApp()->Particle(ipart);
1046 if(pt < 1.e-9) pt=1.e-9;
1048 tanth = TMath::Abs(tanth);
1049 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
1050 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1052 // space-point resolutions
1054 sigmaRphi=SigmaY2(rpad,tanth,pt);
1055 sigmaZ =SigmaZ2(rpad,tanth );
1059 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
1060 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
1062 // temporary protection
1064 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
1065 if(sigmaZ < 0.) sigmaZ=0.4e-3;
1066 if(clRphi < 0.) clRphi=2.5e-3;
1067 if(clZ < 0.) clZ=2.5e-5;
1072 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1073 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1075 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
1076 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
1077 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
1078 Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
1079 fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
1080 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
1081 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
1082 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
1083 if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z();
1084 xyz[2]=sigmaRphi; // fSigmaY2
1085 xyz[3]=sigmaZ; // fSigmaZ2
1086 xyz[4]=tpcHit->fQ; // q
1088 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
1089 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
1091 Int_t tracks[3]={tpcHit->Track(), -1, -1};
1092 AliTPCcluster cluster(tracks,xyz);
1094 clrow->InsertCluster(&cluster); nclusters++;
1096 tpcHit = (AliTPChit*)NextHit();
1099 } // end of loop over hits
1101 } // end of loop over tracks
1103 Int_t nrows=fTPCParam->GetNRow(isec);
1104 for (Int_t irow=0; irow<nrows; irow++) {
1105 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
1106 if (!clrow) continue;
1107 carray.StoreRow(isec,irow);
1108 carray.ClearRow(isec,irow);
1111 } // end of loop over sectors
1113 cerr<<"Number of made clusters : "<<nclusters<<" \n";
1114 fLoader->WriteRecPoints("OVERWRITE");
1117 } // end of function
1119 //_________________________________________________________________
1120 void AliTPC::Hits2ExactClustersSector(Int_t isec)
1122 //--------------------------------------------------------
1123 //calculate exact cross point of track and given pad row
1124 //resulting values are expressed in "digit" coordinata
1125 //--------------------------------------------------------
1127 //-----------------------------------------------------------------
1128 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1129 //-----------------------------------------------------------------
1131 if (fClustersArray==0){
1135 TParticle *particle; // pointer to a given particle
1136 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1137 // Int_t sector,nhits;
1139 const Int_t kcmaxhits=30000;
1140 AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
1141 AliTPCFastVector & xxx = *xxxx;
1142 Int_t maxhits = kcmaxhits;
1143 //construct array for each padrow
1144 for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
1145 fClustersArray->CreateRow(isec,i);
1147 //---------------------------------------------------------------
1148 // Get the access to the tracks
1149 //---------------------------------------------------------------
1151 TTree *tH = TreeH();
1154 Fatal("Hits2Clusters","Can not find TreeH in folder");
1159 Stat_t ntracks = tH->GetEntries();
1160 Int_t npart = gAlice->GetMCApp()->GetNtrack();
1163 if (fHitType>1) branch = tH->GetBranch("TPC2");
1164 else branch = tH->GetBranch("TPC");
1166 //------------------------------------------------------------
1168 //------------------------------------------------------------
1170 for(Int_t track=0;track<ntracks;track++){
1171 Bool_t isInSector=kTRUE;
1173 isInSector = TrackInVolume(isec,track);
1174 if (!isInSector) continue;
1176 branch->GetEntry(track); // get next track
1178 // Get number of the TPC hits and a pointer
1182 Int_t currentIndex=0;
1183 Int_t lastrow=-1; //last writen row
1187 tpcHit = (AliTPChit*)FirstHit(-1);
1190 Int_t sector=tpcHit->fSector; // sector number
1192 tpcHit = (AliTPChit*) NextHit();
1196 ipart=tpcHit->Track();
1197 if (ipart<npart) particle=gAlice->GetMCApp()->Particle(ipart);
1201 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
1202 Int_t index[3]={1,isec,0};
1203 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
1204 if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;}
1205 if (lastrow<0) lastrow=currentrow;
1206 if (currentrow==lastrow){
1207 if ( currentIndex>=maxhits){
1209 xxx.ResizeTo(4*maxhits);
1211 xxx(currentIndex*4)=x[0];
1212 xxx(currentIndex*4+1)=x[1];
1213 xxx(currentIndex*4+2)=x[2];
1214 xxx(currentIndex*4+3)=tpcHit->fQ;
1218 if (currentIndex>2){
1230 for (Int_t index=0;index<currentIndex;index++){
1232 x=x2=x3=x4=xxx(index*4);
1240 sumy+=xxx(index*4+1);
1241 sumxy+=xxx(index*4+1)*x;
1242 sumx2y+=xxx(index*4+1)*x2;
1243 sumz+=xxx(index*4+2);
1244 sumxz+=xxx(index*4+2)*x;
1245 sumx2z+=xxx(index*4+2)*x2;
1246 sumq+=xxx(index*4+3);
1248 Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
1249 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1250 sumx2*(sumx*sumx3-sumx2*sumx2);
1252 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1253 sumx2*(sumxy*sumx3-sumx2y*sumx2);
1254 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1255 sumx2*(sumxz*sumx3-sumx2z*sumx2);
1257 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1258 sumx2*(sumx*sumx2y-sumx2*sumxy);
1259 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1260 sumx2*(sumx*sumx2z-sumx2*sumxz);
1262 if (TMath::Abs(det)<0.00001){
1263 tpcHit = (AliTPChit*)NextHit();
1267 Float_t y=detay/det+centralPad;
1268 Float_t z=detaz/det;
1269 Float_t by=detby/det; //y angle
1270 Float_t bz=detbz/det; //z angle
1271 sumy/=Float_t(currentIndex);
1272 sumz/=Float_t(currentIndex);
1274 AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
1276 AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
1277 // AliComplexCluster cl;
1283 cl->fTracks[0]=ipart;
1287 } //end of calculating cluster for given row
1290 tpcHit = (AliTPChit*)NextHit();
1291 } // end of loop over hits
1292 } // end of loop over tracks
1293 //write padrows to tree
1294 for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1295 fClustersArray->StoreRow(isec,ii);
1296 fClustersArray->ClearRow(isec,ii);
1304 //______________________________________________________________________
1305 AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
1307 return new AliTPCDigitizer(manager);
1310 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
1312 //create digits from summable digits
1313 GenerNoise(500000); //create teble with noise
1315 //conect tree with sSDigits
1316 TTree *t = fLoader->TreeS();
1320 fLoader->LoadSDigits("READ");
1321 t = fLoader->TreeS();
1324 Error("SDigits2Digits2","Can not get input TreeS");
1329 if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
1331 AliSimDigits digarr, *dummy=&digarr;
1332 TBranch* sdb = t->GetBranch("Segment");
1335 Error("SDigits2Digits2","Can not find branch with segments in TreeS.");
1339 sdb->SetAddress(&dummy);
1341 Stat_t nentries = t->GetEntries();
1343 // set zero suppression
1345 fTPCParam->SetZeroSup(2);
1347 // get zero suppression
1349 Int_t zerosup = fTPCParam->GetZeroSup();
1351 //make tree with digits
1353 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1354 arr->SetClass("AliSimDigits");
1355 arr->Setup(fTPCParam);
1356 arr->MakeTree(fLoader->TreeD());
1358 AliTPCParam * par = fTPCParam;
1360 //Loop over segments of the TPC
1362 for (Int_t n=0; n<nentries; n++) {
1365 if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1366 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
1369 if (!IsSectorActive(sec))
1371 // cout<<n<<" NOT Active \n";
1376 // cout<<n<<" Active \n";
1378 AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1379 Int_t nrows = digrow->GetNRows();
1380 Int_t ncols = digrow->GetNCols();
1382 digrow->ExpandBuffer();
1383 digarr.ExpandBuffer();
1384 digrow->ExpandTrackBuffer();
1385 digarr.ExpandTrackBuffer();
1388 Short_t * pamp0 = digarr.GetDigits();
1389 Int_t * ptracks0 = digarr.GetTracks();
1390 Short_t * pamp1 = digrow->GetDigits();
1391 Int_t * ptracks1 = digrow->GetTracks();
1392 Int_t nelems =nrows*ncols;
1393 Int_t saturation = fTPCParam->GetADCSat();
1394 //use internal structure of the AliDigits - for speed reason
1395 //if you cahnge implementation
1396 //of the Alidigits - it must be rewriten -
1397 for (Int_t i= 0; i<nelems; i++){
1398 // Float_t q = *pamp0;
1399 //q/=16.; //conversion faktor
1400 //Float_t noise= GetNoise();
1402 //q= TMath::Nint(q);
1403 Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1405 if (q>saturation) q=saturation;
1407 //if (ptracks0[0]==0)
1410 ptracks1[0]=ptracks0[0];
1411 ptracks1[nelems]=ptracks0[nelems];
1412 ptracks1[2*nelems]=ptracks0[2*nelems];
1420 arr->StoreRow(sec,row);
1421 arr->ClearRow(sec,row);
1422 // cerr<<sec<<"\t"<<row<<"\n";
1427 fLoader->WriteDigits("OVERWRITE");
1431 //__________________________________________________________________
1432 void AliTPC::SetDefaults(){
1434 // setting the defaults
1437 cerr<<"Setting default parameters...\n";
1439 // Set response functions
1442 AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
1444 AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1446 printf("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...\n");
1448 param = new AliTPCParamSR();
1451 param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1454 printf("No TPC parameters found\n");
1459 AliTPCPRF2D * prfinner = new AliTPCPRF2D;
1460 AliTPCPRF2D * prfouter1 = new AliTPCPRF2D;
1461 AliTPCPRF2D * prfouter2 = new AliTPCPRF2D;
1462 AliTPCRF1D * rf = new AliTPCRF1D(kTRUE);
1463 rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1464 rf->SetOffset(3*param->GetZSigma());
1467 TDirectory *savedir=gDirectory;
1468 TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1470 cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
1475 prfinner->Read("prf_07504_Gati_056068_d02");
1476 //PH Set different names
1477 s=prfinner->GetGRF()->GetName();
1479 prfinner->GetGRF()->SetName(s.Data());
1481 prfouter1->Read("prf_10006_Gati_047051_d03");
1482 s=prfouter1->GetGRF()->GetName();
1484 prfouter1->GetGRF()->SetName(s.Data());
1486 prfouter2->Read("prf_15006_Gati_047051_d03");
1487 s=prfouter2->GetGRF()->GetName();
1489 prfouter2->GetGRF()->SetName(s.Data());
1494 param->SetInnerPRF(prfinner);
1495 param->SetOuter1PRF(prfouter1);
1496 param->SetOuter2PRF(prfouter2);
1497 param->SetTimeRF(rf);
1507 //__________________________________________________________________
1508 void AliTPC::Hits2Digits()
1511 // creates digits from hits
1514 fLoader->LoadHits("read");
1515 fLoader->LoadDigits("recreate");
1516 AliRunLoader* runLoader = fLoader->GetRunLoader();
1518 for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1519 runLoader->GetEvent(iEvent);
1521 Hits2Digits(iEvent);
1524 fLoader->UnloadHits();
1525 fLoader->UnloadDigits();
1527 //__________________________________________________________________
1528 void AliTPC::Hits2Digits(Int_t eventnumber)
1530 //----------------------------------------------------
1531 // Loop over all sectors for a single event
1532 //----------------------------------------------------
1533 AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
1534 rl->GetEvent(eventnumber);
1535 if (fLoader->TreeH() == 0x0)
1537 if(fLoader->LoadHits())
1539 Error("Hits2Digits","Can not load hits.");
1544 if (fLoader->TreeD() == 0x0 )
1546 fLoader->MakeTree("D");
1547 if (fLoader->TreeD() == 0x0 )
1549 Error("Hits2Digits","Can not get TreeD");
1554 if(fDefaults == 0) SetDefaults(); // check if the parameters are set
1555 GenerNoise(500000); //create teble with noise
1557 //setup TPCDigitsArray
1559 if(GetDigitsArray()) delete GetDigitsArray();
1561 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1562 arr->SetClass("AliSimDigits");
1563 arr->Setup(fTPCParam);
1565 arr->MakeTree(fLoader->TreeD());
1566 SetDigitsArray(arr);
1568 fDigitsSwitch=0; // standard digits
1570 cerr<<"Digitizing TPC -- normal digits...\n";
1572 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1573 if (IsSectorActive(isec))
1575 if (fDebug) Info("Hits2Digits","Sector %d is active.",isec);
1576 Hits2DigitsSector(isec);
1580 if (fDebug) Info("Hits2Digits","Sector %d is NOT active.",isec);
1583 fLoader->WriteDigits("OVERWRITE");
1585 //this line prevents the crash in the similar one
1586 //on the beginning of this method
1587 //destructor attempts to reset the tree, which is deleted by the loader
1588 //need to be redesign
1589 if(GetDigitsArray()) delete GetDigitsArray();
1590 SetDigitsArray(0x0);
1594 //__________________________________________________________________
1595 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1598 //-----------------------------------------------------------
1599 // summable digits - 16 bit "ADC", no noise, no saturation
1600 //-----------------------------------------------------------
1602 //----------------------------------------------------
1603 // Loop over all sectors for a single event
1604 //----------------------------------------------------
1605 // AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
1607 AliRunLoader* rl = fLoader->GetRunLoader();
1609 rl->GetEvent(eventnumber);
1610 if (fLoader->TreeH() == 0x0)
1612 if(fLoader->LoadHits())
1614 Error("Hits2Digits","Can not load hits.");
1621 if (fLoader->TreeS() == 0x0 )
1623 fLoader->MakeTree("S");
1626 if(fDefaults == 0) SetDefaults();
1628 GenerNoise(500000); //create table with noise
1629 //setup TPCDigitsArray
1631 if(GetDigitsArray()) delete GetDigitsArray();
1634 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1635 arr->SetClass("AliSimDigits");
1636 arr->Setup(fTPCParam);
1637 arr->MakeTree(fLoader->TreeS());
1639 SetDigitsArray(arr);
1641 cerr<<"Digitizing TPC -- summable digits...\n";
1643 fDigitsSwitch=1; // summable digits
1645 // set zero suppression to "0"
1647 fTPCParam->SetZeroSup(0);
1649 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1650 if (IsSectorActive(isec))
1652 // cout<<"Sector "<<isec<<" is active\n";
1653 Hits2DigitsSector(isec);
1656 fLoader->WriteSDigits("OVERWRITE");
1658 //this line prevents the crash in the similar one
1659 //on the beginning of this method
1660 //destructor attempts to reset the tree, which is deleted by the loader
1661 //need to be redesign
1662 if(GetDigitsArray()) delete GetDigitsArray();
1663 SetDigitsArray(0x0);
1665 //__________________________________________________________________
1667 void AliTPC::Hits2SDigits()
1670 //-----------------------------------------------------------
1671 // summable digits - 16 bit "ADC", no noise, no saturation
1672 //-----------------------------------------------------------
1674 fLoader->LoadHits("read");
1675 fLoader->LoadSDigits("recreate");
1676 AliRunLoader* runLoader = fLoader->GetRunLoader();
1678 for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1679 runLoader->GetEvent(iEvent);
1682 Hits2SDigits2(iEvent);
1685 fLoader->UnloadHits();
1686 fLoader->UnloadSDigits();
1688 //_____________________________________________________________________________
1690 void AliTPC::Hits2DigitsSector(Int_t isec)
1692 //-------------------------------------------------------------------
1693 // TPC conversion from hits to digits.
1694 //-------------------------------------------------------------------
1696 //-----------------------------------------------------------------
1697 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1698 //-----------------------------------------------------------------
1700 //-------------------------------------------------------
1701 // Get the access to the track hits
1702 //-------------------------------------------------------
1704 // check if the parameters are set - important if one calls this method
1705 // directly, not from the Hits2Digits
1707 if(fDefaults == 0) SetDefaults();
1709 TTree *tH = TreeH(); // pointer to the hits tree
1712 Fatal("Hits2DigitsSector","Can not find TreeH in folder");
1716 Stat_t ntracks = tH->GetEntries();
1720 //-------------------------------------------
1721 // Only if there are any tracks...
1722 //-------------------------------------------
1726 //printf("*** Processing sector number %d ***\n",isec);
1728 Int_t nrows =fTPCParam->GetNRow(isec);
1730 row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1732 MakeSector(isec,nrows,tH,ntracks,row);
1734 //--------------------------------------------------------
1735 // Digitize this sector, row by row
1736 // row[i] is the pointer to the TObjArray of AliTPCFastVectors,
1737 // each one containing electrons accepted on this
1738 // row, assigned into tracks
1739 //--------------------------------------------------------
1743 if (fDigitsArray->GetTree()==0)
1745 Fatal("Hits2DigitsSector","Tree not set in fDigitsArray");
1748 for (i=0;i<nrows;i++){
1750 AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1752 DigitizeRow(i,isec,row);
1754 fDigitsArray->StoreRow(isec,i);
1756 Int_t ndig = dig->GetDigitSize();
1759 printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1761 fDigitsArray->ClearRow(isec,i);
1764 } // end of the sector digitization
1766 for(i=0;i<nrows+2;i++){
1771 delete [] row; // delete the array of pointers to TObjArray-s
1775 } // end of Hits2DigitsSector
1778 //_____________________________________________________________________________
1779 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1781 //-----------------------------------------------------------
1782 // Single row digitization, coupling from the neighbouring
1783 // rows taken into account
1784 //-----------------------------------------------------------
1786 //-----------------------------------------------------------------
1787 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1788 // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1789 //-----------------------------------------------------------------
1792 Float_t zerosup = fTPCParam->GetZeroSup();
1793 // Int_t nrows =fTPCParam->GetNRow(isec);
1794 fCurrentIndex[1]= isec;
1797 Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1798 Int_t nofTbins = fTPCParam->GetMaxTBin();
1799 Int_t indexRange[4];
1801 // Integrated signal for this row
1802 // and a single track signal
1805 AliTPCFastMatrix *m1 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // integrated
1806 AliTPCFastMatrix *m2 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // single
1808 AliTPCFastMatrix &total = *m1;
1810 // Array of pointers to the label-signal list
1812 Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1813 Float_t **pList = new Float_t* [nofDigits];
1817 for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1821 //Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1822 //Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1825 for (Int_t row= row1;row<=row2;row++){
1826 Int_t nTracks= rows[row]->GetEntries();
1827 for (i1=0;i1<nTracks;i1++){
1828 fCurrentIndex[2]= row;
1829 fCurrentIndex[3]=irow+1;
1831 m2->Zero(); // clear single track signal matrix
1832 Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1833 GetList(trackLabel,nofPads,m2,indexRange,pList);
1835 else GetSignal(rows[row],i1,0,m1,indexRange);
1841 AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1843 Float_t fzerosup = zerosup+0.5;
1844 for(Int_t it=0;it<nofTbins;it++){
1845 Float_t *pq = &(total.UncheckedAt(0,it));
1846 for(Int_t ip=0;ip<nofPads;ip++){
1850 if(fDigitsSwitch == 0){
1852 if(q <=fzerosup) continue; // do not fill zeros
1854 if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation
1859 if(q <= 0.) continue; // do not fill zeros
1860 if(q>2000.) q=2000.;
1866 // "real" signal or electronic noise (list = -1)?
1869 for(Int_t j1=0;j1<3;j1++){
1870 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1875 <A NAME="AliDigits"></A>
1876 using of AliDigits object
1879 dig->SetDigitFast((Short_t)q,it,ip);
1880 if (fDigitsArray->IsSimulated())
1882 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1883 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1884 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1888 } // end of loop over time buckets
1889 } // end of lop over pads
1892 // This row has been digitized, delete nonused stuff
1895 for(lp=0;lp<nofDigits;lp++){
1896 if(pList[lp]) delete [] pList[lp];
1905 } // end of DigitizeRow
1907 //_____________________________________________________________________________
1909 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1910 AliTPCFastMatrix *m1, AliTPCFastMatrix *m2,Int_t *indexRange)
1913 //---------------------------------------------------------------
1914 // Calculates 2-D signal (pad,time) for a single track,
1915 // returns a pointer to the signal matrix and the track label
1916 // No digitization is performed at this level!!!
1917 //---------------------------------------------------------------
1919 //-----------------------------------------------------------------
1920 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1921 // Modified: Marian Ivanov
1922 //-----------------------------------------------------------------
1924 AliTPCFastVector *tv;
1926 tv = (AliTPCFastVector*)p1->At(ntr); // pointer to a track
1927 AliTPCFastVector &v = *tv;
1929 Float_t label = v(0);
1930 Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2;
1932 Int_t nElectrons = (tv->GetNrows()-1)/4;
1933 indexRange[0]=9999; // min pad
1934 indexRange[1]=-1; // max pad
1935 indexRange[2]=9999; //min time
1936 indexRange[3]=-1; // max time
1938 AliTPCFastMatrix &signal = *m1;
1939 AliTPCFastMatrix &total = *m2;
1941 // Loop over all electrons
1943 for(Int_t nel=0; nel<nElectrons; nel++){
1945 Float_t aval = v(idx+4);
1946 Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1947 Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1948 Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
1950 Int_t *index = fTPCParam->GetResBin(0);
1951 Float_t *weight = & (fTPCParam->GetResWeight(0));
1953 if (n>0) for (Int_t i =0; i<n; i++){
1954 Int_t pad=index[1]+centralPad; //in digit coordinates central pad has coordinate 0
1957 Int_t time=index[2];
1958 Float_t qweight = *(weight)*eltoadcfac;
1960 if (m1!=0) signal.UncheckedAt(pad,time)+=qweight;
1961 total.UncheckedAt(pad,time)+=qweight;
1962 if (indexRange[0]>pad) indexRange[0]=pad;
1963 if (indexRange[1]<pad) indexRange[1]=pad;
1964 if (indexRange[2]>time) indexRange[2]=time;
1965 if (indexRange[3]<time) indexRange[3]=time;
1972 } // end of loop over electrons
1974 return label; // returns track label when finished
1977 //_____________________________________________________________________________
1978 void AliTPC::GetList(Float_t label,Int_t np,AliTPCFastMatrix *m,
1979 Int_t *indexRange, Float_t **pList)
1981 //----------------------------------------------------------------------
1982 // Updates the list of tracks contributing to digits for a given row
1983 //----------------------------------------------------------------------
1985 //-----------------------------------------------------------------
1986 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1987 //-----------------------------------------------------------------
1989 AliTPCFastMatrix &signal = *m;
1991 // lop over nonzero digits
1993 for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1994 for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1997 // accept only the contribution larger than 500 electrons (1/2 s_noise)
1999 if(signal(ip,it)<0.5) continue;
2002 Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
2004 if(!pList[globalIndex]){
2007 // Create new list (6 elements - 3 signals and 3 labels),
2010 pList[globalIndex] = new Float_t [6];
2014 *pList[globalIndex] = -1.;
2015 *(pList[globalIndex]+1) = -1.;
2016 *(pList[globalIndex]+2) = -1.;
2017 *(pList[globalIndex]+3) = -1.;
2018 *(pList[globalIndex]+4) = -1.;
2019 *(pList[globalIndex]+5) = -1.;
2022 *pList[globalIndex] = label;
2023 *(pList[globalIndex]+3) = signal(ip,it);
2027 // check the signal magnitude
2029 Float_t highest = *(pList[globalIndex]+3);
2030 Float_t middle = *(pList[globalIndex]+4);
2031 Float_t lowest = *(pList[globalIndex]+5);
2034 // compare the new signal with already existing list
2037 if(signal(ip,it)<lowest) continue; // neglect this track
2041 if (signal(ip,it)>highest){
2042 *(pList[globalIndex]+5) = middle;
2043 *(pList[globalIndex]+4) = highest;
2044 *(pList[globalIndex]+3) = signal(ip,it);
2046 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
2047 *(pList[globalIndex]+1) = *pList[globalIndex];
2048 *pList[globalIndex] = label;
2050 else if (signal(ip,it)>middle){
2051 *(pList[globalIndex]+5) = middle;
2052 *(pList[globalIndex]+4) = signal(ip,it);
2054 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
2055 *(pList[globalIndex]+1) = label;
2058 *(pList[globalIndex]+5) = signal(ip,it);
2059 *(pList[globalIndex]+2) = label;
2063 } // end of loop over pads
2064 } // end of loop over time bins
2069 //___________________________________________________________________
2070 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
2071 Stat_t ntracks,TObjArray **row)
2074 //-----------------------------------------------------------------
2075 // Prepares the sector digitization, creates the vectors of
2076 // tracks for each row of this sector. The track vector
2077 // contains the track label and the position of electrons.
2078 //-----------------------------------------------------------------
2080 //-----------------------------------------------------------------
2081 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2082 //-----------------------------------------------------------------
2084 Float_t gasgain = fTPCParam->GetGasGain();
2088 AliTPChit *tpcHit; // pointer to a sigle TPC hit
2091 if (fHitType>1) branch = TH->GetBranch("TPC2");
2092 else branch = TH->GetBranch("TPC");
2095 //----------------------------------------------
2096 // Create TObjArray-s, one for each row,
2097 // each TObjArray will store the AliTPCFastVectors
2098 // of electrons, one AliTPCFastVectors per each track.
2099 //----------------------------------------------
2101 Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
2102 AliTPCFastVector **tracks = new AliTPCFastVector* [nrows+2]; //pointers to the track vectors
2104 for(i=0; i<nrows+2; i++){
2105 row[i] = new TObjArray;
2112 //--------------------------------------------------------------------
2113 // Loop over tracks, the "track" contains the full history
2114 //--------------------------------------------------------------------
2116 Int_t previousTrack,currentTrack;
2117 previousTrack = -1; // nothing to store so far!
2119 for(Int_t track=0;track<ntracks;track++){
2120 Bool_t isInSector=kTRUE;
2122 isInSector = TrackInVolume(isec,track);
2123 if (!isInSector) continue;
2125 branch->GetEntry(track); // get next track
2129 tpcHit = (AliTPChit*)FirstHit(-1);
2131 //--------------------------------------------------------------
2133 //--------------------------------------------------------------
2138 Int_t sector=tpcHit->fSector; // sector number
2140 tpcHit = (AliTPChit*) NextHit();
2144 currentTrack = tpcHit->Track(); // track number
2147 if(currentTrack != previousTrack){
2149 // store already filled fTrack
2151 for(i=0;i<nrows+2;i++){
2152 if(previousTrack != -1){
2153 if(nofElectrons[i]>0){
2154 AliTPCFastVector &v = *tracks[i];
2155 v(0) = previousTrack;
2156 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2157 row[i]->Add(tracks[i]);
2160 delete tracks[i]; // delete empty AliTPCFastVector
2166 tracks[i] = new AliTPCFastVector(481); // AliTPCFastVectors for the next fTrack
2168 } // end of loop over rows
2170 previousTrack=currentTrack; // update track label
2173 Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2175 //---------------------------------------------------
2176 // Calculate the electron attachment probability
2177 //---------------------------------------------------
2180 Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
2181 /fTPCParam->GetDriftV();
2183 Float_t attProb = fTPCParam->GetAttCoef()*
2184 fTPCParam->GetOxyCont()*time; // fraction!
2186 //-----------------------------------------------
2187 // Loop over electrons
2188 //-----------------------------------------------
2191 for(Int_t nel=0;nel<qI;nel++){
2192 // skip if electron lost due to the attachment
2193 if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2198 // protection for the nonphysical avalanche size (10**6 maximum)
2200 Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2201 xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
2204 TransportElectron(xyz,index);
2206 fTPCParam->GetPadRow(xyz,index);
2207 // row 0 - cross talk from the innermost row
2208 // row fNRow+1 cross talk from the outermost row
2209 rowNumber = index[2]+1;
2210 //transform position to local digit coordinates
2211 //relative to nearest pad row
2212 if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2214 if (isec <fTPCParam->GetNInnerSector()) {
2215 x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2216 y1 = fTPCParam->GetYInner(rowNumber);
2219 x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2220 y1 = fTPCParam->GetYOuter(rowNumber);
2222 // gain inefficiency at the wires edges - linear
2225 if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
2227 nofElectrons[rowNumber]++;
2228 //----------------------------------
2229 // Expand vector if necessary
2230 //----------------------------------
2231 if(nofElectrons[rowNumber]>120){
2232 Int_t range = tracks[rowNumber]->GetNrows();
2233 if((nofElectrons[rowNumber])>(range-1)/4){
2235 tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
2239 AliTPCFastVector &v = *tracks[rowNumber];
2240 Int_t idx = 4*nofElectrons[rowNumber]-3;
2241 Real_t * position = &(((AliTPCFastVector&)v).UncheckedAt(idx)); //make code faster
2242 memcpy(position,xyz,4*sizeof(Float_t));
2244 } // end of loop over electrons
2246 tpcHit = (AliTPChit*)NextHit();
2248 } // end of loop over hits
2249 } // end of loop over tracks
2252 // store remaining track (the last one) if not empty
2255 for(i=0;i<nrows+2;i++){
2256 if(nofElectrons[i]>0){
2257 AliTPCFastVector &v = *tracks[i];
2258 v(0) = previousTrack;
2259 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2260 row[i]->Add(tracks[i]);
2269 delete [] nofElectrons;
2272 } // end of MakeSector
2275 //_____________________________________________________________________________
2279 // Initialise TPC detector after definition of geometry
2284 printf("\n%s: ",ClassName());
2285 for(i=0;i<35;i++) printf("*");
2286 printf(" TPC_INIT ");
2287 for(i=0;i<35;i++) printf("*");
2288 printf("\n%s: ",ClassName());
2290 for(i=0;i<80;i++) printf("*");
2295 //_____________________________________________________________________________
2296 void AliTPC::MakeBranch(Option_t* option)
2299 // Create Tree branches for the TPC.
2301 if(GetDebug()) Info("MakeBranch","");
2302 Int_t buffersize = 4000;
2303 char branchname[10];
2304 sprintf(branchname,"%s",GetName());
2306 const char *h = strstr(option,"H");
2308 if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2310 AliDetector::MakeBranch(option);
2312 const char *d = strstr(option,"D");
2314 if (fDigits && fLoader->TreeD() && d)
2316 MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
2319 if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
2322 //_____________________________________________________________________________
2323 void AliTPC::ResetDigits()
2326 // Reset number of digits and the digits array for this detector
2329 if (fDigits) fDigits->Clear();
2332 //_____________________________________________________________________________
2333 void AliTPC::SetSecAL(Int_t sec)
2335 //---------------------------------------------------
2336 // Activate/deactivate selection for lower sectors
2337 //---------------------------------------------------
2339 //-----------------------------------------------------------------
2340 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2341 //-----------------------------------------------------------------
2345 //_____________________________________________________________________________
2346 void AliTPC::SetSecAU(Int_t sec)
2348 //----------------------------------------------------
2349 // Activate/deactivate selection for upper sectors
2350 //---------------------------------------------------
2352 //-----------------------------------------------------------------
2353 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2354 //-----------------------------------------------------------------
2358 //_____________________________________________________________________________
2359 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2361 //----------------------------------------
2362 // Select active lower sectors
2363 //----------------------------------------
2365 //-----------------------------------------------------------------
2366 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2367 //-----------------------------------------------------------------
2377 //_____________________________________________________________________________
2378 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2379 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2380 Int_t s11 , Int_t s12)
2382 //--------------------------------
2383 // Select active upper sectors
2384 //--------------------------------
2386 //-----------------------------------------------------------------
2387 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2388 //-----------------------------------------------------------------
2404 //_____________________________________________________________________________
2405 void AliTPC::SetSens(Int_t sens)
2408 //-------------------------------------------------------------
2409 // Activates/deactivates the sensitive strips at the center of
2410 // the pad row -- this is for the space-point resolution calculations
2411 //-------------------------------------------------------------
2413 //-----------------------------------------------------------------
2414 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2415 //-----------------------------------------------------------------
2421 void AliTPC::SetSide(Float_t side=0.)
2423 // choice of the TPC side
2428 //____________________________________________________________________________
2429 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2430 Float_t p2,Float_t p3)
2433 // gax mixture definition
2447 //_____________________________________________________________________________
2449 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2452 // electron transport taking into account:
2454 // 2.ExB at the wires
2455 // 3. nonisochronity
2457 // xyz and index must be already transformed to system 1
2460 fTPCParam->Transform1to2(xyz,index);
2463 Float_t driftl=xyz[2];
2464 if(driftl<0.01) driftl=0.01;
2465 driftl=TMath::Sqrt(driftl);
2466 Float_t sigT = driftl*(fTPCParam->GetDiffT());
2467 Float_t sigL = driftl*(fTPCParam->GetDiffL());
2468 xyz[0]=gRandom->Gaus(xyz[0],sigT);
2469 xyz[1]=gRandom->Gaus(xyz[1],sigT);
2470 xyz[2]=gRandom->Gaus(xyz[2],sigL);
2474 if (fTPCParam->GetMWPCReadout()==kTRUE){
2475 Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2476 xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2478 //add nonisochronity (not implemented yet)
2481 ClassImp(AliTPCdigit)
2483 //_____________________________________________________________________________
2484 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2488 // Creates a TPC digit object
2490 fSector = digits[0];
2491 fPadRow = digits[1];
2494 fSignal = digits[4];
2500 //_____________________________________________________________________________
2501 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2505 // Creates a TPC hit object
2515 //________________________________________________________________________
2516 // Additional code because of the AliTPCTrackHitsV2
2518 void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
2521 // Create a new branch in the current Root Tree
2522 // The branch of fHits is automatically split
2523 // MI change 14.09.2000
2524 if(GetDebug()) Info("MakeBranch2","");
2525 if (fHitType<2) return;
2526 char branchname[10];
2527 sprintf(branchname,"%s2",GetName());
2529 // Get the pointer to the header
2530 const char *cH = strstr(option,"H");
2532 if (fTrackHits && TreeH() && cH && fHitType&4)
2534 if(GetDebug()) Info("MakeBranch2","Making branch for Type 4 Hits");
2535 TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2538 if (fTrackHitsOld && TreeH() && cH && fHitType&2)
2540 if(GetDebug()) Info("MakeBranch2","Making branch for Type 2 Hits");
2541 AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2542 TreeH(),fBufferSize,99);
2543 TreeH()->GetListOfBranches()->Add(branch);
2547 void AliTPC::SetTreeAddress()
2549 //Sets tree address for hits
2552 if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2553 AliDetector::SetTreeAddress();
2555 if (fHitType>1) SetTreeAddress2();
2558 void AliTPC::SetTreeAddress2()
2561 // Set branch address for the TrackHits Tree
2563 if(GetDebug()) Info("SetTreeAddress2","");
2566 char branchname[20];
2567 sprintf(branchname,"%s2",GetName());
2569 // Branch address for hit tree
2570 TTree *treeH = TreeH();
2571 if ((treeH)&&(fHitType&4)) {
2572 branch = treeH->GetBranch(branchname);
2575 branch->SetAddress(&fTrackHits);
2576 if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Setting");
2579 if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Failed (can not find branch)");
2582 if ((treeH)&&(fHitType&2)) {
2583 branch = treeH->GetBranch(branchname);
2586 branch->SetAddress(&fTrackHitsOld);
2587 if (GetDebug()) Info("SetTreeAddress2","fHitType&2 Setting");
2589 else if (GetDebug())
2590 Info("SetTreeAddress2","fHitType&2 Failed (can not find branch)");
2592 //set address to TREETR
2594 TTree *treeTR = TreeTR();
2595 if (treeTR && fTrackReferences) {
2596 branch = treeTR->GetBranch(GetName());
2597 if (branch) branch->SetAddress(&fTrackReferences);
2602 void AliTPC::FinishPrimary()
2604 if (fTrackHits &&fHitType&4) fTrackHits->FlushHitStack();
2605 if (fTrackHitsOld && fHitType&2) fTrackHitsOld->FlushHitStack();
2609 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2612 // add hit to the list
2615 int primary = gAlice->GetMCApp()->GetPrimary(track);
2616 gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2620 gAlice->GetMCApp()->FlagTrack(track);
2622 //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
2623 //if (hit->fTrack!=rtrack)
2624 // cout<<"bad track number\n";
2625 if (fTrackHits && fHitType&4)
2626 fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2627 hits[1],hits[2],(Int_t)hits[3]);
2628 if (fTrackHitsOld &&fHitType&2 )
2629 fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2630 hits[1],hits[2],(Int_t)hits[3]);
2634 void AliTPC::ResetHits()
2636 if (fHitType&1) AliDetector::ResetHits();
2637 if (fHitType>1) ResetHits2();
2640 void AliTPC::ResetHits2()
2644 if (fTrackHits && fHitType&4) fTrackHits->Clear();
2645 if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2649 AliHit* AliTPC::FirstHit(Int_t track)
2651 if (fHitType>1) return FirstHit2(track);
2652 return AliDetector::FirstHit(track);
2654 AliHit* AliTPC::NextHit()
2659 if (fHitType>1) return NextHit2();
2661 return AliDetector::NextHit();
2664 AliHit* AliTPC::FirstHit2(Int_t track)
2667 // Initialise the hit iterator
2668 // Return the address of the first hit for track
2669 // If track>=0 the track is read from disk
2670 // while if track<0 the first hit of the current
2671 // track is returned
2674 gAlice->ResetHits();
2675 TreeH()->GetEvent(track);
2678 if (fTrackHits && fHitType&4) {
2679 fTrackHits->First();
2680 return fTrackHits->GetHit();
2682 if (fTrackHitsOld && fHitType&2) {
2683 fTrackHitsOld->First();
2684 return fTrackHitsOld->GetHit();
2690 AliHit* AliTPC::NextHit2()
2693 //Return the next hit for the current track
2696 if (fTrackHitsOld && fHitType&2) {
2697 fTrackHitsOld->Next();
2698 return fTrackHitsOld->GetHit();
2702 return fTrackHits->GetHit();
2708 void AliTPC::LoadPoints(Int_t)
2712 /* if(fHitType==1) return AliDetector::LoadPoints(a);
2715 if(fHitType==1) AliDetector::LoadPoints(a);
2716 else LoadPoints2(a);
2723 void AliTPC::RemapTrackHitIDs(Int_t *map)
2728 if (!fTrackHits) return;
2730 if (fTrackHitsOld && fHitType&2){
2731 AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2732 for (UInt_t i=0;i<arr->GetSize();i++){
2733 AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2734 info->fTrackID = map[info->fTrackID];
2737 if (fTrackHitsOld && fHitType&4){
2738 TClonesArray * arr = fTrackHits->GetArray();;
2739 for (Int_t i=0;i<arr->GetEntriesFast();i++){
2740 AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2741 info->fTrackID = map[info->fTrackID];
2746 Bool_t AliTPC::TrackInVolume(Int_t id,Int_t track)
2748 //return bool information - is track in given volume
2749 //load only part of the track information
2750 //return true if current track is in volume
2753 if (fTrackHitsOld && fHitType&2) {
2754 TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
2755 br->GetEvent(track);
2756 AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2757 for (UInt_t j=0;j<ar->GetSize();j++){
2758 if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2762 if (fTrackHits && fHitType&4) {
2763 TBranch * br1 = TreeH()->GetBranch("fVolumes");
2764 TBranch * br2 = TreeH()->GetBranch("fNVolumes");
2765 br2->GetEvent(track);
2766 br1->GetEvent(track);
2767 Int_t *volumes = fTrackHits->GetVolumes();
2768 Int_t nvolumes = fTrackHits->GetNVolumes();
2769 if (!volumes && nvolumes>0) {
2770 printf("Problematic track\t%d\t%d",track,nvolumes);
2773 for (Int_t j=0;j<nvolumes; j++)
2774 if (volumes[j]==id) return kTRUE;
2778 TBranch * br = TreeH()->GetBranch("fSector");
2779 br->GetEvent(track);
2780 for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2781 if ( ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2788 //_____________________________________________________________________________
2789 void AliTPC::LoadPoints2(Int_t)
2792 // Store x, y, z of all hits in memory
2794 if (fTrackHits == 0 && fTrackHitsOld==0) return;
2797 if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
2798 if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
2800 if (nhits == 0) return;
2801 Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2802 if (fPoints == 0) fPoints = new TObjArray(tracks);
2805 Int_t *ntrk=new Int_t[tracks];
2806 Int_t *limi=new Int_t[tracks];
2807 Float_t **coor=new Float_t*[tracks];
2808 for(Int_t i=0;i<tracks;i++) {
2814 AliPoints *points = 0;
2817 Int_t chunk=nhits/4+1;
2819 // Loop over all the hits and store their position
2821 ahit = FirstHit2(-1);
2823 trk=ahit->GetTrack();
2824 if(ntrk[trk]==limi[trk]) {
2826 // Initialise a new track
2827 fp=new Float_t[3*(limi[trk]+chunk)];
2829 memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2830 delete [] coor[trk];
2837 fp[3*ntrk[trk] ] = ahit->X();
2838 fp[3*ntrk[trk]+1] = ahit->Y();
2839 fp[3*ntrk[trk]+2] = ahit->Z();
2847 for(trk=0; trk<tracks; ++trk) {
2849 points = new AliPoints();
2850 points->SetMarkerColor(GetMarkerColor());
2851 points->SetMarkerSize(GetMarkerSize());
2852 points->SetDetector(this);
2853 points->SetParticle(trk);
2854 points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2855 fPoints->AddAt(points,trk);
2856 delete [] coor[trk];
2866 //_____________________________________________________________________________
2867 void AliTPC::LoadPoints3(Int_t)
2870 // Store x, y, z of all hits in memory
2871 // - only intersection point with pad row
2872 if (fTrackHits == 0) return;
2874 Int_t nhits = fTrackHits->GetEntriesFast();
2875 if (nhits == 0) return;
2876 Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2877 if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2878 fPoints->Expand(2*tracks);
2881 Int_t *ntrk=new Int_t[tracks];
2882 Int_t *limi=new Int_t[tracks];
2883 Float_t **coor=new Float_t*[tracks];
2884 for(Int_t i=0;i<tracks;i++) {
2890 AliPoints *points = 0;
2893 Int_t chunk=nhits/4+1;
2895 // Loop over all the hits and store their position
2897 ahit = FirstHit2(-1);
2898 //for (Int_t hit=0;hit<nhits;hit++) {
2902 // ahit = (AliHit*)fHits->UncheckedAt(hit);
2903 trk=ahit->GetTrack();
2904 Float_t x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2905 Int_t index[3]={1,((AliTPChit*)ahit)->fSector,0};
2906 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
2907 if (currentrow!=lastrow){
2908 lastrow = currentrow;
2909 //later calculate intersection point
2910 if(ntrk[trk]==limi[trk]) {
2912 // Initialise a new track
2913 fp=new Float_t[3*(limi[trk]+chunk)];
2915 memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2916 delete [] coor[trk];
2923 fp[3*ntrk[trk] ] = ahit->X();
2924 fp[3*ntrk[trk]+1] = ahit->Y();
2925 fp[3*ntrk[trk]+2] = ahit->Z();
2932 for(trk=0; trk<tracks; ++trk) {
2934 points = new AliPoints();
2935 points->SetMarkerColor(GetMarkerColor()+1);
2936 points->SetMarkerStyle(5);
2937 points->SetMarkerSize(0.2);
2938 points->SetDetector(this);
2939 points->SetParticle(trk);
2940 // points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20);
2941 points->SetPolyMarker(ntrk[trk],coor[trk],30);
2942 fPoints->AddAt(points,tracks+trk);
2943 delete [] coor[trk];
2954 void AliTPC::FindTrackHitsIntersection(TClonesArray * /*arr*/)
2958 //fill clones array with intersection of current point with the
2963 const Int_t kcmaxhits=30000;
2964 AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
2965 AliTPCFastVector & xxx = *xxxx;
2966 Int_t maxhits = kcmaxhits;
2969 AliTPChit * tpcHit=0;
2970 tpcHit = (AliTPChit*)FirstHit2(-1);
2971 Int_t currentIndex=0;
2972 Int_t lastrow=-1; //last writen row
2975 if (tpcHit==0) continue;
2976 sector=tpcHit->fSector; // sector number
2977 ipart=tpcHit->Track();
2981 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
2982 Int_t index[3]={1,sector,0};
2983 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
2984 if (currentrow<0) continue;
2985 if (lastrow<0) lastrow=currentrow;
2986 if (currentrow==lastrow){
2987 if ( currentIndex>=maxhits){
2989 xxx.ResizeTo(4*maxhits);
2991 xxx(currentIndex*4)=x[0];
2992 xxx(currentIndex*4+1)=x[1];
2993 xxx(currentIndex*4+2)=x[2];
2994 xxx(currentIndex*4+3)=tpcHit->fQ;
2998 if (currentIndex>2){
3010 for (Int_t index=0;index<currentIndex;index++){
3012 x=x2=x3=x4=xxx(index*4);
3020 sumy+=xxx(index*4+1);
3021 sumxy+=xxx(index*4+1)*x;
3022 sumx2y+=xxx(index*4+1)*x2;
3023 sumz+=xxx(index*4+2);
3024 sumxz+=xxx(index*4+2)*x;
3025 sumx2z+=xxx(index*4+2)*x2;
3026 sumq+=xxx(index*4+3);
3028 Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
3029 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
3030 sumx2*(sumx*sumx3-sumx2*sumx2);
3032 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
3033 sumx2*(sumxy*sumx3-sumx2y*sumx2);
3034 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
3035 sumx2*(sumxz*sumx3-sumx2z*sumx2);
3037 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
3038 sumx2*(sumx*sumx2y-sumx2*sumxy);
3039 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
3040 sumx2*(sumx*sumx2z-sumx2*sumxz);
3042 Float_t y=detay/det+centralPad;
3043 Float_t z=detaz/det;
3044 Float_t by=detby/det; //y angle
3045 Float_t bz=detbz/det; //z angle
3046 sumy/=Float_t(currentIndex);
3047 sumz/=Float_t(currentIndex);
3049 AliComplexCluster cl;
3055 cl.fTracks[0]=ipart;
3057 AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
3058 if (row!=0) row->InsertCluster(&cl);
3061 } //end of calculating cluster for given row
3063 } // end of loop over hits
3067 //_______________________________________________________________________________
3068 void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
3070 // produces rec points from digits and writes them on the root file
3071 // AliTPCclusters.root
3073 TDirectory *cwd = gDirectory;
3076 AliTPCParamSR *dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60");
3078 printf("You are running 2 pad-length geom hits with 3 pad-length geom digits\n");
3080 dig = new AliTPCParamSR();
3084 dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
3087 printf("No TPC parameters found\n");
3092 cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<<endl;
3094 if(!gSystem->Getenv("CONFIG_FILE")){
3095 out=TFile::Open("AliTPCclusters.root","recreate");
3101 tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX");
3102 tmp2=gSystem->Getenv("CONFIG_OUTDIR");
3103 sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2);
3104 out=TFile::Open(tmp3,"recreate");
3108 cout<<"AliTPC::Digits2Reco - determination of rec points begins"<<endl;
3111 for(Int_t iev=firstevent;iev<lastevent+1;iev++){
3113 printf("Processing event %d\n",iev);
3114 Digits2Clusters(iev);
3116 cout<<"AliTPC::Digits2Reco - determination of rec points ended"<<endl;
3125 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
3128 fLoader = new AliTPCLoader(GetName(),topfoldername);
3132 ////////////////////////////////////////////////////////////////////////
3133 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
3135 // load TPC paarmeters from a given file or create new if the object
3136 // is not found there
3137 // 12/05/2003 This method should be moved to the AliTPCLoader
3138 // and one has to decide where to store the TPC parameters
3141 sprintf(paramName,"75x40_100x60_150x60");
3142 AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
3144 cout<<"TPC parameters "<<paramName<<" found."<<endl;
3146 cerr<<"TPC parameters not found. Create new (they may be incorrect)."
3148 paramTPC = new AliTPCParamSR;
3152 // the older version of parameters can be accessed with this code.
3153 // In some cases, we have old parameters saved in the file but
3154 // digits were created with new parameters, it can be distinguish
3155 // by the name of TPC TreeD. The code here is just for the case
3156 // we would need to compare with old data, uncomment it if needed.
3158 // char paramName[50];
3159 // sprintf(paramName,"75x40_100x60");
3160 // AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
3162 // cout<<"TPC parameters "<<paramName<<" found."<<endl;
3164 // sprintf(paramName,"75x40_100x60_150x60");
3165 // paramTPC=(AliTPCParam*)in->Get(paramName);
3167 // cout<<"TPC parameters "<<paramName<<" found."<<endl;
3169 // cerr<<"TPC parameters not found. Create new (they may be incorrect)."
3171 // paramTPC = new AliTPCParamSR;
3179 //____________________________________________________________________________
3180 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
3183 // Parametrised error of the cluster reconstruction (pad direction)
3186 const Float_t kArphi=0.41818e-2;
3187 const Float_t kBrphi=0.17460e-4;
3188 const Float_t kCrphi=0.30993e-2;
3189 const Float_t kDrphi=0.41061e-3;
3191 pt=TMath::Abs(pt)*1000.;
3193 tgl=TMath::Abs(tgl);
3194 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
3195 if (s<0.4e-3) s=0.4e-3;
3196 s*=1.3; //Iouri Belikov
3202 //____________________________________________________________________________
3203 Double_t SigmaZ2(Double_t r, Double_t tgl)
3206 // Parametrised error of the cluster reconstruction (drift direction)
3209 const Float_t kAz=0.39614e-2;
3210 const Float_t kBz=0.22443e-4;
3211 const Float_t kCz=0.51504e-1;
3214 tgl=TMath::Abs(tgl);
3215 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
3216 if (s<0.4e-3) s=0.4e-3;
3217 s*=1.3; //Iouri Belikov