X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPC.cxx;h=8148aacaf8ec9e02569bc7ea1f128dfe9ca904e4;hb=f8cf550cef855758e7f3fabb986a54910576fd35;hp=3fd91e2fb7a0c08aa2e93f5d7760e74aa107755e;hpb=1578254f491848aff213be55a70a61d195dd001a;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPC.cxx b/TPC/AliTPC.cxx index 3fd91e2fb7a..8148aacaf8e 100644 --- a/TPC/AliTPC.cxx +++ b/TPC/AliTPC.cxx @@ -1,3 +1,111 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* +$Log$ +Revision 1.36 2001/05/16 14:57:25 alibrary +New files for folders and Stack + +Revision 1.35 2001/05/08 16:02:22 kowal2 +Updated material specifications + +Revision 1.34 2001/05/08 15:00:15 hristov +Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov) + +Revision 1.33 2001/04/03 12:40:43 kowal2 +Removed printouts + +Revision 1.32 2001/03/12 17:47:36 hristov +Changes needed on Sun with CC 5.0 + +Revision 1.31 2001/03/12 08:21:50 kowal2 +Corrected C++ bug in the material definitions + +Revision 1.30 2001/03/01 17:34:47 kowal2 +Correction due to the accuracy problem + +Revision 1.29 2001/02/28 16:34:40 kowal2 +Protection against nonphysical values of the avalanche size, +10**6 is the maximum + +Revision 1.28 2001/01/26 19:57:19 hristov +Major upgrade of AliRoot code + +Revision 1.27 2001/01/13 17:29:33 kowal2 +Sun compiler correction + +Revision 1.26 2001/01/10 07:59:43 kowal2 +Corrections to load points from the noncompressed hits. + +Revision 1.25 2000/11/02 07:25:31 kowal2 +Changes due to the new hit structure. +Memory leak removed. + +Revision 1.24 2000/10/05 16:06:09 kowal2 +Forward declarations. Changes due to a new class AliComplexCluster. + +Revision 1.23 2000/10/02 21:28:18 fca +Removal of useless dependecies via forward declarations + +Revision 1.22 2000/07/10 20:57:39 hristov +Update of TPC code and macros by M.Kowalski + +Revision 1.19.2.4 2000/06/26 07:39:42 kowal2 +Changes to obey the coding rules + +Revision 1.19.2.3 2000/06/25 08:38:41 kowal2 +Splitted from AliTPCtracking + +Revision 1.19.2.2 2000/06/16 12:59:28 kowal2 +Changed parameter settings + +Revision 1.19.2.1 2000/06/09 07:15:07 kowal2 + +Defaults loaded automatically (hard-wired) +Optional parameters can be set via macro called in the constructor + +Revision 1.19 2000/04/18 19:00:59 fca +Small bug fixes to TPC files + +Revision 1.18 2000/04/17 09:37:33 kowal2 +removed obsolete AliTPCDigitsDisplay.C + +Revision 1.17.2.2 2000/04/10 08:15:12 kowal2 + +New, experimental data structure from M. Ivanov +New tracking algorithm +Different pad geometry for different sectors +Digitization rewritten + +Revision 1.17.2.1 2000/04/10 07:56:53 kowal2 +Not used anymore - removed + +Revision 1.17 2000/01/19 17:17:30 fca +Introducing a list of lists of hits -- more hits allowed for detector now + +Revision 1.16 1999/11/05 09:29:23 fca +Accept only signals > 0 + +Revision 1.15 1999/10/08 06:26:53 fca +Removed ClustersIndex - not used anymore + +Revision 1.14 1999/09/29 09:24:33 fca +Introduction of the Copyright and cvs Log + +*/ + /////////////////////////////////////////////////////////////////////////////// // // // Time Projection Chamber // @@ -14,6 +122,8 @@ // // /////////////////////////////////////////////////////////////////////////////// +// + #include #include #include @@ -24,16 +134,38 @@ #include #include "TParticle.h" #include "AliTPC.h" +#include #include "AliRun.h" #include +#include #include #include "AliMC.h" +#include "AliMagF.h" + -//MI change -#include "AliTPCParam.h" -#include "AliTPCD.h" +#include "AliTPCParamSR.h" #include "AliTPCPRF2D.h" #include "AliTPCRF1D.h" +#include "AliDigits.h" +#include "AliSimDigits.h" +#include "AliTPCTrackHits.h" +#include "AliPoints.h" +#include "AliArrayBranch.h" + + +#include "AliTPCDigitsArray.h" +#include "AliComplexCluster.h" +#include "AliClusters.h" +#include "AliTPCClustersRow.h" +#include "AliTPCClustersArray.h" + +#include "AliTPCcluster.h" +#include "AliTPCclusterer.h" +#include "AliTPCtracker.h" + +#include +#include + ClassImp(AliTPC) @@ -45,16 +177,16 @@ AliTPC::AliTPC() // Default constructor // fIshunt = 0; - fClusters = 0; fHits = 0; fDigits = 0; - fTracks = 0; fNsectors = 0; - fNtracks = 0; - fNclusters= 0; //MI changes - fDigParam= new AliTPCD(); - fDigits = fDigParam->GetArray(); + fDigitsArray = 0; + fClustersArray = 0; + fTPCParam=0; + fTrackHits = 0; + fHitType = 2; + fTPCParam = 0; } //_____________________________________________________________________________ @@ -68,24 +200,38 @@ AliTPC::AliTPC(const char *name, const char *title) // // Initialise arrays of hits and digits fHits = new TClonesArray("AliTPChit", 176); - // fDigits = new TClonesArray("AliTPCdigit",10000); - //MI change - fDigParam= new AliTPCD; - fDigits = fDigParam->GetArray(); + gAlice->AddHitList(fHits); + //MI change + fDigitsArray = 0; + fClustersArray= 0; + // + fTrackHits = new AliTPCTrackHits; //MI - 13.09.2000 + fTrackHits->SetHitPrecision(0.002); + fTrackHits->SetStepPrecision(0.003); + fTrackHits->SetMaxDistance(100); + + fHitType = 2; // // Initialise counters - fClusters = 0; - fTracks = 0; - fNsectors = 72; - fNtracks = 0; - fNclusters= 0; - fDigitsIndex = new Int_t[fNsectors+1]; - fClustersIndex = new Int_t[fNsectors+1]; + fNsectors = 0; + // fIshunt = 0; // // Initialise color attributes SetMarkerColor(kYellow); + + // + // Set TPC parameters + // + + if (!strcmp(title,"Default")) { + fTPCParam = new AliTPCParamSR; + } else { + cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n"; + fTPCParam=0; + } + } //_____________________________________________________________________________ @@ -94,126 +240,115 @@ AliTPC::~AliTPC() // // TPC destructor // + fIshunt = 0; delete fHits; delete fDigits; - delete fClusters; - delete fTracks; - delete fDigParam; - if (fDigitsIndex) delete [] fDigitsIndex; - if (fClustersIndex) delete [] fClustersIndex; + delete fTPCParam; + delete fTrackHits; //MI 15.09.2000 } -//_____________________________________________________________________________ -void AliTPC::AddCluster(Float_t *hits, Int_t *tracks) -{ - // - // Add a simulated cluster to the list - // - if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000); - TClonesArray &lclusters = *fClusters; - new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks); -} - -//_____________________________________________________________________________ -void AliTPC::AddCluster(const AliTPCcluster &c) -{ - // - // Add a simulated cluster copy to the list - // - if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000); - TClonesArray &lclusters = *fClusters; - new(lclusters[fNclusters++]) AliTPCcluster(c); -} - -//_____________________________________________________________________________ -void AliTPC::AddDigit(Int_t *tracks, Int_t *digits) -{ - // - // Add a TPC digit to the list - // - // TClonesArray &ldigits = *fDigits; - //MI change - TClonesArray &ldigits = *fDigParam->GetArray(); - new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits); -} - //_____________________________________________________________________________ void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits) { // // Add a hit to the list // - TClonesArray &lhits = *fHits; - new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits); + // TClonesArray &lhits = *fHits; + // new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits); + if (fHitType&1){ + TClonesArray &lhits = *fHits; + new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits); + } + if (fHitType&2) + AddHit2(track,vol,hits); } -//_____________________________________________________________________________ -void AliTPC::AddTrack(Float_t *hits) -{ - // - // Add a track to the list of tracks - // - TClonesArray <racks = *fTracks; - new(ltracks[fNtracks++]) AliTPCtrack(hits); -} - -//_____________________________________________________________________________ -void AliTPC::AddTrack(const AliTPCtrack& t) -{ - // - // Add a track copy to the list of tracks - // - if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000); - TClonesArray <racks = *fTracks; - new(ltracks[fNtracks++]) AliTPCtrack(t); -} - //_____________________________________________________________________________ void AliTPC::BuildGeometry() { + // // Build TPC ROOT TNode geometry for the event display // - TNode *Node, *Top; + TNode *nNode, *nTop; TTUBS *tubs; Int_t i; const int kColorTPC=19; - char name[5], title[20]; + char name[5], title[25]; const Double_t kDegrad=TMath::Pi()/180; - const Double_t loAng=30; - const Double_t hiAng=15; - const Int_t nLo = Int_t (360/loAng+0.5); - const Int_t nHi = Int_t (360/hiAng+0.5); - const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad); - const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad); + const Double_t kRaddeg=180./TMath::Pi(); + + + Float_t innerOpenAngle = fTPCParam->GetInnerAngle(); + Float_t outerOpenAngle = fTPCParam->GetOuterAngle(); + + Float_t innerAngleShift = fTPCParam->GetInnerAngleShift(); + Float_t outerAngleShift = fTPCParam->GetOuterAngleShift(); + + Int_t nLo = fTPCParam->GetNInnerSector()/2; + Int_t nHi = fTPCParam->GetNOuterSector()/2; + + const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg); + const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg); + const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg); + const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg); + + + const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad); + const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad); + + Double_t rl,ru; + + // // Get ALICE top node - Top=gAlice->GetGeometry()->GetNode("alice"); // - // Inner sectors + + nTop=gAlice->GetGeometry()->GetNode("alice"); + + // inner sectors + + rl = fTPCParam->GetInnerRadiusLow(); + ru = fTPCParam->GetInnerRadiusUp(); + + for(i=0;iSetNumberOfDivisions(1); - Top->cd(); - Node = new TNode(name,title,name,0,0,0,""); - Node->SetLineColor(kColorTPC); - fNodes->Add(Node); + nTop->cd(); + nNode = new TNode(name,title,name,0,0,0,""); + nNode->SetLineColor(kColorTPC); + fNodes->Add(nNode); } + // Outer sectors + + rl = fTPCParam->GetOuterRadiusLow(); + ru = fTPCParam->GetOuterRadiusUp(); + for(i=0;iSetNumberOfDivisions(1); - Top->cd(); - Node = new TNode(name,title,name,0,0,0,""); - Node->SetLineColor(kColorTPC); - fNodes->Add(Node); + nTop->cd(); + nNode = new TNode(name,title,name,0,0,0,""); + nNode->SetLineColor(kColorTPC); + fNodes->Add(nNode); } -} + +} + //_____________________________________________________________________________ Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) { @@ -224,654 +359,435 @@ Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) return 9999; } -//_____________________________________________________________________________ -//const int MAX_CLUSTER=nrow_low+nrow_up; -const int MAX_CLUSTER=200; -const int S_MAXSEC=24; -const int L_MAXSEC=48; -const int ROWS_TO_SKIP=21; -const Float_t MAX_CHI2=12.; - - -//_____________________________________________________________________________ -static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt) -{ - // - // Calculate spread in Y - // - pt=TMath::Abs(pt)*1000.; - Double_t x=r/pt; - tgl=TMath::Abs(tgl); - Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x; - if (s<0.4e-3) s=0.4e-3; - return s; +void AliTPC::Clusters2Tracks(TFile *of) { + //----------------------------------------------------------------- + // This is a track finder. + //----------------------------------------------------------------- + AliTPCtracker tracker(fTPCParam); + tracker.Clusters2Tracks(gFile,of); } //_____________________________________________________________________________ -static Double_t SigmaZ2(Double_t r, Double_t tgl) +void AliTPC::CreateMaterials() { - // - // Calculate spread in Z - // - tgl=TMath::Abs(tgl); - Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl; - if (s<0.4e-3) s=0.4e-3; - return s; -} + //----------------------------------------------- + // Create Materials for for TPC simulations + //----------------------------------------------- -//_____________________________________________________________________________ -inline Double_t f1(Double_t x1,Double_t y1, //C - Double_t x2,Double_t y2, - Double_t x3,Double_t y3) -{ - // - // Function f1 - // - Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); - Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- - (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); - Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- - (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); + //----------------------------------------------------------------- + // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl + //----------------------------------------------------------------- - Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); + Int_t iSXFLD=gAlice->Field()->Integ(); + Float_t sXMGMX=gAlice->Field()->Max(); - return -xr*yr/sqrt(xr*xr+yr*yr); -} + Float_t amat[5]; // atomic numbers + Float_t zmat[5]; // z + Float_t wmat[5]; // proportions + Float_t density; + Float_t apure[2]; -//_____________________________________________________________________________ -inline Double_t f2(Double_t x1,Double_t y1, //eta=C*x0 - Double_t x2,Double_t y2, - Double_t x3,Double_t y3) -{ - // - // Function f2 - // - Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); - Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- - (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); - Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- - (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); - Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); + //***************** Gases ************************* - return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); -} + //------------------------------------------------- + // pure gases + //------------------------------------------------- -//_____________________________________________________________________________ -inline Double_t f3(Double_t x1,Double_t y1, //tgl - Double_t x2,Double_t y2, - Double_t z1,Double_t z2) -{ - // - // Function f3 - // - return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); -} + // Neon -//_____________________________________________________________________________ -static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec, - int s, int ri, int rf=0) -{ - // - // Propagate track - // - int try_again=ROWS_TO_SKIP; - Double_t alpha=sec->GetAlpha(); - int ns=int(2*TMath::Pi()/alpha)+1; - for (int nr=ri; nr>=rf; nr--) { - Double_t x=sec[s].GetX(nr), ymax=sec[s].GetMaxY(nr); - if (!t.PropagateTo(x)) return -1; + amat[0]= 20.18; + zmat[0]= 10.; + density = 0.0009; + + apure[0]=amat[0]; - const AliTPCcluster *cl=0; - Double_t max_chi2=MAX_CHI2; - const AliTPCRow& row=sec[s][nr]; - Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); - Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl()); - Double_t road=3.*sqrt(t.GetSigmaY2() + 4*sy2), y=t.GetY(), z=t.GetZ(); + AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.); - if (road>30) { - if (t>3) cerr<fY > y+road) break; - if (c->IsUsed()) continue; - if ((c->fZ - z)*(c->fZ - z) > 9.*(t.GetSigmaZ2() + 4*sz2)) continue; - Double_t chi2=t.GetPredictedChi2(c); - if (chi2 > max_chi2) continue; - max_chi2=chi2; - cl=c; - } - } - if (cl) { - t.Update(cl,max_chi2); - try_again=ROWS_TO_SKIP; - } else { - if (try_again==0) break; - if (y > ymax) { - s = (s+1) % ns; - if (!t.Rotate(alpha)) return -1; - } else if (y <-ymax) { - s = (s-1+ns) % ns; - if (!t.Rotate(-alpha)) return -1; - } - try_again--; - } - } + amat[0]= 39.948; + zmat[0]= 18.; + density = 0.001782; - return s; -} + apure[1]=amat[0]; + AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.); + -//_____________________________________________________________________________ -static void MakeSeeds(TObjArray& seeds,const AliTPCSector* sec,int i1,int i2, -const AliTPCParam *p) -{ - // - // Find seed for tracking - // - TMatrix C(5,5); TVector x(5); - int max_sec=L_MAXSEC/2; - for (int ns=0; nsfY, z1=r1[is]->fZ; - for (int js=0; js < nl+nm+nu; js++) { - const AliTPCcluster *cl; - Double_t cs,sn; - int ks; - - if (jsfY, z2=cl->fZ; - //if (z1*z2 < 0) continue; - //if (TMath::Abs(z1) < TMath::Abs(z2)) continue; - - Double_t tmp= x2*cs+y2*sn; - y2 =-x2*sn+y2*cs; - x2=tmp; - - x(0)=y1; - x(1)=z1; - x(2)=f1(x1,y1,x2,y2,0.,0.); - x(3)=f2(x1,y1,x2,y2,0.,0.); - x(4)=f3(x1,y1,x2,y2,z1,z2); - - if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue; - - if (TMath::Abs(x(4)) > 1.2) continue; - - Double_t a=asin(x(3)); - Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3))); - if (TMath::Abs(zv)>33.) continue; - - TMatrix X(6,6); X=0.; - X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2; - X(2,2)=cl->fSigmaY2; X(3,3)=cl->fSigmaZ2; - X(4,4)=3./12.; X(5,5)=3./12.; - TMatrix F(5,6); F.UnitMatrix(); - Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1)); - F(2,0)=(f1(x1,y1+sy,x2,y2,0.,0.)-x(2))/sy; - F(2,2)=(f1(x1,y1,x2,y2+sy,0.,0.)-x(2))/sy; - F(2,4)=(f1(x1,y1,x2,y2,0.,0.+sy)-x(2))/sy; - F(3,0)=(f2(x1,y1+sy,x2,y2,0.,0.)-x(3))/sy; - F(3,2)=(f2(x1,y1,x2,y2+sy,0.,0.)-x(3))/sy; - F(3,4)=(f2(x1,y1,x2,y2,0.,0.+sy)-x(3))/sy; - F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy; - F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz; - F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy; - F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz; - F(4,4)=0; - F(3,3)=0; - - TMatrix t(F,TMatrix::kMult,X); - C.Mult(t,TMatrix(TMatrix::kTransposed,F)); + //-------------------------------------------------------------- + // gases - compounds + //-------------------------------------------------------------- - TrackSeed *track=new TrackSeed(*(r1[is]),x,C,p); - int rc=FindProlongation(*track,sec,ns,i1-1,i2); - if (rc<0 || *track<(i1-i2)/2) delete track; - else seeds.AddLast(track); - } - } - } -} + Float_t amol[3]; -//_____________________________________________________________________________ -void AliTPC::Clusters2Tracks() -{ - // - // TPC Track finder from clusters. - // - if (!fClusters) return; + // CO2 - AliTPCParam *p=&fDigParam->GetParam(); - Int_t nrow_low=p->GetNRowLow(); - Int_t nrow_up=p->GetNRowUp(); + amat[0]=12.011; + amat[1]=15.9994; - AliTPCSSector ssec[S_MAXSEC/2]; - for (int i=0; iGetEntriesFast(); - while (ncl--) { - AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl); + density=0.001977; - int sec=int(c->fSector)-1, row=int(c->fPadRow)-1; - - if (sec<24) { - if (row<0 || row>nrow_low) {cerr<<"low !!!"<nrow_up ) {cerr<<"up !!!"< 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); - if (alpha < 0. ) alpha += 2.*TMath::Pi(); - int ns=int(alpha/lsec->GetAlpha() + 0.5); - - Double_t x=t.GetX(); - int nr; - if (xGetPadRowRadiiUp(nrow_up-1-4-7)) nr=nrow_up-1-4-8; - else if (xGetPadRowRadiiUp(nrow_up-1-7)) nr=nrow_up-1-8; - else {cerr<Field()->Integ(); - Float_t SXMGMX=gAlice->Field()->Max(); - - Float_t absl, radl, a, d, z; - Float_t dg; - Float_t x0ne; - Float_t buf[1]; - Int_t nbuf; - - // --- Methane (CH4) --- - Float_t am[2] = { 12.,1. }; - Float_t zm[2] = { 6.,1. }; - Float_t wm[2] = { 1.,4. }; - Float_t dm = 7.17e-4; - // --- The Neon CO2 90/10 mixture --- - Float_t ag[2] = { 20.18 }; - Float_t zg[2] = { 10. }; - Float_t wg[2] = { .8,.2 }; - Float_t dne = 9e-4; // --- Neon density in g/cm3 --- - - // --- Mylar (C5H4O2) --- - Float_t amy[3] = { 12.,1.,16. }; - Float_t zmy[3] = { 6.,1.,8. }; - Float_t wmy[3] = { 5.,4.,2. }; - Float_t dmy = 1.39; - // --- CO2 --- - Float_t ac[2] = { 12.,16. }; - Float_t zc[2] = { 6.,8. }; - Float_t wc[2] = { 1.,2. }; - Float_t dc = .001977; - // --- Carbon density and radiation length --- - Float_t densc = 2.265; - Float_t radlc = 18.8; - // --- Silicon --- - Float_t asi = 28.09; - Float_t zsi = 14.; - Float_t desi = 2.33; - Float_t radsi = 9.36; - - // --- Define the various materials for GEANT --- - AliMaterial(0, "Al $", 26.98, 13., 2.7, 8.9, 37.2); - x0ne = 28.94 / dne; - AliMaterial(1, "Ne $", 20.18, 10., dne, x0ne, 999.); - - // -- Methane, defined by the proportions of atoms - - AliMixture(2, "Methane$", am, zm, dm, -2, wm); - - // --- CO2, defined by the proportion of atoms - - AliMixture(7, "CO2$", ac, zc, dc, -2, wc); - - // -- Get A,Z etc. for CO2 - - char namate[21]; - pMC->Gfmate((*fIdmate)[7], namate, a, z, d, radl, absl, buf, nbuf); - ag[1] = a; - zg[1] = z; - dg = dne * .9 + dc * .1; - - // -- Create Ne/CO2 90/10 mixture - - AliMixture(3, "Gas-mixt $", ag, zg, dg, 2, wg); - AliMixture(4, "Gas-mixt $", ag, zg, dg, 2, wg); - - AliMaterial(5, "G10$", 20., 10., 1.7, 19.4, 999.); - AliMixture(6, "Mylar$", amy, zmy, dmy, -3, wmy); - - a = ac[0]; - z = zc[0]; - AliMaterial(8, "Carbon", a, z, densc, radlc, 999.); - - AliMaterial(9, "Silicon", asi, zsi, desi, radsi, 999.); - AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500.); - - AliMedium(400, "Al wall$", 0, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(402, "Gas mix1$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01); - AliMedium(403, "Gas mix2$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01); - AliMedium(404, "Gas mix3$", 4, 1, ISXFLD, SXMGMX, 10., .01,.1, .001, .01); - AliMedium(405, "G10 pln$", 5, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 ); - AliMedium(406, "Mylar $", 6, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01); - AliMedium(407, "CO2 $", 7, 0, ISXFLD, SXMGMX, 10., .01,.1, .01, .01); - AliMedium(408, "Carbon $", 8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 ); - AliMedium(409, "Silicon$", 9, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 ); - AliMedium(499, "Air gap$", 99, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 ); -} + amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1]; -//_____________________________________________________________________________ -struct Bin { - const AliTPCdigit *dig; - int idx; - Bin() {dig=0; idx=-1;} -}; - -struct PreCluster : public AliTPCcluster { - const AliTPCdigit* summit; - int idx; - int cut; - int npeaks; - PreCluster(); -}; -PreCluster::PreCluster() : AliTPCcluster() {cut=npeaks=0;} + AliMixture(11,"CF4",amat,zmat,density,-2,wmat); -//_____________________________________________________________________________ -static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c) -{ - // - // Find clusters - // - Bin& b=bins[i][j]; - double q=double(b.dig->fSignal); + // CH4 - if (q<0) { // digit is at the edge of the pad row - q=-q; - c.cut=1; - } - if (b.idx >= 0 && b.idx != c.idx) { - c.idx=b.idx; - c.npeaks++; - } - - if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig; - - c.fY += i*q; - c.fZ += j*q; - c.fSigmaY2 += i*i*q; - c.fSigmaZ2 += j*j*q; - c.fQ += q; + amat[0]=12.011; + amat[1]=1.; - b.dig = 0; b.idx = c.idx; - - if (bins[i-1][j].dig) FindCluster(i-1,j,bins,c); - if (bins[i][j-1].dig) FindCluster(i,j-1,bins,c); - if (bins[i+1][j].dig) FindCluster(i+1,j,bins,c); - if (bins[i][j+1].dig) FindCluster(i,j+1,bins,c); - -} + zmat[0]=6.; + zmat[1]=1.; -//_____________________________________________________________________________ -void AliTPC::Digits2Clusters() -{ - // - // simple TPC cluster finder from digits. - // - // - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - - const Int_t MAX_PAD=200+2, MAX_BUCKET=MAXTPCTBK+2; - const Int_t Q_min=60; - const Int_t THRESHOLD=20; - - TTree *t=(TTree*)gDirectory->Get("TreeD0_Param1"); - t->GetBranch("Digits")->SetAddress(&fDigits); - Int_t sectors_by_rows=(Int_t)t->GetEntries(); - - int ncls=0; - - for (Int_t n=0; nGetEvent(n)) continue; - Bin bins[MAX_PAD][MAX_BUCKET]; - AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0); - Int_t nsec=dig->fSector, nrow=dig->fPadRow; - Int_t ndigits=fDigits->GetEntriesFast(); - - int npads; int sign_z; - if (nsec<25) { - sign_z=(nsec<13) ? 1 : -1; - npads=fTPCParam->GetNPadsLow(nrow-1); - } else { - sign_z=(nsec<49) ? 1 : -1; - npads=fTPCParam->GetNPadsUp(nrow-1); - } + wmat[0]=1.; + wmat[1]=4.; + + density=0.000717; + + amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1]; + + AliMixture(12,"CH4",amat,zmat,density,-2,wmat); + + //---------------------------------------------------------------- + // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds + //---------------------------------------------------------------- + + char namate[21]; + density = 0.; + Float_t am=0; + Int_t nc; + Float_t rho,absl,X0,buf[1]; + Int_t nbuf; + Float_t a,z; + + for(nc = 0;ncUncheckedAt(ndig); - int i=dig->fPad, j=dig->fTime; - if (dig->fSignal >= THRESHOLD) bins[i][j].dig=dig; - if (i==1 || i==npads || j==1 || j==MAXTPCTBK) dig->fSignal*=-1; + // retrive material constants + + gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf); + + amat[nc] = a; + zmat[nc] = z; + + Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10; + + am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]); + density += fMixtProp[nc]*rho; // density of the mixture + } - - int ncl=0; - int i,j; - - for (i=1; iGetPadPitchWidth()* - fTPCParam->GetPadPitchWidth(); - if (s2 != 0.) c.fSigmaY2 *= 0.022*8*4; - - s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ; - c.fSigmaZ2 = s2 + 1./12.; - c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth(); - if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*4; - - c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth(); - c.fZ = fTPCParam->GetZWidth()*(c.fZ+1); - c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay - c.fZ = sign_z*(z_end - c.fZ); - //c.fZ += 0.023; - c.fSector=nsec; - c.fPadRow=nrow; - c.fTracks[0]=c.summit->fTracks[0]; - c.fTracks[1]=c.summit->fTracks[1]; - c.fTracks[2]=c.summit->fTracks[2]; - - if (c.cut) { - c.fSigmaY2 *= 25.; - c.fSigmaZ2 *= 4.; - } - - AddCluster(c); ncls++; ncl++; - } - } - - for (ndig=0; ndigUncheckedAt(ndig); - if (TMath::Abs(dig->fSignal) >= 0) - bins[dig->fPad][dig->fTime].dig=dig; - } - - for (i=1; i1) continue; //overlapped cluster - c.fY /= c.fQ; - c.fZ /= c.fQ; - - double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY; - c.fSigmaY2 = s2 + 1./12.; - c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()* - fTPCParam->GetPadPitchWidth(); - if (s2 != 0.) c.fSigmaY2 *= 0.022*4*0.6*4; - - s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ; - c.fSigmaZ2 = s2 + 1./12.; - c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth(); - if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*0.4; - - c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth(); - c.fZ = fTPCParam->GetZWidth()*(c.fZ+1); - c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay - c.fZ = sign_z*(z_end - c.fZ); - //c.fZ += 0.023; - c.fSector=nsec; - c.fPadRow=nrow; - c.fTracks[0]=c.summit->fTracks[0]; - c.fTracks[1]=c.summit->fTracks[1]; - c.fTracks[2]=c.summit->fTracks[2]; - - if (c.cut) { - c.fSigmaY2 *= 25.; - c.fSigmaZ2 *= 4.; - } - - if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;} - else { - new ((*fClusters)[c.idx]) AliTPCcluster(c); - } - } - } - - cerr<<"sector, row, digits, clusters: " - <Clear(); - - } -} -//_____________________________________________________________________________ -void AliTPC::ElDiff(Float_t *xyz) -{ - //-------------------------------------------------- - // calculates the diffusion of a single electron - //-------------------------------------------------- + // mixture proportions by weight! + + for(nc = 0;nc=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10; + + wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? + apure[nnc] : amol[nnc])/am; + + } + + // Drift gases 1 - nonsensitive, 2 - sensitive + + AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat); + AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat); + + + // Air + + amat[0] = 14.61; + zmat[0] = 7.3; + density = 0.001205; + + AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.); + + + //---------------------------------------------------------------------- + // solid materials + //---------------------------------------------------------------------- + + + // Kevlar C14H22O2N2 + + amat[0] = 12.011; + amat[1] = 1.; + amat[2] = 15.999; + amat[3] = 14.006; + + zmat[0] = 6.; + zmat[1] = 1.; + zmat[2] = 8.; + zmat[3] = 7.; + + wmat[0] = 14.; + wmat[1] = 22.; + wmat[2] = 2.; + wmat[3] = 2.; + + density = 1.45; + + AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat); + + // NOMEX + + amat[0] = 12.011; + amat[1] = 1.; + amat[2] = 15.999; + amat[3] = 14.006; + + zmat[0] = 6.; + zmat[1] = 1.; + zmat[2] = 8.; + zmat[3] = 7.; + + wmat[0] = 14.; + wmat[1] = 22.; + wmat[2] = 2.; + wmat[3] = 2.; + + density = 0.03; - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t driftl; - Float_t z0=xyz[2]; + AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat); - driftl=z_end-TMath::Abs(xyz[2]); + // Makrolon C16H18O3 - if(driftl<0.01) driftl=0.01; + amat[0] = 12.011; + amat[1] = 1.; + amat[2] = 15.999; - // check the attachment + zmat[0] = 6.; + zmat[1] = 1.; + zmat[2] = 8.; - driftl=TMath::Sqrt(driftl); + wmat[0] = 16.; + wmat[1] = 18.; + wmat[2] = 3.; + + density = 1.2; - // Float_t sig_t = driftl*diff_t; - //Float_t sig_l = driftl*diff_l; - Float_t sig_t = driftl*fTPCParam->GetDiffT(); - Float_t sig_l = driftl*fTPCParam->GetDiffL(); - xyz[0]=gRandom->Gaus(xyz[0],sig_t); - xyz[1]=gRandom->Gaus(xyz[1],sig_t); - xyz[2]=gRandom->Gaus(xyz[2],sig_l); + AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat); - if (TMath::Abs(xyz[2])>z_end){ - xyz[2]=TMath::Sign(z_end,z0); - } - if(xyz[2]*z0 < 0.){ - Float_t eps = 0.0001; - xyz[2]=TMath::Sign(eps,z0); - } + // Mylar C5H4O2 + + amat[0]=12.011; + amat[1]=1.; + amat[2]=15.9994; + + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; + + wmat[0]=5.; + wmat[1]=4.; + wmat[2]=2.; + + density = 1.39; + + AliMixture(37, "Mylar",amat,zmat,density,-3,wmat); + + // SiO2 - used later for the glass fiber + + amat[0]=28.086; + amat[1]=15.9994; + + zmat[0]=14.; + zmat[1]=8.; + + wmat[0]=1.; + wmat[1]=2.; + + + AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2) + + // Al + + amat[0] = 26.98; + zmat[0] = 13.; + + density = 2.7; + + AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.); + + // Si + + amat[0] = 28.086; + zmat[0] = 14.; + + density = 2.33; + + AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.); + + // Cu + + amat[0] = 63.546; + zmat[0] = 29.; + + density = 8.96; + + AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.); + + // Tedlar C2H3F + + amat[0] = 12.011; + amat[1] = 1.; + amat[2] = 18.998; + + zmat[0] = 6.; + zmat[1] = 1.; + zmat[2] = 9.; + + wmat[0] = 2.; + wmat[1] = 3.; + wmat[2] = 1.; + + density = 1.71; + + AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat); + + + // Plexiglas C5H8O2 + + amat[0]=12.011; + amat[1]=1.; + amat[2]=15.9994; + + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; + + wmat[0]=5.; + wmat[1]=8.; + wmat[2]=2.; + + density=1.18; + + AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat); + + // Epoxy - C14 H20 O3 + + + amat[0]=12.011; + amat[1]=1.; + amat[2]=15.9994; + + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; + + wmat[0]=14.; + wmat[1]=20.; + wmat[2]=3.; + + density=1.25; + + AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat); + + // Carbon + + amat[0]=12.011; + zmat[0]=6.; + density= 2.265; + + AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.); + + // get epoxy + + gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf); + + // Carbon fiber + + wmat[0]=0.644; // by weight! + wmat[1]=0.356; + + density=0.5*(1.25+2.265); + + AliMixture(47,"Cfiber",amat,zmat,density,2,wmat); + + // get SiO2 + + gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf); + + wmat[0]=0.725; // by weight! + wmat[1]=0.275; + + density=1.7; + + AliMixture(39,"G10",amat,zmat,density,2,wmat); + + + + + //---------------------------------------------------------- + // tracking media for gases + //---------------------------------------------------------- + + AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1); + AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); + AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); + AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001); + + //----------------------------------------------------------- + // tracking media for solids + //----------------------------------------------------------- + + AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + } + +void AliTPC::Digits2Clusters(TFile *of) +{ + //----------------------------------------------------------------- + // This is a simple cluster finder. + //----------------------------------------------------------------- + AliTPCclusterer::Digits2Clusters(fTPCParam,of); +} + +extern Double_t SigmaY2(Double_t, Double_t, Double_t); +extern Double_t SigmaZ2(Double_t, Double_t); //_____________________________________________________________________________ -void AliTPC::Hits2Clusters() +void AliTPC::Hits2Clusters(TFile *of) { //-------------------------------------------------------- // TPC simple cluster generator from hits @@ -882,13 +798,44 @@ void AliTPC::Hits2Clusters() //----------------------------------------------------------------- // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t sigma_rphi,sigma_z,cl_rphi,cl_z; + // Adopted to Marian's cluster data structure by I.Belikov, CERN, + // Jouri.Belikov@cern.ch + //---------------------------------------------------------------- + + ///////////////////////////////////////////////////////////////////////////// + // + //--------------------------------------------------------------------- + // ALICE TPC Cluster Parameters + //-------------------------------------------------------------------- + + + + // Cluster width in rphi + const Float_t kACrphi=0.18322; + const Float_t kBCrphi=0.59551e-3; + const Float_t kCCrphi=0.60952e-1; + // Cluster width in z + const Float_t kACz=0.19081; + const Float_t kBCz=0.55938e-3; + const Float_t kCCz=0.30428; + + TDirectory *savedir=gDirectory; + + if (!of->IsOpen()) { + cerr<<"AliTPC::Hits2Clusters(): output file not open !\n"; + return; + } + + if(fTPCParam == 0){ + printf("AliTPCParam MUST be created firstly\n"); + return; + } + + Float_t sigmaRphi,sigmaZ,clRphi,clZ; // TParticle *particle; // pointer to a given particle AliTPChit *tpcHit; // pointer to a sigle TPC hit - TClonesArray *Particles; //pointer to the particle list - Int_t sector,nhits; + Int_t sector; Int_t ipart; Float_t xyz[5]; Float_t pl,pt,tanth,rpad,ratio; @@ -898,24 +845,32 @@ void AliTPC::Hits2Clusters() // Get the access to the tracks //--------------------------------------------------------------- - TTree *TH = gAlice->TreeH(); - Stat_t ntracks = TH->GetEntries(); + TTree *tH = gAlice->TreeH(); + Stat_t ntracks = tH->GetEntries(); + + //Switch to the output file + of->cd(); + + fTPCParam->Write(fTPCParam->GetTitle()); + AliTPCClustersArray carray; + carray.Setup(fTPCParam); + carray.SetClusterType("AliTPCcluster"); + carray.MakeTree(); + + Int_t nclusters=0; //cluster counter //------------------------------------------------------------ - // Loop over all sectors (72 sectors) - // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0 - // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0 + // Loop over all sectors (72 sectors for 20 deg + // segmentation for both lower and upper sectors) + // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0 + // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0 // - // First cluster for sector 1 starts at "0" + // First cluster for sector 0 starts at "0" //------------------------------------------------------------ - - - fClustersIndex[0] = 0; - - // - for(Int_t isec=1;isecGetNSector();isec++){ //MI change - fTPCParam->AdjustAngles(isec,cph,sph); + fTPCParam->AdjustCosSin(isec,cph,sph); //------------------------------------------------------------ // Loop over tracks @@ -923,46 +878,61 @@ void AliTPC::Hits2Clusters() for(Int_t track=0;trackGetEvent(track); + tH->GetEvent(track); // - // Get number of the TPC hits and a pointer - // to the particles + // Get number of the TPC hits // - nhits=fHits->GetEntriesFast(); - Particles=gAlice->Particles(); + // nhits=fHits->GetEntriesFast(); // + + tpcHit = (AliTPChit*)FirstHit(-1); + // Loop over hits // - for(Int_t hit=0;hitUncheckedAt(hit); + // for(Int_t hit=0;hitUncheckedAt(hit); + + while(tpcHit){ + + if (tpcHit->fQ == 0.) { + tpcHit = (AliTPChit*) NextHit(); + continue; //information about track (I.Belikov) + } sector=tpcHit->fSector; // sector number - if(sector != isec) continue; //terminate iteration - ipart=tpcHit->fTrack; - particle=(TParticle*)Particles->UncheckedAt(ipart); + + + // if(sector != isec) continue; //terminate iteration + + if(sector != isec){ + tpcHit = (AliTPChit*) NextHit(); + continue; + } + ipart=tpcHit->Track(); + particle=gAlice->Particle(ipart); pl=particle->Pz(); pt=particle->Pt(); if(pt < 1.e-9) pt=1.e-9; tanth=pl/pt; tanth = TMath::Abs(tanth); - rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY); + rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y()); ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason - + // space-point resolutions - sigma_rphi=SigmaY2(rpad,tanth,pt); - sigma_z =SigmaZ2(rpad,tanth ); + sigmaRphi=SigmaY2(rpad,tanth,pt); + sigmaZ =SigmaZ2(rpad,tanth ); // cluster widths - cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio; - cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth; + clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio; + clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth; // temporary protection - if(sigma_rphi < 0.) sigma_rphi=0.4e-3; - if(sigma_z < 0.) sigma_z=0.4e-3; - if(cl_rphi < 0.) cl_rphi=2.5e-3; - if(cl_z < 0.) cl_z=2.5e-5; + if(sigmaRphi < 0.) sigmaRphi=0.4e-3; + if(sigmaZ < 0.) sigmaZ=0.4e-3; + if(clRphi < 0.) clRphi=2.5e-3; + if(clZ < 0.) clZ=2.5e-5; // @@ -970,47 +940,311 @@ void AliTPC::Hits2Clusters() // smearing --> rotate to the 1 (13) or to the 25 (49) sector, // then the inaccuracy in a X-Y plane is only along Y (pad row)! // - Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph; - Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph; - xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi)); // y - Double_t alpha=(sector<25) ? alpha_low : alpha_up; - if (TMath::Abs(xyz[0]/xprim) > TMath::Tan(0.5*alpha)) xyz[0]=yprim; - xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z - if (TMath::Abs(xyz[1]) > 250) xyz[1]=tpcHit->fZ; - xyz[2]=tpcHit->fQ; // q - xyz[3]=sigma_rphi; // fSigmaY2 - xyz[4]=sigma_z; // fSigmaZ2 - - //find row number - //MI we must change - Int_t row = fTPCParam->GetPadRow(sector,xprim) ; - // and finally add the cluster - Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, row+1}; - AddCluster(xyz,tracks); - + Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph; + Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph; + xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y + Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ? + fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle(); + Float_t ymax=xprim*TMath::Tan(0.5*alpha); + if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim; + xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z + if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z(); + xyz[2]=sigmaRphi; // fSigmaY2 + xyz[3]=sigmaZ; // fSigmaZ2 + xyz[4]=tpcHit->fQ; // q + + AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow); + if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow); + + Int_t tracks[3]={tpcHit->Track(), -1, -1}; + AliTPCcluster cluster(tracks,xyz); + + clrow->InsertCluster(&cluster); nclusters++; + + tpcHit = (AliTPChit*)NextHit(); + + } // end of loop over hits - } // end of loop over tracks - - fClustersIndex[isec] = fNclusters; // update clusters index - - } // end of loop over sectors + + } // end of loop over tracks + + Int_t nrows=fTPCParam->GetNRow(isec); + for (Int_t irow=0; irowWrite(); + + savedir->cd(); //switch back to the input file - fClustersIndex[fNsectors]--; // set end of the clusters buffer +} // end of function + +//_________________________________________________________________ +void AliTPC::Hits2ExactClustersSector(Int_t isec) +{ + //-------------------------------------------------------- + //calculate exact cross point of track and given pad row + //resulting values are expressed in "digit" coordinata + //-------------------------------------------------------- + + //----------------------------------------------------------------- + // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de + //----------------------------------------------------------------- + // + if (fClustersArray==0){ + return; + } + // + TParticle *particle; // pointer to a given particle + AliTPChit *tpcHit; // pointer to a sigle TPC hit + Int_t sector,nhits; + Int_t ipart; + const Int_t kcmaxhits=30000; + TVector * xxxx = new TVector(kcmaxhits*4); + TVector & xxx = *xxxx; + Int_t maxhits = kcmaxhits; + //construct array for each padrow + for (Int_t i=0; iGetNRow(isec);i++) + fClustersArray->CreateRow(isec,i); + //--------------------------------------------------------------- + // Get the access to the tracks + //--------------------------------------------------------------- + + TTree *tH = gAlice->TreeH(); + Stat_t ntracks = tH->GetEntries(); + Int_t npart = gAlice->GetNtrack(); + + //------------------------------------------------------------ + // Loop over tracks + //------------------------------------------------------------ + + for(Int_t track=0;trackGetEvent(track); + // + // Get number of the TPC hits and a pointer + // to the particles + // + nhits=fHits->GetEntriesFast(); + // + // Loop over hits + // + Int_t currentIndex=0; + Int_t lastrow=-1; //last writen row + for(Int_t hit=0;hitUncheckedAt(hit); + if (tpcHit==0) continue; + sector=tpcHit->fSector; // sector number + if(sector != isec) continue; + ipart=tpcHit->Track(); + if (ipartParticle(ipart); + + //find row number + + Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()}; + Int_t index[3]={1,isec,0}; + Int_t currentrow = fTPCParam->GetPadRow(x,index) ; + if (currentrow<0) continue; + if (lastrow<0) lastrow=currentrow; + if (currentrow==lastrow){ + if ( currentIndex>=maxhits){ + maxhits+=kcmaxhits; + xxx.ResizeTo(4*maxhits); + } + xxx(currentIndex*4)=x[0]; + xxx(currentIndex*4+1)=x[1]; + xxx(currentIndex*4+2)=x[2]; + xxx(currentIndex*4+3)=tpcHit->fQ; + currentIndex++; + } + else + if (currentIndex>2){ + Float_t sumx=0; + Float_t sumx2=0; + Float_t sumx3=0; + Float_t sumx4=0; + Float_t sumy=0; + Float_t sumxy=0; + Float_t sumx2y=0; + Float_t sumz=0; + Float_t sumxz=0; + Float_t sumx2z=0; + Float_t sumq=0; + for (Int_t index=0;indexGetNPads(isec,lastrow)-1)/2; + Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+ + sumx2*(sumx*sumx3-sumx2*sumx2); + + Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+ + sumx2*(sumxy*sumx3-sumx2y*sumx2); + Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+ + sumx2*(sumxz*sumx3-sumx2z*sumx2); + + Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+ + sumx2*(sumx*sumx2y-sumx2*sumxy); + Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+ + sumx2*(sumx*sumx2z-sumx2*sumxz); + + Float_t y=detay/det+centralPad; + Float_t z=detaz/det; + Float_t by=detby/det; //y angle + Float_t bz=detbz/det; //z angle + sumy/=Float_t(currentIndex); + sumz/=Float_t(currentIndex); + AliComplexCluster cl; + cl.fX=z; + cl.fY=y; + cl.fQ=sumq; + cl.fSigmaX2=bz; + cl.fSigmaY2=by; + cl.fTracks[0]=ipart; + + AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow)); + if (row!=0) row->InsertCluster(&cl); + currentIndex=0; + lastrow=currentrow; + } //end of calculating cluster for given row + + + + } // end of loop over hits + } // end of loop over tracks + //write padrows to tree + for (Int_t ii=0; iiGetNRow(isec);ii++) { + fClustersArray->StoreRow(isec,ii); + fClustersArray->ClearRow(isec,ii); + } + xxxx->Delete(); + } +//___________________________________________ +void AliTPC::SDigits2Digits() +{ + AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60"); + AliTPCPRF2D * prfinner = new AliTPCPRF2D; + AliTPCPRF2D * prfouter = new AliTPCPRF2D; + AliTPCRF1D * rf = new AliTPCRF1D(kTRUE); + + TDirectory *cwd = gDirectory; + rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.); + rf->SetOffset(3*param->GetZSigma()); + rf->Update(); + TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root"); + if (!f->IsOpen()) { + cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n"; + exit(3); + } + prfinner->Read("prf_07504_Gati_056068_d02"); + prfouter->Read("prf_10006_Gati_047051_d03"); + f->Close(); + cwd->cd(); + + param->SetInnerPRF(prfinner); + param->SetOuterPRF(prfouter); + param->SetTimeRF(rf); + + SetParam(param); + + cerr<<"Digitizing TPC...\n"; + + //setup TPCDigitsArray + AliTPCDigitsArray *arr = new AliTPCDigitsArray; + arr->SetClass("AliSimDigits"); + arr->Setup(param); + SetParam(param); + + arr->MakeTree(fDigitsFile); + + SetDigitsArray(arr); + + Hits2Digits(); + + // Hits2DigitsSector(1); + // Hits2DigitsSector(2); + // Hits2DigitsSector(3); + // Hits2DigitsSector(1+18); + // Hits2DigitsSector(2+18); + // Hits2DigitsSector(3+18); + + // Hits2DigitsSector(36+1); + // Hits2DigitsSector(36+2); + // Hits2DigitsSector(36+3); + // Hits2DigitsSector(36+1+18); + // Hits2DigitsSector(36+2+18); + // Hits2DigitsSector(36+3+18); + + //write results + char treeName[100]; + sprintf(treeName,"TreeD_%s",param->GetTitle()); + GetDigitsArray()->GetTree()->Write(treeName,TObject::kOverwrite); +} +//__________________________________________________________________ void AliTPC::Hits2Digits() { - //---------------------------------------------------- - // Loop over all sectors (72 sectors) - // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0 - // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0 - //---- - for(Int_t isec=1;isecGetNSector();isec++) Hits2DigitsSector(isec); + } +//__________________________________________________________________ +void AliTPC::Hits2SDigits() +{ + + //----------------------------------------------------------- + // summable digits - 16 bit "ADC", no noise, no saturation + //----------------------------------------------------------- + + //---------------------------------------------------- + // Loop over all sectors + //---------------------------------------------------- + + if(fTPCParam == 0){ + printf("AliTPCParam MUST be created firstly\n"); + return; + } + + fDigitsSwitch=1; + + for(Int_t isec=0;isecGetNSector();isec++) Hits2DigitsSector(isec); + +} //_____________________________________________________________________________ void AliTPC::Hits2DigitsSector(Int_t isec) @@ -1027,9 +1261,9 @@ void AliTPC::Hits2DigitsSector(Int_t isec) // Get the access to the track hits //------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - TTree *TH = gAlice->TreeH(); // pointer to the hits tree - Stat_t ntracks = TH->GetEntries(); + + TTree *tH = gAlice->TreeH(); // pointer to the hits tree + Stat_t ntracks = tH->GetEntries(); if( ntracks > 0){ @@ -1037,28 +1271,15 @@ void AliTPC::Hits2DigitsSector(Int_t isec) // Only if there are any tracks... //------------------------------------------- - - // TObjArrays for three neighouring pad-rows - - TObjArray **rowTriplet = new TObjArray* [3]; - - // TObjArray-s for each pad-row - TObjArray **row; - - for (Int_t trip=0;trip<3;trip++){ - rowTriplet[trip]=new TObjArray; - } - - - printf("*** Processing sector number %d ***\n",isec); + //printf("*** Processing sector number %d ***\n",isec); Int_t nrows =fTPCParam->GetNRow(isec); row= new TObjArray* [nrows]; - MakeSector(isec,nrows,TH,ntracks,row); + MakeSector(isec,nrows,tH,ntracks,row); //-------------------------------------------------------- // Digitize this sector, row by row @@ -1069,203 +1290,40 @@ void AliTPC::Hits2DigitsSector(Int_t isec) Int_t i; - for (i=0;iGetTree()==0) fDigitsArray->MakeTree(fDigitsFile); - // Triplets for i = 0 and i=1 are identical! - // The same for i = nrows-1 and nrows! - - if(i != 1 && i != nrows-1){ - MakeTriplet(i,rowTriplet,row); - } + for (i=0;iCreateRow(isec,i); - fDigParam->Fill(); + DigitizeRow(i,isec,row); - Int_t ndig=fDigParam->GetArray()->GetEntriesFast(); + fDigitsArray->StoreRow(isec,i); - printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig); + Int_t ndig = dig->GetDigitSize(); + + + //printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig); + + fDigitsArray->ClearRow(isec,i); - ResetDigits(); // reset digits for this row after storing them - + } // end of the sector digitization - - // delete the last triplet - for (i=0;i<3;i++) rowTriplet[i]->Delete(); - + for(i=0;iDelete(); + delete row[i]; + } + delete [] row; // delete the array of pointers to TObjArray-s } // ntracks >0 -} // end of Hits2Digits -//_____________________________________________________________________________ -void AliTPC::MakeTriplet(Int_t row, - TObjArray **rowTriplet, TObjArray **prow) -{ - //------------------------------------------------------------------ - // Makes the "triplet" of the neighbouring pad-row for the - // digitization including the cross-talk between the pad-rows - //------------------------------------------------------------------ - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t gasgain = fTPCParam->GetGasGain(); - Int_t nTracks[3]; - - Int_t nElements,nElectrons; - - TVector *pv; - TVector *tv; - - //------------------------------------------------------------------- - // pv is an "old" track, i.e. label + triplets of (x,y,z) - // for each electron - // - //------------------------------------------------------------------- - - - Int_t i1,i2; - Int_t nel,nt; - - if(row == 0 || row == 1){ - - // create entire triplet for the first AND the second row - - nTracks[0] = prow[0]->GetEntries(); - nTracks[1] = prow[1]->GetEntries(); - nTracks[2] = prow[2]->GetEntries(); - - for(i1=0;i1<3;i1++){ - nt = nTracks[i1]; // number of tracks for this row - - for(i2=0;i2At(i2); - TVector &v1 = *pv; - nElements = pv->GetNrows(); - nElectrons = (nElements-1)/3; - tv = new TVector(4*nElectrons+1); // create TVector for a modified track - TVector &v2 = *tv; - v2(0)=v1(0); //track label +} // end of Hits2DigitsSector - for(nel=0;nelRndm())); - v2(idx2+1)= v1(idx1+1); - v2(idx2+2)= v1(idx1+2); - v2(idx2+3)= v1(idx1+3); - v2(idx2+4)= (Float_t)aval; // in number of electrons! - } // end of loop over electrons - // - // Add this track to a row - // - - rowTriplet[i1]->Add(tv); - - - } // end of loop over tracks for this row - - prow[i1]->Delete(); // remove "old" tracks - delete prow[i1]; // delete TObjArray for this row - prow[i1]=0; // set pointer to NULL - - } // end of loop over row triplets - - - } - else{ - - rowTriplet[0]->Delete(); // remove old lower row - - nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row - nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row - nTracks[2]=prow[row+1]->GetEntries(); // next row - - - //------------------------------------------- - // shift new tracks downwards - //------------------------------------------- - - for(i1=0;i1At(i1); - rowTriplet[0]->Add(pv); - } - - rowTriplet[1]->Clear(); // leave tracks on the heap!!! - - for(i1=0;i1At(i1); - rowTriplet[1]->Add(pv); - } - - rowTriplet[2]->Clear(); // leave tracks on the heap!!! - //--------------------------------------------- - // Create new upper row - //--------------------------------------------- - - - - for(i1=0;i1At(i1); - TVector &v1 = *pv; - nElements = pv->GetNrows(); - nElectrons = (nElements-1)/3; - - tv = new TVector(4*nElectrons+1); // create TVector for a modified track - TVector &v2 = *tv; - v2(0)=v1(0); //track label - - for(nel=0;nelRndm())); - - v2(idx2+1)= v1(idx1+1); - v2(idx2+2)= v1(idx1+2); - v2(idx2+3)= v1(idx1+3); - v2(idx2+4)= (Float_t)aval; // in number of electrons! - } // end of loop over electrons - - rowTriplet[2]->Add(tv); - - } // end of loop over tracks - - prow[row+1]->Delete(); // delete tracks for this row - delete prow[row+1]; // delete TObjArray for this row - prow[row+1]=0; // set a pointer to NULL - - } - -} // end of MakeTriplet -//_____________________________________________________________________________ -void AliTPC::ExB(Float_t *xyz) -{ - //------------------------------------------------ - // ExB at the wires and wire number calulation - //------------------------------------------------ - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - - Float_t x1=xyz[0]; - fTPCParam->GetWire(x1); //calculate nearest wire position - Float_t dx=xyz[0]-x1; - xyz[1]+=dx*fTPCParam->GetOmegaTau(); - -} // end of ExB //_____________________________________________________________________________ -void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet) +void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) { //----------------------------------------------------------- // Single row digitization, coupling from the neighbouring @@ -1274,139 +1332,79 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet) //----------------------------------------------------------------- // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl + // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t chipgain= fTPCParam->GetChipGain(); Float_t zerosup = fTPCParam->GetZeroSup(); Int_t nrows =fTPCParam->GetNRow(isec); + fCurrentIndex[1]= isec; - Int_t nTracks[3]; - Int_t n_of_pads[3]; - Int_t IndexRange[4]; - Int_t i1; - Int_t iFlag; - - // - // iFlag = 0 -> inner row, iFlag = 1 -> the middle one, iFlag = 2 -> the outer one - // - - nTracks[0]=rowTriplet[0]->GetEntries(); // lower row - nTracks[1]=rowTriplet[1]->GetEntries(); // middle row - nTracks[2]=rowTriplet[2]->GetEntries(); // upper row - - if(irow == 0){ - iFlag=0; - n_of_pads[0]=fTPCParam->GetNPads(isec,0); - n_of_pads[1]=fTPCParam->GetNPads(isec,1); - } - else if(irow == nrows-1){ - iFlag=2; - n_of_pads[1]=fTPCParam->GetNPads(isec,irow-1); - n_of_pads[2]=fTPCParam->GetNPads(isec,irow); - } - else { - iFlag=1; - for(i1=0;i1<3;i1++){ - n_of_pads[i1]=fTPCParam->GetNPads(isec,irow-1+i1); - } - } - + Int_t nofPads = fTPCParam->GetNPads(isec,irow); + Int_t nofTbins = fTPCParam->GetMaxTBin(); + Int_t indexRange[4]; // // Integrated signal for this row // and a single track signal - // - - TMatrix *m1 = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // integrated - TMatrix *m2 = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // single - + // + TMatrix *m1 = new TMatrix(0,nofPads,0,nofTbins); // integrated + TMatrix *m2 = new TMatrix(0,nofPads,0,nofTbins); // single // - - TMatrix &Total = *m1; + TMatrix &total = *m1; // Array of pointers to the label-signal list - Int_t NofDigits = n_of_pads[iFlag]*MAXTPCTBK; // number of digits for this row - - Float_t **pList = new Float_t* [NofDigits]; + Int_t nofDigits = nofPads*nofTbins; // number of digits for this row + Float_t **pList = new Float_t* [nofDigits]; Int_t lp; - - for(lp=0;lpZero(); // clear single track signal matrix - - Float_t TrackLabel = - GetSignal(rowTriplet[iFlag],i1,n_of_pads[iFlag],m2,m1,IndexRange); - - GetList(TrackLabel,n_of_pads[iFlag],m2,IndexRange,pList); - - } - - // - // Cross talk from the neighbouring pad-rows - // - - TMatrix *m3 = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // cross-talk - - TMatrix &Cross = *m3; - - if(iFlag == 0){ - - // cross-talk from the outer row only (first pad row) - - GetCrossTalk(0,rowTriplet[1],nTracks[1],n_of_pads,m3); - - } - else if(iFlag == 2){ - - // cross-talk from the inner row only (last pad row) - - GetCrossTalk(2,rowTriplet[1],nTracks[1],n_of_pads,m3); - + Int_t i1; + for(lp=0;lpGetNCrossRows(),0); + Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1); + for (Int_t row= row1;row<=row2;row++){ + Int_t nTracks= rows[row]->GetEntries(); + for (i1=0;i1Zero(); // clear single track signal matrix + Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange); + GetList(trackLabel,nofPads,m2,indexRange,pList); + } + else GetSignal(rows[row],i1,0,m1,indexRange); + } } - else{ - - // cross-talk from both inner and outer rows + + Int_t tracks[3]; - GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner - GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer - } + AliDigits *dig = fDigitsArray->GetRow(isec,irow); + for(Int_t ip=0;ipGaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); - Float_t q = Total(ip,it); + q = (Int_t)q; - Int_t gi =(it-1)*n_of_pads[iFlag]+ip-1; // global index + if(q <=zerosup) continue; // do not fill zeros + if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation - q = gRandom->Gaus(q,fTPCParam->GetNoise()); // apply noise - q *= (q_el*1.e15); // convert to fC - q *= chipgain; // convert to mV - q *= (adc_sat/dyn_range); // convert to ADC counts + } - if(q adc_sat) q = adc_sat; // saturation + else { + q *= 16.; + q = (Int_t)q; + } // // "real" signal or electronic noise (list = -1)? @@ -1416,15 +1414,20 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet) tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1; } - digits[0]=isec; - digits[1]=irow+1; - digits[2]=ip; - digits[3]=it; - digits[4]= (Int_t)q; - - // Add this digit - - AddDigit(tracks,digits); +//Begin_Html +/* + + using of AliDigits object +*/ +//End_Html + dig->SetDigitFast((Short_t)q,it,ip); + if (fDigitsArray->IsSimulated()) + { + ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0); + ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1); + ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2); + } + } // end of loop over time buckets } // end of lop over pads @@ -1433,7 +1436,7 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet) // This row has been digitized, delete nonused stuff // - for(lp=0;lpGetParam()); - AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D()); - AliTPCRF1D * fRF = &(fDigParam->GetRF()); - //to make the code faster we put parameters to the stack - - Float_t zwidth = fTPCParam->GetZWidth(); - Float_t zwidthm1 =1./zwidth; - tv = (TVector*)p1->At(ntr); // pointer to a track TVector &v = *tv; - Float_t label = v(0); - - Int_t CentralPad = (np+1)/2; - Int_t PadNumber; - Int_t nElectrons = (tv->GetNrows()-1)/4; - Float_t range=((np-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); // pad range - range -= 0.5; // dead zone, 5mm from the edge, according to H.G. Fischer - - Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above - - - Float_t PadSignal[7]; // signal from a single electron - - TMatrix &signal = *m1; - TMatrix &total = *m2; - - - IndexRange[0]=9999; // min pad - IndexRange[1]=-1; // max pad - IndexRange[2]=9999; //min time - IndexRange[3]=-1; // max time - - // - // Loop over all electrons - // - - for(Int_t nel=0; nelGetPadPitchWidth()){ - PadNumber=CentralPad; - } - else if (absy < range){ - PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.); - PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad); - } - else continue; // electron out of pad-range , lost at the sector edge - - Float_t aval = (absyGetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2; - Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth(); - for (Int_t i=0;i<7;i++){ - PadSignal[i]=fPRF2D->GetPRF(dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval; - PadSignal[i] *= fTPCParam->GetPadCoupling(); - } + Int_t nElectrons = (tv->GetNrows()-1)/4; + indexRange[0]=9999; // min pad + indexRange[1]=-1; // max pad + indexRange[2]=9999; //min time + indexRange[3]=-1; // max time - Int_t LeftPad = TMath::Max(1,PadNumber-3); - Int_t RightPad = TMath::Min(np,PadNumber+3); + // Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above - Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector - Int_t pmax=RightPad-PadNumber+3; // upper index - - Float_t z_drift = (z_end-z)*zwidthm1; - Float_t z_offset = z_drift-(Int_t)z_drift; - //distance to the centre of nearest time bin (in time bin units) - Int_t FirstBucket = (Int_t)z_drift+1; - - - // loop over time bins (4 bins is enough - 3 sigma truncated Gaussian) - for (Int_t i2=0;i2<4;i2++){ - Int_t TrueTime = FirstBucket+i2; // current time bucket - Float_t dz = (Float_t(i2)+z_offset)*zwidth; - Float_t ampl = fRF->GetRF(dz); - if( (TrueTime>MAXTPCTBK) ) break; // beyond the time range - - IndexRange[2]=TMath::Min(IndexRange[2],TrueTime); // min time - IndexRange[3]=TMath::Max(IndexRange[3],TrueTime); // max time - - // loop over pads, from pmin to pmax - for(Int_t i3=pmin;i3<=pmax;i3++){ - Int_t TruePad = LeftPad+i3-pmin; - IndexRange[0]=TMath::Min(IndexRange[0],TruePad); // min pad - IndexRange[1]=TMath::Max(IndexRange[1],TruePad); // max pad - signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!! - total(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!! - } // end of pads loop - } // end of time loop + TMatrix &signal = *m1; + TMatrix &total = *m2; + // + // Loop over all electrons + // + for(Int_t nel=0; nelGetTotalNormFac(); + Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)}; + Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]); + + if (n>0) for (Int_t i =0; iGetResBin(i); + Int_t pad=index[1]+centralPad; //in digit coordinates central pad has coordinate 0 + if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) { + Int_t time=index[2]; + Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel + weight *= eltoadcfac; + + if (m1!=0) signal(pad,time)+=weight; + total(pad,time)+=weight; + indexRange[0]=TMath::Min(indexRange[0],pad); + indexRange[1]=TMath::Max(indexRange[1],pad); + indexRange[2]=TMath::Min(indexRange[2],time); + indexRange[3]=TMath::Max(indexRange[3],time); + } + } } // end of loop over electrons - + return label; // returns track label when finished } //_____________________________________________________________________________ -void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange, +void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange, Float_t **pList) { //---------------------------------------------------------------------- @@ -1577,43 +1530,45 @@ void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange, // lop over nonzero digits - for(Int_t it=IndexRange[2];it 0. // - if(signal(ip,it)>0.){ - - pList[GlobalIndex] = new Float_t [6]; + pList[globalIndex] = new Float_t [6]; // set list to -1 - *pList[GlobalIndex] = -1.; - *(pList[GlobalIndex]+1) = -1.; - *(pList[GlobalIndex]+2) = -1.; - *(pList[GlobalIndex]+3) = -1.; - *(pList[GlobalIndex]+4) = -1.; - *(pList[GlobalIndex]+5) = -1.; + *pList[globalIndex] = -1.; + *(pList[globalIndex]+1) = -1.; + *(pList[globalIndex]+2) = -1.; + *(pList[globalIndex]+3) = -1.; + *(pList[globalIndex]+4) = -1.; + *(pList[globalIndex]+5) = -1.; - *pList[GlobalIndex] = label; - *(pList[GlobalIndex]+3) = signal(ip,it);} + *pList[globalIndex] = label; + *(pList[globalIndex]+3) = signal(ip,it); } else{ // check the signal magnitude - Float_t highest = *(pList[GlobalIndex]+3); - Float_t middle = *(pList[GlobalIndex]+4); - Float_t lowest = *(pList[GlobalIndex]+5); + Float_t highest = *(pList[globalIndex]+3); + Float_t middle = *(pList[globalIndex]+4); + Float_t lowest = *(pList[globalIndex]+5); // // compare the new signal with already existing list @@ -1624,24 +1579,24 @@ void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange, // if (signal(ip,it)>highest){ - *(pList[GlobalIndex]+5) = middle; - *(pList[GlobalIndex]+4) = highest; - *(pList[GlobalIndex]+3) = signal(ip,it); + *(pList[globalIndex]+5) = middle; + *(pList[globalIndex]+4) = highest; + *(pList[globalIndex]+3) = signal(ip,it); - *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1); - *(pList[GlobalIndex]+1) = *pList[GlobalIndex]; - *pList[GlobalIndex] = label; + *(pList[globalIndex]+2) = *(pList[globalIndex]+1); + *(pList[globalIndex]+1) = *pList[globalIndex]; + *pList[globalIndex] = label; } else if (signal(ip,it)>middle){ - *(pList[GlobalIndex]+5) = middle; - *(pList[GlobalIndex]+4) = signal(ip,it); + *(pList[globalIndex]+5) = middle; + *(pList[globalIndex]+4) = signal(ip,it); - *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1); - *(pList[GlobalIndex]+1) = label; + *(pList[globalIndex]+2) = *(pList[globalIndex]+1); + *(pList[globalIndex]+1) = label; } else{ - *(pList[GlobalIndex]+5) = signal(ip,it); - *(pList[GlobalIndex]+2) = label; + *(pList[globalIndex]+5) = signal(ip,it); + *(pList[globalIndex]+2) = label; } } @@ -1650,7 +1605,6 @@ void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange, - }//end of GetList //___________________________________________________________________ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, @@ -1667,11 +1621,16 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); + Float_t gasgain = fTPCParam->GetGasGain(); Int_t i; - Float_t xyz[3]; + Float_t xyz[4]; AliTPChit *tpcHit; // pointer to a sigle TPC hit + //MI change + TBranch * branch=0; + if (fHitType&2) branch = TH->GetBranch("TPC2"); + else branch = TH->GetBranch("TPC"); + //---------------------------------------------- // Create TObjArray-s, one for each row, @@ -1679,11 +1638,15 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, // of electrons, one TVector per each track. //---------------------------------------------- + Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row + TVector **tracks = new TVector* [nrows]; //pointers to the track vectors for(i=0; iGetEvent(track); // get next track - Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track + if (fHitType&2) { + isInSector=kFALSE; + TBranch * br = TH->GetBranch("fTrackHitsInfo"); + br->GetEvent(track); + AliObjectArray * ar = fTrackHits->fTrackHitsInfo; + for (UInt_t j=0;jGetSize();j++){ + if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE; + } + } + if (!isInSector) continue; + //MI change + branch->GetEntry(track); // get next track + + //M.I. changes - if(nhits == 0) continue; // no hits in the TPC for this track + tpcHit = (AliTPChit*)FirstHit(-1); //-------------------------------------------------------------- // Loop over hits //-------------------------------------------------------------- - for(Int_t hit=0;hitUncheckedAt(hit); // get a pointer to a hit + while(tpcHit){ Int_t sector=tpcHit->fSector; // sector number - if(sector != isec) continue; //terminate iteration + // if(sector != isec) continue; + if(sector != isec){ + tpcHit = (AliTPChit*) NextHit(); + continue; + } + + currentTrack = tpcHit->Track(); // track number + - currentTrack = tpcHit->fTrack; // track number if(currentTrack != previousTrack){ // store already filled fTrack for(i=0;i0){ - TVector &v = *tr[i]; + if(nofElectrons[i]>0){ + TVector &v = *tracks[i]; v(0) = previousTrack; - tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary - row[i]->Add(tr[i]); + tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary + row[i]->Add(tracks[i]); } else{ - delete tr[i]; // delete empty TVector - tr[i]=0; + delete tracks[i]; // delete empty TVector + tracks[i]=0; } } - n_of_electrons[i]=0; - tr[i] = new TVector(361); // TVectors for the next fTrack + nofElectrons[i]=0; + tracks[i] = new TVector(481); // TVectors for the next fTrack } // end of loop over rows previousTrack=currentTrack; // update track label } - Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons) + Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons) //--------------------------------------------------- // Calculate the electron attachment probability //--------------------------------------------------- - Float_t time = 1.e6*(z_end-TMath::Abs(tpcHit->fZ))/fTPCParam->GetDriftV(); + + Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z())) + /fTPCParam->GetDriftV(); // in microseconds! - Float_t AttProb = fTPCParam->GetAttCoef()* + Float_t attProb = fTPCParam->GetAttCoef()* fTPCParam->GetOxyCont()*time; // fraction! //----------------------------------------------- // Loop over electrons //----------------------------------------------- - - for(Int_t nel=0;nelRndm(0)) < AttProb) continue; // electron lost! - xyz[0]=tpcHit->fX; - xyz[1]=tpcHit->fY; - xyz[2]=tpcHit->fZ; - ElDiff(xyz); // Appply the diffusion - Int_t row_number; - fTPCParam->XYZtoCRXYZ(xyz,isec,row_number,3); - - //transform position to local coordinates - //option 3 means that we calculate x position relative to - //nearest pad row - - if ((row_number<0)||row_number>=fTPCParam->GetNRow(isec)) continue; - ExB(xyz); // ExB effect at the sense wires - n_of_electrons[row_number]++; + if((gRandom->Rndm(0)) < attProb) continue; // electron lost! + xyz[0]=tpcHit->X(); + xyz[1]=tpcHit->Y(); + xyz[2]=tpcHit->Z(); + // + // protection for the nonphysical avalanche size (10**6 maximum) + // + Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22); + xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); + index[0]=1; + + TransportElectron(xyz,index); //MI change -august + Int_t rowNumber; + fTPCParam->GetPadRow(xyz,index); //MI change august + rowNumber = index[2]; + //transform position to local digit coordinates + //relative to nearest pad row + if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue; + nofElectrons[rowNumber]++; //---------------------------------- // Expand vector if necessary //---------------------------------- - if(n_of_electrons[row_number]>120){ - Int_t range = tr[row_number]->GetNrows(); - if(n_of_electrons[row_number] > (range-1)/3){ - tr[row_number]->ResizeTo(range+150); // Add 50 electrons + if(nofElectrons[rowNumber]>120){ + Int_t range = tracks[rowNumber]->GetNrows(); + if((nofElectrons[rowNumber])>(range-1)/4){ + + tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons } } - TVector &v = *tr[row_number]; - Int_t idx = 3*n_of_electrons[row_number]-2; + TVector &v = *tracks[rowNumber]; + Int_t idx = 4*nofElectrons[rowNumber]-3; - v(idx)= xyz[0]; // X - v(idx+1)=xyz[1]; // Y (along the pad-row) - v(idx+2)=xyz[2]; // Z - + v(idx)= xyz[0]; // X - pad row coordinate + v(idx+1)=xyz[1]; // Y - pad coordinate (along the pad-row) + v(idx+2)=xyz[2]; // Z - time bin coordinate + v(idx+3)=xyz[3]; // avalanche size } // end of loop over electrons + + tpcHit = (AliTPChit*)NextHit(); } // end of loop over hits } // end of loop over tracks @@ -1798,171 +1788,23 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, // for(i=0;i0){ - TVector &v = *tr[i]; + if(nofElectrons[i]>0){ + TVector &v = *tracks[i]; v(0) = previousTrack; - tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary - row[i]->Add(tr[i]); + tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary + row[i]->Add(tracks[i]); } else{ - delete tr[i]; - tr[i]=0; + delete tracks[i]; + tracks[i]=0; } } - delete [] tr; - delete [] n_of_electrons; - -} // end of MakeSector -//_____________________________________________________________________________ -void AliTPC::GetCrossTalk (Int_t iFlag,TObjArray *p,Int_t ntracks,Int_t *npads, - TMatrix *m) -{ - - //------------------------------------------------------------- - // Calculates the cross-talk from one row (inner or outer) - //------------------------------------------------------------- - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - - // - // iFlag=2 & 3 --> cross-talk from the inner row - // iFlag=0 & 4 --> cross-talk from the outer row - // - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D()); - AliTPCRF1D * fRF = &(fDigParam->GetRF()); + delete [] tracks; + delete [] nofElectrons; - //to make code faster - - Float_t zwidth = fTPCParam->GetZWidth(); - Float_t zwidthm1 =1/fTPCParam->GetZWidth(); - - Int_t PadNumber; - Float_t xwire; - - Int_t nPadsSignal; // for this pads the signal is calculated - Float_t range; // sense wire range - Int_t nPadsDiff; - - Float_t IneffFactor=0.5; // gain inefficiency close to the edges - - if(iFlag == 0){ - // 1-->0 - nPadsSignal = npads[1]; - range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); - nPadsDiff = (npads[1]-npads[0])/2; - } - else if (iFlag == 2){ - // 1-->2 - nPadsSignal = npads[2]; - range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); - nPadsDiff = 0; - } - else if (iFlag == 3){ - // 0-->1 - nPadsSignal = npads[1]; - range = ((npads[0]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); - nPadsDiff = 0; - } - else{ - // 2-->1 - nPadsSignal = npads[2]; - range = ((npads[2]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); - nPadsDiff = (npads[2]-npads[1])/2; - } - - range-=0.5; // dead zone close to the edges - - TVector *tv; - TMatrix &signal = *m; - - Int_t CentralPad = (nPadsSignal+1)/2; - Float_t PadSignal[7]; // signal from a single electron - // Loop over tracks - for(Int_t nt=0;ntAt(nt); // pointer to a track - TVector &v = *tv; - Int_t nElectrons = (tv->GetNrows()-1)/4; - // Loop over electrons - for(Int_t nel=0; nelGetPadPitchLength(); - if (iFlag==2) xwire-=fTPCParam->GetPadPitchLength(); - if (iFlag==3) xwire-=fTPCParam->GetPadPitchLength(); - if (iFlag==4) xwire+=fTPCParam->GetPadPitchLength(); - - // electron acceptance for the cross-talk (only the closest wire) - - Float_t dxMax = fTPCParam->GetPadPitchLength()*0.5+fTPCParam->GetWWPitch(); - if(TMath::Abs(xwire)>dxMax) continue; - - Float_t y = v(idx+2); - Float_t z = v(idx+3); - Float_t absy=TMath::Abs(y); - - if(absy < 0.5*fTPCParam->GetPadPitchWidth()){ - PadNumber=CentralPad; - } - else if (absy < range){ - PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.); - PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad); - } - else continue; // electron out of sense wire range, lost at the sector edge - - Float_t aval = (absyGetPadPitchWidth(); - - for (Int_t i=0;i<7;i++){ - PadSignal[i]=fPRF2D->GetPRF(dist+(3-i)*fTPCParam->GetPadPitchWidth(),xwire)*aval; - - PadSignal[i] *= fTPCParam->GetPadCoupling(); - } - // real pad range - - Int_t LeftPad = TMath::Max(1,PadNumber-3); - Int_t RightPad = TMath::Min(nPadsSignal,PadNumber+3); - - Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector - Int_t pmax=RightPad-PadNumber+3; // upper index - - - Float_t z_drift = (z_end-z)*zwidthm1; - Float_t z_offset = z_drift-(Int_t)z_drift; - //distance to the centre of nearest time bin (in time bin units) - Int_t FirstBucket = (Int_t)z_drift+1; - // MI check it --time offset - for (Int_t i2=0;i2<4;i2++){ - Int_t TrueTime = FirstBucket+i2; // current time bucket - Float_t dz = (Float_t(i2)+z_offset)*zwidth; - Float_t ampl = fRF->GetRF(dz); - if((TrueTime>MAXTPCTBK)) break; // beyond the time range - - - // loop over pads, from pmin to pmax - - for(Int_t i3=pmin;i3 nPadsSignal-nPadsDiff) continue; - - TruePad -= nPadsDiff; - signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge! - - } // end of loop over pads - } // end of loop over time bins - - } // end of loop over electrons - - } // end of loop over tracks - -} // end of GetCrossTalk +} // end of MakeSector //_____________________________________________________________________________ @@ -1973,18 +1815,20 @@ void AliTPC::Init() // Int_t i; // - printf("\n"); - for(i=0;i<35;i++) printf("*"); - printf(" TPC_INIT "); - for(i=0;i<35;i++) printf("*"); - printf("\n"); - // - for(i=0;i<80;i++) printf("*"); - printf("\n"); + if(fDebug) { + printf("\n%s: ",ClassName()); + for(i=0;i<35;i++) printf("*"); + printf(" TPC_INIT "); + for(i=0;i<35;i++) printf("*"); + printf("\n%s: ",ClassName()); + // + for(i=0;i<80;i++) printf("*"); + printf("\n"); + } } //_____________________________________________________________________________ -void AliTPC::MakeBranch(Option_t* option) +void AliTPC::MakeBranch(Option_t* option, const char *file) { // // Create Tree branches for the TPC. @@ -1993,21 +1837,16 @@ void AliTPC::MakeBranch(Option_t* option) char branchname[10]; sprintf(branchname,"%s",GetName()); - AliDetector::MakeBranch(option); + AliDetector::MakeBranch(option,file); - char *D = strstr(option,"D"); + const char *d = strstr(option,"D"); - if (fDigits && gAlice->TreeD() && D) { - gAlice->TreeD()->Branch(branchname,&fDigits, buffersize); - printf("Making Branch %s for digits\n",branchname); + if (fDigits && gAlice->TreeD() && d) { + MakeBranchInTree(gAlice->TreeD(), + branchname, &fDigits, buffersize, file); } - char *R = strstr(option,"R"); - - if (fClusters && gAlice->TreeR() && R) { - gAlice->TreeR()->Branch(branchname,&fClusters, buffersize); - printf("Making Branch %s for Clusters\n",branchname); - } + if (fHitType&2) MakeBranch2(option,file); // MI change 14.09.2000 } //_____________________________________________________________________________ @@ -2015,13 +1854,9 @@ void AliTPC::ResetDigits() { // // Reset number of digits and the digits array for this detector - // reset clusters // fNdigits = 0; - // if (fDigits) fDigits->Clear(); - if (fDigParam->GetArray()!=0) fDigParam->GetArray()->Clear(); - fNclusters = 0; - if (fClusters) fClusters->Clear(); + if (fDigits) fDigits->Clear(); } //_____________________________________________________________________________ @@ -2114,82 +1949,70 @@ void AliTPC::SetSens(Int_t sens) fSens = sens; } -//_____________________________________________________________________________ -void AliTPC::Streamer(TBuffer &R__b) -{ - // - // Stream an object of class AliTPC. - // - if (R__b.IsReading()) { - Version_t R__v = R__b.ReadVersion(); if (R__v) { } - AliDetector::Streamer(R__b); - if (R__v < 2) return; - R__b >> fNsectors; - R__b >> fNclusters; - R__b >> fNtracks; - fClustersIndex = new Int_t[fNsectors+1]; - fDigitsIndex = new Int_t[fNsectors+1]; - } else { - R__b.WriteVersion(AliTPC::IsA()); - AliDetector::Streamer(R__b); - R__b << fNsectors; - R__b << fNclusters; - R__b << fNtracks; - } -} - -ClassImp(AliTPCcluster) -//_____________________________________________________________________________ -AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab) +void AliTPC::SetSide(Float_t side=0.) { - // - // Creates a simulated cluster for the TPC - // - fTracks[0] = lab[0]; - fTracks[1] = lab[1]; - fTracks[2] = lab[2]; - fSector = lab[3]; - fPadRow = lab[4]; - fY = hits[0]; - fZ = hits[1]; - fQ = hits[2]; - fSigmaY2 = hits[3]; - fSigmaZ2 = hits[4]; + // choice of the TPC side + + fSide = side; + } +//____________________________________________________________________________ +void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1, + Float_t p2,Float_t p3) +{ + + // gax mixture definition + + fNoComp = nc; + + fMixtComp[0]=c1; + fMixtComp[1]=c2; + fMixtComp[2]=c3; + + fMixtProp[0]=p1; + fMixtProp[1]=p2; + fMixtProp[2]=p3; + +} //_____________________________________________________________________________ -void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const + +void AliTPC::TransportElectron(Float_t *xyz, Int_t *index) { // - // Transformation from local to global coordinate system + // electron transport taking into account: + // 1. diffusion, + // 2.ExB at the wires + // 3. nonisochronity // - x[0]=par->GetPadRowRadii(fSector,fPadRow-1); - x[1]=fY; - x[2]=fZ; - par->CRXYZtoXYZ(x,fSector,fPadRow-1,1); -} - -//_____________________________________________________________________________ -Int_t AliTPCcluster::Compare(TObject * o) -{ + // xyz and index must be already transformed to system 1 // - // compare two clusters according y coordinata - AliTPCcluster *cl= (AliTPCcluster *)o; - if (fYfY) return -1; - if (fY==cl->fY) return 0; - return 1; -} -Bool_t AliTPCcluster::IsSortable() const -{ - // - //make AliTPCcluster sortabale - return kTRUE; + fTPCParam->Transform1to2(xyz,index); + + //add diffusion + Float_t driftl=xyz[2]; + if(driftl<0.01) driftl=0.01; + driftl=TMath::Sqrt(driftl); + Float_t sigT = driftl*(fTPCParam->GetDiffT()); + Float_t sigL = driftl*(fTPCParam->GetDiffL()); + xyz[0]=gRandom->Gaus(xyz[0],sigT); + xyz[1]=gRandom->Gaus(xyz[1],sigT); + xyz[2]=gRandom->Gaus(xyz[2],sigL); + + // ExB + + if (fTPCParam->GetMWPCReadout()==kTRUE){ + Float_t x1=xyz[0]; + fTPCParam->Transform2to2NearestWire(xyz,index); + Float_t dx=xyz[0]-x1; + xyz[1]+=dx*(fTPCParam->GetOmegaTau()); + } + //add nonisochronity (not implemented yet) + } - - - + ClassImp(AliTPCdigit) //_____________________________________________________________________________ @@ -2224,389 +2047,457 @@ AliHit(shunt,track) fQ = hits[3]; } - -ClassImp(AliTPCtrack) - -//_____________________________________________________________________________ -AliTPCtrack::AliTPCtrack(Float_t *hits) + +//________________________________________________________________________ +// Additional code because of the AliTPCTrackHits + +void AliTPC::MakeBranch2(Option_t *option,const char *file) { // - // Default creator for a TPC reconstructed track object - // - ref=hits[0]; // This is dummy code ! + // Create a new branch in the current Root Tree + // The branch of fHits is automatically split + // MI change 14.09.2000 + if (fHitType&2==0) return; + char branchname[10]; + sprintf(branchname,"%s2",GetName()); + // + // Get the pointer to the header + const char *cH = strstr(option,"H"); + // + if (fTrackHits && gAlice->TreeH() && cH) { + AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHits, + gAlice->TreeH(),fBufferSize,1); + gAlice->TreeH()->GetListOfBranches()->Add(branch); + if (GetDebug()>1) + printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname); + const char folder [] = "RunMC/Event/Data"; + if (GetDebug()) + printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder); + Publish(folder,&fTrackHits,branchname); + if (file) { + TBranch *b = gAlice->TreeH()->GetBranch(branchname); + TDirectory *wd = gDirectory; + b->SetFile(file); + TIter next( b->GetListOfBranches()); + while ((b=(TBranch*)next())) { + b->SetFile(file); + } + wd->cd(); + if (GetDebug()>1) + cout << "Diverting branch " << branchname << " to file " << file << endl; + } + } } -AliTPCtrack::AliTPCtrack(const AliTPCcluster& c,const TVector& xx, - const TMatrix& CC, const AliTPCParam *p): - x(xx),C(CC),clusters(MAX_CLUSTER) +void AliTPC::SetTreeAddress() { - // - // Standard creator for a TPC reconstructed track object - // - chi2=0.; - int sec=c.fSector-1, row=c.fPadRow-1; - ref=p->GetPadRowRadii(sec+1,row); + if (fHitType&1) AliDetector::SetTreeAddress(); + if (fHitType&2) SetTreeAddress2(); +} - if (sec<24) { - fAlpha=(sec%12)*alpha_low; - } else { - fAlpha=(sec%24)*alpha_up; +void AliTPC::SetTreeAddress2() +{ + // + // Set branch address for the TrackHits Tree + // + TBranch *branch; + char branchname[20]; + sprintf(branchname,"%s2",GetName()); + // + // Branch address for hit tree + TTree *treeH = gAlice->TreeH(); + if (treeH) { + branch = treeH->GetBranch(branchname); + if (branch) branch->SetAddress(&fTrackHits); } - clusters.AddLast((AliTPCcluster*)(&c)); } -//_____________________________________________________________________________ -AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C), - clusters(t.clusters.GetEntriesFast()) +void AliTPC::FinishPrimary() { + if (fTrackHits) fTrackHits->FlushHitStack(); +} + + +void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits) +{ // - // Copy creator for a TPC reconstructed track - // - ref=t.ref; - chi2=t.chi2; - fAlpha=t.fAlpha; - int n=t.clusters.GetEntriesFast(); - for (int i=0; iGetPrimary(track); + gAlice->Particle(primary)->SetBit(kKeepBit); + rtrack=primary; + } else { + rtrack=track; + gAlice->FlagTrack(track); + } + //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1); + //if (hit->fTrack!=rtrack) + // cout<<"bad track number\n"; + if (fTrackHits) + fTrackHits->AddHitKartez(vol[0],rtrack, hits[0], + hits[1],hits[2],(Int_t)hits[3]); } -//_____________________________________________________________________________ -Double_t AliTPCtrack::GetY(Double_t xk) const +void AliTPC::ResetHits() { - // - // - // - Double_t c2=x(2)*xk - x(3); - if (TMath::Abs(c2) >= 0.999) { - if (*this>10) cerr<<*this<<" AliTPCtrack warning: No y for this x !\n"; - return 0.; - } - Double_t c1=x(2)*ref - x(3); - Double_t r1=sqrt(1.-c1*c1), r2=sqrt(1.-c2*c2); - Double_t dx=xk-ref; - return x(0) + dx*(c1+c2)/(r1+r2); + if (fHitType&1) AliDetector::ResetHits(); + if (fHitType&2) ResetHits2(); } -//_____________________________________________________________________________ -int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) +void AliTPC::ResetHits2() { // - // Propagate a TPC reconstructed track - // - if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) { - if (*this>3) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n"; - return 0; - } + //reset hits + if (fTrackHits) fTrackHits->Clear(); +} - Double_t x1=ref, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1); - Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1); - Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2); - - x(0) += dx*(c1+c2)/(r1+r2); - x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4); - - TMatrix F(5,5); F.UnitMatrix(); - Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2; - F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr); - F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr); - Double_t cr=c1*r2+c2*r1; - F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr); - F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr); - F(1,4)= dx*cc/cr; - TMatrix tmp(F,TMatrix::kMult,C); - C.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); - - ref=x2; - - //Multiple scattering****************** - Double_t ey=x(2)*ref - x(3); - Double_t ex=sqrt(1-ey*ey); - Double_t ez=x(4); - TMatrix Q(5,5); Q=0.; - Q(2,2)=ez*ez+ey*ey; Q(2,3)=-ex*ey; Q(2,4)=-ex*ez; - Q(3,2)=Q(2,3); Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez; - Q(4,2)=Q(2,4); Q(4,3)= Q(3,4); Q(4,4)=1.; - - F=0; - F(2,2)=-x(2)*ex; F(2,3)=-x(2)*ey; - F(3,2)=-ex*(x(2)*ref-ey); F(3,3)=-(1.+ x(2)*ref*ey - ey*ey); - F(4,2)=-ez*ex; F(4,3)=-ez*ey; F(4,4)=1.; - - tmp.Mult(F,Q); - Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); - - Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4)); - Double_t beta2=p2/(p2 + pm*pm); - Double_t d=sqrt((x1-ref)*(x1-ref)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1))); - d*=2.; - Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho; - Q*=theta2; - C+=Q; - - //Energy losses************************ - Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho; - if (x1 < x2) dE=-dE; - x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE); - //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE); - - x1=ref; x2=xk; y1=x(0); z1=x(1); - c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1); - c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2); - - x(0) += dx*(c1+c2)/(r1+r2); - x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4); - - F.UnitMatrix(); - rr=r1+r2; cc=c1+c2; xx=x1+x2; - F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr); - F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr); - cr=c1*r2+c2*r1; - F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr); - F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr); - F(1,4)= dx*cc/cr; - tmp.Mult(F,C); - C.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); - - ref=x2; - - return 1; +AliHit* AliTPC::FirstHit(Int_t track) +{ + if (fHitType&2) return FirstHit2(track); + return AliDetector::FirstHit(track); } - -//_____________________________________________________________________________ -void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) +AliHit* AliTPC::NextHit() { - // - // Propagate a reconstructed track from the vertex - // - Double_t c=x(2)*ref - x(3); - Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c)); - Double_t snf=tgf/sqrt(1.+ tgf*tgf); - Double_t xv=(x(3)+snf)/x(2); - PropagateTo(xv,x0,rho,pm); + if (fHitType&2) return NextHit2(); + return AliDetector::NextHit(); } -//_____________________________________________________________________________ -void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq) +AliHit* AliTPC::FirstHit2(Int_t track) { // - // Update statistics for a reconstructed TPC track + // Initialise the hit iterator + // Return the address of the first hit for track + // If track>=0 the track is read from disk + // while if track<0 the first hit of the current + // track is returned + // + if(track>=0) { + gAlice->ResetHits(); + gAlice->TreeH()->GetEvent(track); + } // - TMatrix H(2,5); H.UnitMatrix(); - TMatrix Ht(TMatrix::kTransposed,H); - TVector m(2); m(0)=c->fY; m(1)=c->fZ; - TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2; - - TMatrix tmp(H,TMatrix::kMult,C); - TMatrix R(tmp,TMatrix::kMult,Ht); R+=V; - - Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0); - R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1); - R(1,0)*=-1; R(0,1)=R(1,0); - R*=1./det; - - //R.Invert(); - - TMatrix K(C,TMatrix::kMult,Ht); K*=R; - - TVector savex=x; - x*=H; x-=m; x*=-1; x*=K; x+=savex; - if (TMath::Abs(x(2)*ref-x(3)) >= 0.999) { - if (*this>3) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n"; - x=savex; - return; + if (fTrackHits) { + fTrackHits->First(); + return fTrackHits->GetHit(); } - - TMatrix saveC=C; - C.Mult(K,tmp); C-=saveC; C*=-1; - - clusters.AddLast((AliTPCcluster*)c); - chi2 += chisq; + else return 0; } -//_____________________________________________________________________________ -int AliTPCtrack::Rotate(Double_t alpha) +AliHit* AliTPC::NextHit2() { // - // Rotate a reconstructed TPC track - // - fAlpha += alpha; - - Double_t x1=ref, y1=x(0); - Double_t ca=cos(alpha), sa=sin(alpha); - Double_t r1=x(2)*ref - x(3); - - ref = x1*ca + y1*sa; - x(0)=-x1*sa + y1*ca; - x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa; - - Double_t r2=x(2)*ref - x(3); - if (TMath::Abs(r2) >= 0.999) { - if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n"; - return 0; + //Return the next hit for the current track + + if (fTrackHits) { + fTrackHits->Next(); + return fTrackHits->GetHit(); } - - Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2); - if ((x(0)-y0)*x(2) >= 0.) { - if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n"; + else return 0; - } - - TMatrix F(5,5); F.UnitMatrix(); - F(0,0)=ca; - F(3,0)=x(2)*sa; - F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa; - F(3,3)= ca + sa*r1/sqrt(1.- r1*r1); - TMatrix tmp(F,TMatrix::kMult,C); - // Double_t dy2=C(0,0); - C.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); - // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1); - // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1); - - return 1; } -//_____________________________________________________________________________ -void AliTPCtrack::UseClusters() const +void AliTPC::LoadPoints(Int_t) { // - // - // - int num_of_clusters=clusters.GetEntriesFast(); - for (int i=0; iUse(); - } + Int_t a = 0; + /* if(fHitType==1) return AliDetector::LoadPoints(a); + LoadPoints2(a); + */ + if(fHitType==1) AliDetector::LoadPoints(a); + else LoadPoints2(a); + + // LoadPoints3(a); + } -//_____________________________________________________________________________ -Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const + +void AliTPC::RemapTrackHitIDs(Int_t *map) { - // - // Calculate chi2 for a reconstructed TPC track - // - TMatrix H(2,5); H.UnitMatrix(); - TVector m(2); m(0)=c->fY; m(1)=c->fZ; - TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2; - TVector res=x; res*=H; res-=m; //res*=-1; - TMatrix tmp(H,TMatrix::kMult,C); - TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V; - - Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0); - if (TMath::Abs(det) < 1.e-10) { - if (*this>3) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n"; - return 1e10; + if (!fTrackHits) return; + AliObjectArray * arr = fTrackHits->fTrackHitsInfo; + for (UInt_t i=0;iGetSize();i++){ + AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i)); + info->fTrackID = map[info->fTrackID]; } - R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1); - R(1,0)*=-1; R(0,1)=R(1,0); - R*=1./det; - //R.Invert(); - - TVector r=res; - res*=R; - return r*res; } + //_____________________________________________________________________________ -int AliTPCtrack::GetLab() const +void AliTPC::LoadPoints2(Int_t) { // + // Store x, y, z of all hits in memory // + if (fTrackHits == 0) return; // - int lab=123456789; - struct { - int lab; - int max; - } s[MAX_CLUSTER]={{0,0}}; - - int i; - int num_of_clusters=clusters.GetEntriesFast(); - for (i=0; ifTracks[0]); - int j; - for (j=0; jmax) {max=s[i].max; lab=s[i].lab;} - - for (i=0; ifTracks[1]) == lab || - TMath::Abs(c->fTracks[2]) == lab ) max++; + Int_t nhits = fTrackHits->GetEntriesFast(); + if (nhits == 0) return; + Int_t tracks = gAlice->GetNtrack(); + if (fPoints == 0) fPoints = new TObjArray(tracks); + AliHit *ahit; + // + Int_t *ntrk=new Int_t[tracks]; + Int_t *limi=new Int_t[tracks]; + Float_t **coor=new Float_t*[tracks]; + for(Int_t i=0;i 0.10) return -lab; - - if (num_of_clusters < 6) return lab; - - max=0; - for (i=1; i<=6; i++) { - AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(num_of_clusters-i); - if (lab == TMath::Abs(c->fTracks[0]) || - lab == TMath::Abs(c->fTracks[1]) || - lab == TMath::Abs(c->fTracks[2])) max++; + // + AliPoints *points = 0; + Float_t *fp=0; + Int_t trk; + Int_t chunk=nhits/4+1; + // + // Loop over all the hits and store their position + // + ahit = FirstHit2(-1); + //for (Int_t hit=0;hitUncheckedAt(hit); + trk=ahit->GetTrack(); + if(ntrk[trk]==limi[trk]) { + // + // Initialise a new track + fp=new Float_t[3*(limi[trk]+chunk)]; + if(coor[trk]) { + memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]); + delete [] coor[trk]; + } + limi[trk]+=chunk; + coor[trk] = fp; + } else { + fp = coor[trk]; + } + fp[3*ntrk[trk] ] = ahit->X(); + fp[3*ntrk[trk]+1] = ahit->Y(); + fp[3*ntrk[trk]+2] = ahit->Z(); + ntrk[trk]++; + ahit = NextHit2(); } - if (max<3) return -lab; - - return lab; -} - -//_____________________________________________________________________________ -void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const -{ // - // Get reconstructed TPC track momentum - // - Double_t pt=0.3*FIELD/TMath::Abs(x(2))/100; // GeV/c - Double_t r=x(2)*ref-x(3); - Double_t y0=x(0) + sqrt(1.- r*r)/x(2); - px=-pt*(x(0)-y0)*x(2); //cos(phi); - py=-pt*(x(3)-ref*x(2)); //sin(phi); - pz=pt*x(4); - Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha); - py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha); - px=tmp; + for(trk=0; trkSetMarkerColor(GetMarkerColor()); + points->SetMarkerSize(GetMarkerSize()); + points->SetDetector(this); + points->SetParticle(trk); + points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()); + fPoints->AddAt(points,trk); + delete [] coor[trk]; + coor[trk]=0; + } + } + delete [] coor; + delete [] ntrk; + delete [] limi; } -//_____________________________________________________________________________ -// -// Classes for internal tracking use -// //_____________________________________________________________________________ -void AliTPCRow::InsertCluster(const AliTPCcluster* c) +void AliTPC::LoadPoints3(Int_t) { // - // Insert a cluster in the list + // Store x, y, z of all hits in memory + // - only intersection point with pad row + if (fTrackHits == 0) return; + // + Int_t nhits = fTrackHits->GetEntriesFast(); + if (nhits == 0) return; + Int_t tracks = gAlice->GetNtrack(); + if (fPoints == 0) fPoints = new TObjArray(2*tracks); + fPoints->Expand(2*tracks); + AliHit *ahit; + // + Int_t *ntrk=new Int_t[tracks]; + Int_t *limi=new Int_t[tracks]; + Float_t **coor=new Float_t*[tracks]; + for(Int_t i=0;iUncheckedAt(hit); + trk=ahit->GetTrack(); + Float_t x[3]={ahit->X(),ahit->Y(),ahit->Z()}; + Int_t index[3]={1,((AliTPChit*)ahit)->fSector,0}; + Int_t currentrow = fTPCParam->GetPadRow(x,index) ; + if (currentrow!=lastrow){ + lastrow = currentrow; + //later calculate intersection point + if(ntrk[trk]==limi[trk]) { + // + // Initialise a new track + fp=new Float_t[3*(limi[trk]+chunk)]; + if(coor[trk]) { + memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]); + delete [] coor[trk]; + } + limi[trk]+=chunk; + coor[trk] = fp; + } else { + fp = coor[trk]; + } + fp[3*ntrk[trk] ] = ahit->X(); + fp[3*ntrk[trk]+1] = ahit->Y(); + fp[3*ntrk[trk]+2] = ahit->Z(); + ntrk[trk]++; + } + ahit = NextHit2(); + } + // - if (num_of_clusters==MAX_CLUSTER_PER_ROW) { - cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return; + for(trk=0; trkSetMarkerColor(GetMarkerColor()+1); + points->SetMarkerStyle(5); + points->SetMarkerSize(0.2); + points->SetDetector(this); + points->SetParticle(trk); + // points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20); + points->SetPolyMarker(ntrk[trk],coor[trk],30); + fPoints->AddAt(points,tracks+trk); + delete [] coor[trk]; + coor[trk]=0; + } } - if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;} - int i=Find(c->fY); - memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*)); - clusters[i]=c; num_of_clusters++; + delete [] coor; + delete [] ntrk; + delete [] limi; } -//_____________________________________________________________________________ -int AliTPCRow::Find(Double_t y) const + + +void AliTPC::FindTrackHitsIntersection(TClonesArray * arr) { + // + //fill clones array with intersection of current point with the + //middle of the row + Int_t sector; + Int_t ipart; + + const Int_t kcmaxhits=30000; + TVector * xxxx = new TVector(kcmaxhits*4); + TVector & xxx = *xxxx; + Int_t maxhits = kcmaxhits; + // - // - if (y <= clusters[0]->fY) return 0; - if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters; - int b=0, e=num_of_clusters-1, m=(b+e)/2; - for (; b clusters[m]->fY) b=m+1; - else e=m; - } - return m; -} + AliTPChit * tpcHit=0; + tpcHit = (AliTPChit*)FirstHit2(-1); + Int_t currentIndex=0; + Int_t lastrow=-1; //last writen row + + while (tpcHit){ + if (tpcHit==0) continue; + sector=tpcHit->fSector; // sector number + ipart=tpcHit->Track(); + + //find row number + + Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()}; + Int_t index[3]={1,sector,0}; + Int_t currentrow = fTPCParam->GetPadRow(x,index) ; + if (currentrow<0) continue; + if (lastrow<0) lastrow=currentrow; + if (currentrow==lastrow){ + if ( currentIndex>=maxhits){ + maxhits+=kcmaxhits; + xxx.ResizeTo(4*maxhits); + } + xxx(currentIndex*4)=x[0]; + xxx(currentIndex*4+1)=x[1]; + xxx(currentIndex*4+2)=x[2]; + xxx(currentIndex*4+3)=tpcHit->fQ; + currentIndex++; + } + else + if (currentIndex>2){ + Float_t sumx=0; + Float_t sumx2=0; + Float_t sumx3=0; + Float_t sumx4=0; + Float_t sumy=0; + Float_t sumxy=0; + Float_t sumx2y=0; + Float_t sumz=0; + Float_t sumxz=0; + Float_t sumx2z=0; + Float_t sumq=0; + for (Int_t index=0;indexGetNPads(sector,lastrow)-1)/2; + Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+ + sumx2*(sumx*sumx3-sumx2*sumx2); + + Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+ + sumx2*(sumxy*sumx3-sumx2y*sumx2); + Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+ + sumx2*(sumxz*sumx3-sumx2z*sumx2); + + Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+ + sumx2*(sumx*sumx2y-sumx2*sumxy); + Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+ + sumx2*(sumx*sumx2z-sumx2*sumxz); + + Float_t y=detay/det+centralPad; + Float_t z=detaz/det; + Float_t by=detby/det; //y angle + Float_t bz=detbz/det; //z angle + sumy/=Float_t(currentIndex); + sumz/=Float_t(currentIndex); + + AliComplexCluster cl; + cl.fX=z; + cl.fY=y; + cl.fQ=sumq; + cl.fSigmaX2=bz; + cl.fSigmaY2=by; + cl.fTracks[0]=ipart; + + AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow)); + if (row!=0) row->InsertCluster(&cl); + currentIndex=0; + lastrow=currentrow; + } //end of calculating cluster for given row + + } // end of loop over hits + xxxx->Delete(); + + +}