X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPC.cxx;h=fd11f41d0c7af43193ae54ce9410db1c3ed5afde;hb=b96715746b109eb6dd6231a91c1abc8bac34ef89;hp=3fd91e2fb7a0c08aa2e93f5d7760e74aa107755e;hpb=1578254f491848aff213be55a70a61d195dd001a;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPC.cxx b/TPC/AliTPC.cxx index 3fd91e2fb7a..fd11f41d0c7 100644 --- a/TPC/AliTPC.cxx +++ b/TPC/AliTPC.cxx @@ -1,3 +1,20 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* $Id$ */ + /////////////////////////////////////////////////////////////////////////////// // // // Time Projection Chamber // @@ -14,30 +31,55 @@ // // /////////////////////////////////////////////////////////////////////////////// +// + +#include +#include + +#include +#include +#include #include -#include +#include #include -#include -#include #include -#include #include -#include "TParticle.h" -#include "AliTPC.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliDigits.h" +#include "AliMagF.h" +#include "AliPoints.h" #include "AliRun.h" -#include -#include -#include "AliMC.h" - -//MI change -#include "AliTPCParam.h" -#include "AliTPCD.h" +#include "AliRunLoader.h" +#include "AliSimDigits.h" +#include "AliTPC.h" +#include "AliTPC.h" +#include "AliTPCDigitsArray.h" +#include "AliTPCLoader.h" #include "AliTPCPRF2D.h" +#include "AliTPCParamSR.h" #include "AliTPCRF1D.h" - +//#include "AliTPCTrackHits.h" +#include "AliTPCTrackHitsV2.h" +#include "AliTrackReference.h" +#include "AliMC.h" +#include "AliTPCDigitizer.h" +#include "AliTPCBuffer.h" +#include "AliTPCDDLRawData.h" +#include "AliLog.h" +#include "AliRawReader.h" +#include "AliTPCRawStream.h" ClassImp(AliTPC) - //_____________________________________________________________________________ AliTPC::AliTPC() { @@ -45,16 +87,23 @@ 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; + fDefaults = 0; + fTrackHits = 0; + // fTrackHitsOld = 0; +#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1) + fHitType = 4; // ROOT containers +#else + fHitType = 2; //default CONTAINERS - based on ROOT structure +#endif + fTPCParam = 0; + fNoiseTable = 0; + fActiveSectors =0; + fSens = 0; + } //_____________________________________________________________________________ @@ -68,1204 +117,1345 @@ 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->GetMCApp()->AddHitList(fHits); + fDigitsArray = 0; + fDefaults = 0; + // + fTrackHits = new AliTPCTrackHitsV2; + fTrackHits->SetHitPrecision(0.002); + fTrackHits->SetStepPrecision(0.003); + fTrackHits->SetMaxDistance(100); + + //fTrackHitsOld = new AliTPCTrackHits; //MI - 13.09.2000 + //fTrackHitsOld->SetHitPrecision(0.002); + //fTrackHitsOld->SetStepPrecision(0.003); + //fTrackHitsOld->SetMaxDistance(100); + + fNoiseTable =0; + +#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1) + fHitType = 4; // ROOT containers +#else + fHitType = 2; +#endif + fActiveSectors = 0; // // 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); -} -//_____________________________________________________________________________ -AliTPC::~AliTPC() -{ // - // TPC destructor + // Set TPC parameters // - fIshunt = 0; - delete fHits; - delete fDigits; - delete fClusters; - delete fTracks; - delete fDigParam; - if (fDigitsIndex) delete [] fDigitsIndex; - if (fClustersIndex) delete [] fClustersIndex; -} -//_____________________________________________________________________________ -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); + + if (!strcmp(title,"Default")) { + fTPCParam = new AliTPCParamSR; + } else { + AliWarning("In Config.C you must set non-default parameters."); + fTPCParam=0; + } + fSens = 0; + } - + //_____________________________________________________________________________ -void AliTPC::AddCluster(const AliTPCcluster &c) -{ +AliTPC::AliTPC(const AliTPC& t):AliDetector(t){ // - // Add a simulated cluster copy to the list + // dummy copy constructor // - if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000); - TClonesArray &lclusters = *fClusters; - new(lclusters[fNclusters++]) AliTPCcluster(c); } - -//_____________________________________________________________________________ -void AliTPC::AddDigit(Int_t *tracks, Int_t *digits) +AliTPC::~AliTPC() { // - // Add a TPC digit to the list + // TPC destructor // - // TClonesArray &ldigits = *fDigits; - //MI change - TClonesArray &ldigits = *fDigParam->GetArray(); - new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits); + + fIshunt = 0; + delete fHits; + delete fDigits; + delete fTPCParam; + delete fTrackHits; //MI 15.09.2000 + // delete fTrackHitsOld; //MI 10.12.2001 + + fDigitsArray = 0x0; + delete [] fNoiseTable; + delete [] fActiveSectors; + } - + //_____________________________________________________________________________ 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); -} - -//_____________________________________________________________________________ -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); + if (fHitType&1){ + TClonesArray &lhits = *fHits; + new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits); + } + if (fHitType>1) + AddHit2(track,vol,hits); } //_____________________________________________________________________________ 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 ) +void AliTPC::CreateMaterials() { - // - // Calculate distance from TPC to mouse on the display - // Dummy procedure - // - return 9999; -} + //----------------------------------------------- + // Create Materials for for TPC simulations + //----------------------------------------------- -//_____________________________________________________________________________ -//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.; + //----------------------------------------------------------------- + // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl + //----------------------------------------------------------------- + Int_t iSXFLD=gAlice->Field()->Integ(); + Float_t sXMGMX=gAlice->Field()->Max(); -//_____________________________________________________________________________ -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; -} + Float_t amat[5]; // atomic numbers + Float_t zmat[5]; // z + Float_t wmat[5]; // proportions -//_____________________________________________________________________________ -static Double_t SigmaZ2(Double_t r, Double_t tgl) -{ - // - // 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; -} + Float_t density; + -//_____________________________________________________________________________ -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)); - Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); + //***************** Gases ************************* - return -xr*yr/sqrt(xr*xr+yr*yr); -} + + //-------------------------------------------------------------- + // gases - air and CO2 + //-------------------------------------------------------------- + // CO2 -//_____________________________________________________________________________ -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)); + amat[0]=12.011; + amat[1]=15.9994; - Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); - - return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); -} + zmat[0]=6.; + zmat[1]=8.; -//_____________________________________________________________________________ -inline Double_t f3(Double_t x1,Double_t y1, //tgl - Double_t x2,Double_t y2, - Double_t z1,Double_t z2) -{ + wmat[0]=0.2729; + wmat[1]=0.7271; + + density=0.001977; + + + AliMixture(10,"CO2",amat,zmat,density,2,wmat); // - // Function f3 + // Air // - return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); -} - -//_____________________________________________________________________________ -static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec, - int s, int ri, int rf=0) -{ + amat[0]=15.9994; + amat[1]=14.007; + // + zmat[0]=8.; + zmat[1]=7.; // - // Propagate track + wmat[0]=0.233; + wmat[1]=0.767; // - int try_again=ROWS_TO_SKIP; - Double_t alpha=sec->GetAlpha(); - int ns=int(2*TMath::Pi()/alpha)+1; + density=0.001205; - 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; + AliMixture(11,"Air",amat,zmat,density,2,wmat); + + //---------------------------------------------------------------- + // drift gases + //---------------------------------------------------------------- - 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(); - 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]= 20.18; + amat[1]=12.011; + amat[2]=15.9994; + amat[3]=14.007; - return s; -} + zmat[0]= 10.; + zmat[1]=6.; + zmat[2]=8.; + zmat[3]=7.; + wmat[0]=0.7707; + wmat[1]=0.0539; + wmat[2]=0.1438; + wmat[3]=0.0316; + + density=0.0010252; -//_____________________________________________________________________________ -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)); + AliMixture(12,"Ne-CO2-N-1",amat,zmat,density,4,wmat); + AliMixture(13,"Ne-CO2-N-2",amat,zmat,density,4,wmat); - 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); - } - } - } -} + //---------------------------------------------------------------------- + // solid materials + //---------------------------------------------------------------------- -//_____________________________________________________________________________ -void AliTPC::Clusters2Tracks() -{ - // - // TPC Track finder from clusters. - // - if (!fClusters) return; - AliTPCParam *p=&fDigParam->GetParam(); - Int_t nrow_low=p->GetNRowLow(); - Int_t nrow_up=p->GetNRowUp(); + // Kevlar C14H22O2N2 - AliTPCSSector ssec[S_MAXSEC/2]; - for (int i=0; iGetEntriesFast(); - while (ncl--) { - AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl); + wmat[0] = 14.; + wmat[1] = 22.; + wmat[2] = 2.; + wmat[3] = 2.; - 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); + AliMedium(4,"Al",23,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(5,"Kevlar",14,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(6,"Nomex",15,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(7,"Makrolon",16,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(8,"Mylar",18,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(9,"Tedlar",17,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + // + AliMedium(10,"Prepreg1",19,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(11,"Prepreg2",20,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(12,"Prepreg3",21,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + AliMedium(13,"Epoxy",26,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); + + AliMedium(14,"Cu",25,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(15,"Si",24,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(16,"G10",22,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(17,"Plexiglas",27,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(18,"Steel",29,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - 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.); + // + //Generate table with noise + // + if (fTPCParam==0) { + // error message + fNoiseDepth=0; + return; + } + if (fNoiseTable) delete[] fNoiseTable; + fNoiseTable = new Float_t[tablesize]; + fNoiseDepth = tablesize; + fCurrentNoise =0; //!index of the noise in the noise table - 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 ); + Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac(); + for (Int_t i=0;iGaus(0,norm); } -//_____________________________________________________________________________ -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;} +Float_t AliTPC::GetNoise() +{ + // get noise from table + // if ((fCurrentNoise%10)==0) + // fCurrentNoise= gRandom->Rndm()*fNoiseDepth; + if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0; + return fNoiseTable[fCurrentNoise++]; + //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); +} -//_____________________________________________________________________________ -static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c) +Bool_t AliTPC::IsSectorActive(Int_t sec) const { // - // Find clusters - // - Bin& b=bins[i][j]; - double q=double(b.dig->fSignal); + // check if the sector is active + if (!fActiveSectors) return kTRUE; + else return fActiveSectors[sec]; +} + +void AliTPC::SetActiveSectors(Int_t * sectors, Int_t n) +{ + // activate interesting sectors + SetTreeAddress();//just for security + if (fActiveSectors) delete [] fActiveSectors; + fActiveSectors = new Bool_t[fTPCParam->GetNSector()]; + for (Int_t i=0;iGetNSector();i++) fActiveSectors[i]=kFALSE; + for (Int_t i=0;i=0) && sectors[i]GetNSector()) fActiveSectors[sectors[i]]=kTRUE; + +} - 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++; +void AliTPC::SetActiveSectors(Int_t flag) +{ + // + // activate sectors which were hitted by tracks + //loop over tracks + SetTreeAddress();//just for security + if (fHitType==0) return; // if Clones hit - not short volume ID information + if (fActiveSectors) delete [] fActiveSectors; + fActiveSectors = new Bool_t[fTPCParam->GetNSector()]; + if (flag) { + for (Int_t i=0;iGetNSector();i++) fActiveSectors[i]=kTRUE; + return; } + for (Int_t i=0;iGetNSector();i++) fActiveSectors[i]=kFALSE; + TBranch * branch=0; + if (TreeH() == 0x0) + { + AliFatal("Can not find TreeH in folder"); + return; + } + if (fHitType>1) branch = TreeH()->GetBranch("TPC2"); + else branch = TreeH()->GetBranch("TPC"); + Stat_t ntracks = TreeH()->GetEntries(); + // loop over all hits + AliDebug(1,Form("Got %d tracks",ntracks)); - 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; + for(Int_t track=0;trackGetBranch("fVolumes"); + TBranch * br2 = TreeH()->GetBranch("fNVolumes"); + br1->GetEvent(track); + br2->GetEvent(track); + Int_t *volumes = fTrackHits->GetVolumes(); + for (Int_t j=0;jGetNVolumes(); j++) + fActiveSectors[volumes[j]]=kTRUE; + } + + // +// if (fTrackHitsOld && fHitType&2) { +// TBranch * br = TreeH()->GetBranch("fTrackHitsInfo"); +// br->GetEvent(track); +// AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo; +// for (UInt_t j=0;jGetSize();j++){ +// fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE; +// } +// } + } +} + + - 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); - -} //_____________________________________________________________________________ -void AliTPC::Digits2Clusters() +void AliTPC::Digits2Raw() { - // - // 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; +// convert digits of the current event to raw data + + static const Int_t kThreshold = 0; + + fLoader->LoadDigits(); + TTree* digits = fLoader->TreeD(); + if (!digits) { + AliError("No digits tree"); + return; + } + + AliSimDigits digarr; + AliSimDigits* digrow = &digarr; + digits->GetBranch("Segment")->SetAddress(&digrow); + + const char* fileName = "AliTPCDDL.dat"; + AliTPCBuffer* buffer = new AliTPCBuffer(fileName); + //Verbose level + // 0: Silent + // 1: cout messages + // 2: txt files with digits + //BE CAREFUL, verbose level 2 MUST be used only for debugging and + //it is highly suggested to use this mode only for debugging digits files + //reasonably small, because otherwise the size of the txt files can reach + //quickly several MB wasting time and disk space. + buffer->SetVerbose(0); + + Int_t nEntries = Int_t(digits->GetEntries()); + Int_t previousSector = -1; + Int_t subSector = 0; + for (Int_t i = 0; i < nEntries; i++) { + digits->GetEntry(i); + Int_t sector, row; + fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row); + if(previousSector != sector) { + subSector = 0; + previousSector = sector; + } + + if (sector < 36) { //inner sector [0;35] + if (row != 30) { + //the whole row is written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0, + sector, subSector, row); + } else { + //only the pads in the range [37;48] are written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1, + sector, subSector, row); + subSector = 1; + //only the pads outside the range [37;48] are written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2, + sector, subSector, row); + }//end else + + } else { //outer sector [36;71] + if (row == 54) subSector = 2; + if ((row != 27) && (row != 76)) { + buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0, + sector, subSector, row); + } else if (row == 27) { + //only the pads outside the range [43;46] are written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2, + sector, subSector, row); + subSector = 1; + //only the pads in the range [43;46] are written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1, + sector, subSector, row); + } else if (row == 76) { + //only the pads outside the range [33;88] are written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2, + sector, subSector, row); + subSector = 3; + //only the pads in the range [33;88] are written into the output file + buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1, + sector, subSector, row); + } + }//end else + }//end for + + delete buffer; + fLoader->UnloadDigits(); + + AliTPCDDLRawData rawWriter; + rawWriter.SetVerbose(0); + + rawWriter.RawData(fileName); + gSystem->Unlink(fileName); + +} + + +//_____________________________________________________________________________ +Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){ + // Converts the TPC raw data into summable digits + // The method is used for merging simulated and + // real data events + if (fLoader->TreeS() == 0x0 ) { + fLoader->MakeTree("S"); + } + + if(fDefaults == 0) SetDefaults(); // check if the parameters are set + + //setup TPCDigitsArray + if(GetDigitsArray()) delete GetDigitsArray(); + + AliTPCDigitsArray *arr = new AliTPCDigitsArray; + arr->SetClass("AliSimDigits"); + arr->Setup(fTPCParam); + arr->MakeTree(fLoader->TreeS()); + + SetDigitsArray(arr); + + // set zero suppression to "0" + fTPCParam->SetZeroSup(0); + + // Loop over sectors + const Int_t kmaxTime = fTPCParam->GetMaxTBin(); + const Int_t kNIS = fTPCParam->GetNInnerSector(); + const Int_t kNOS = fTPCParam->GetNOuterSector(); + const Int_t kNS = kNIS + kNOS; + + Short_t** allBins = NULL; //array which contains the data for one sector - 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(); + for(Int_t iSector = 0; iSector < kNS; iSector++) { - 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); + Int_t nRows = fTPCParam->GetNRow(iSector); + Int_t nDDLs = 0, indexDDL = 0; + if (iSector < kNIS) { + nDDLs = 2; + indexDDL = iSector * 2; } - - int ndig; - for (ndig=0; ndigUncheckedAt(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; + else { + nDDLs = 4; + indexDDL = (iSector-kNIS) * 4 + kNIS * 2; } - - 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++; - } + + // Loas the raw data for corresponding DDLs + rawReader->Reset(); + AliTPCRawStream input(rawReader); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + + // Alocate and init the array with the sector data + allBins = new Short_t*[nRows]; + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad = fTPCParam->GetNPads(iSector,iRow); + Int_t maxBin = kmaxTime*maxPad; + allBins[iRow] = new Short_t[maxBin]; + memset(allBins[iRow],0,sizeof(Short_t)*maxBin); } + + // Begin loop over altro data + while (input.Next()) { + + if (input.GetSector() != iSector) + AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSector,input.GetSector())); + + Int_t iRow = input.GetRow(); + if (iRow < 0 || iRow >= nRows) + AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !", + iRow, 0, nRows -1)); + Int_t iPad = input.GetPad(); + + Int_t maxPad = fTPCParam->GetNPads(iSector,iRow); + + if (iPad < 0 || iPad >= maxPad) + AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !", + iPad, 0, maxPad -1)); + + Int_t iTimeBin = input.GetTime(); + if ( iTimeBin < 0 || iTimeBin >= kmaxTime) + AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !", + iTimeBin, 0, kmaxTime -1)); + + Int_t maxBin = kmaxTime*maxPad; + + if (((iPad*kmaxTime+iTimeBin) >= maxBin) || + ((iPad*kmaxTime+iTimeBin) < 0)) + AliFatal(Form("Index outside the allowed range" + " Sector=%d Row=%d Pad=%d Timebin=%d" + " (Max.index=%d)",iSector,iRow,iPad,iTimeBin,maxBin)); + + allBins[iRow][iPad*kmaxTime+iTimeBin] = input.GetSignal(); + + } // End loop over altro data - for (ndig=0; ndigUncheckedAt(ndig); - if (TMath::Abs(dig->fSignal) >= 0) - bins[dig->fPad][dig->fTime].dig=dig; + // Now fill the digits array + if (fDigitsArray->GetTree()==0) { + AliFatal("Tree not set in fDigitsArray"); } - - 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); + + for (Int_t iRow = 0; iRow < nRows; iRow++) { + AliDigits * dig = fDigitsArray->CreateRow(iSector,iRow); + + Int_t maxPad = fTPCParam->GetNPads(iSector,iRow); + for(Int_t iPad = 0; iPad < maxPad; iPad++) { + for(Int_t iTimeBin = 0; iTimeBin < kmaxTime; iTimeBin++) { + Short_t q = allBins[iRow][iPad*kmaxTime + iTimeBin]; + if (q <= 0) continue; + q *= 16; + dig->SetDigitFast((Short_t)q,iTimeBin,iPad); } } - } - - cerr<<"sector, row, digits, clusters: " - <Clear(); - - } -} + fDigitsArray->StoreRow(iSector,iRow); + Int_t ndig = dig->GetDigitSize(); + + AliDebug(10, + Form("*** Sector, row, compressed digits %d %d %d ***\n", + iSector,iRow,ndig)); + + fDigitsArray->ClearRow(iSector,iRow); -//_____________________________________________________________________________ -void AliTPC::ElDiff(Float_t *xyz) -{ - //-------------------------------------------------- - // calculates the diffusion of a single electron - //-------------------------------------------------- + } // end of the sector digitization - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t driftl; - - Float_t z0=xyz[2]; + for (Int_t iRow = 0; iRow < nRows; iRow++) + delete [] allBins[iRow]; - driftl=z_end-TMath::Abs(xyz[2]); + delete [] allBins; - if(driftl<0.01) driftl=0.01; + } - // check the attachment + fLoader->WriteSDigits("OVERWRITE"); - driftl=TMath::Sqrt(driftl); + if(GetDigitsArray()) delete GetDigitsArray(); + SetDigitsArray(0x0); - // 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); - - 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); - } + return kTRUE; } -//_____________________________________________________________________________ -void AliTPC::Hits2Clusters() +//______________________________________________________________________ +AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const { - //-------------------------------------------------------- - // TPC simple cluster generator from hits - // obtained from the TPC Fast Simulator - // The point errors are taken from the parametrization - //-------------------------------------------------------- - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t sigma_rphi,sigma_z,cl_rphi,cl_z; - // - 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 ipart; - Float_t xyz[5]; - Float_t pl,pt,tanth,rpad,ratio; - Float_t cph,sph; - - //--------------------------------------------------------------- - // Get the access to the tracks - //--------------------------------------------------------------- - - TTree *TH = gAlice->TreeH(); - Stat_t ntracks = TH->GetEntries(); + return new AliTPCDigitizer(manager); +} +//__ +void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/) +{ + //create digits from summable digits + GenerNoise(500000); //create teble with noise + + //conect tree with sSDigits + TTree *t = fLoader->TreeS(); + + if (t == 0x0) { + fLoader->LoadSDigits("READ"); + t = fLoader->TreeS(); + if (t == 0x0) { + AliError("Can not get input TreeS"); + return; + } + } - //------------------------------------------------------------ - // 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 - // - // First cluster for sector 1 starts at "0" - //------------------------------------------------------------ + if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D"); + AliSimDigits digarr, *dummy=&digarr; + TBranch* sdb = t->GetBranch("Segment"); + if (sdb == 0x0) { + AliError("Can not find branch with segments in TreeS."); + return; + } + + sdb->SetAddress(&dummy); + + Stat_t nentries = t->GetEntries(); + + // set zero suppression + + fTPCParam->SetZeroSup(2); + + // get zero suppression + + Int_t zerosup = fTPCParam->GetZeroSup(); + + //make tree with digits - fClustersIndex[0] = 0; + AliTPCDigitsArray *arr = new AliTPCDigitsArray; + arr->SetClass("AliSimDigits"); + arr->Setup(fTPCParam); + arr->MakeTree(fLoader->TreeD()); - // - for(Int_t isec=1;isecAdjustAngles(isec,cph,sph); - - //------------------------------------------------------------ - // Loop over tracks - //------------------------------------------------------------ + AliTPCParam * par = fTPCParam; + + //Loop over segments of the TPC + + for (Int_t n=0; nGetEvent(n); + Int_t sec, row; + if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) { + AliWarning(Form("Invalid segment ID ! %d",digarr.GetID())); + continue; + } + if (!IsSectorActive(sec)) continue; - for(Int_t track=0;trackGetEvent(track); - // - // Get number of the TPC hits and a pointer - // to the particles - // - nhits=fHits->GetEntriesFast(); - Particles=gAlice->Particles(); - // - // Loop over hits - // - for(Int_t hit=0;hitUncheckedAt(hit); - sector=tpcHit->fSector; // sector number - if(sector != isec) continue; //terminate iteration - ipart=tpcHit->fTrack; - particle=(TParticle*)Particles->UncheckedAt(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); - 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 ); - - // 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; - - // 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; - - // - - // - // 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); - - } // end of loop over hits - } // end of loop over tracks + AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row); + Int_t nrows = digrow->GetNRows(); + Int_t ncols = digrow->GetNCols(); + + digrow->ExpandBuffer(); + digarr.ExpandBuffer(); + digrow->ExpandTrackBuffer(); + digarr.ExpandTrackBuffer(); + - fClustersIndex[isec] = fNclusters; // update clusters index + Short_t * pamp0 = digarr.GetDigits(); + Int_t * ptracks0 = digarr.GetTracks(); + Short_t * pamp1 = digrow->GetDigits(); + Int_t * ptracks1 = digrow->GetTracks(); + Int_t nelems =nrows*ncols; + Int_t saturation = fTPCParam->GetADCSat() - 1; + //use internal structure of the AliDigits - for speed reason + //if you cahnge implementation + //of the Alidigits - it must be rewriten - + for (Int_t i= 0; izerosup){ + if (q>saturation) q=saturation; + *pamp1=(Short_t)q; + + ptracks1[0]=ptracks0[0]; + ptracks1[nelems]=ptracks0[nelems]; + ptracks1[2*nelems]=ptracks0[2*nelems]; + } + pamp0++; + pamp1++; + ptracks0++; + ptracks1++; + } + + arr->StoreRow(sec,row); + arr->ClearRow(sec,row); + } + - } // end of loop over sectors - - fClustersIndex[fNsectors]--; // set end of the clusters buffer - + //write results + fLoader->WriteDigits("OVERWRITE"); + + delete arr; } +//__________________________________________________________________ +void AliTPC::SetDefaults(){ + // + // setting the defaults + // + + // Set response functions + // + AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName()); + rl->CdGAFile(); + AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60"); + if(param){ + AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits..."); + delete param; + param = new AliTPCParamSR(); + } + else { + param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60"); + } + if(!param){ + AliFatal("No TPC parameters found"); + } -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;isecSetGauss(param->GetZSigma(),param->GetZWidth(),1.); + rf->SetOffset(3*param->GetZSigma()); + rf->Update(); + + TDirectory *savedir=gDirectory; + TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root"); + if (!f->IsOpen()) + AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !"); + TString s; + prfinner->Read("prf_07504_Gati_056068_d02"); + //PH Set different names + s=prfinner->GetGRF()->GetName(); + s+="in"; + prfinner->GetGRF()->SetName(s.Data()); -//_____________________________________________________________________________ -void AliTPC::Hits2DigitsSector(Int_t isec) -{ - //------------------------------------------------------------------- - // TPC conversion from hits to digits. - //------------------------------------------------------------------- + prfouter1->Read("prf_10006_Gati_047051_d03"); + s=prfouter1->GetGRF()->GetName(); + s+="out1"; + prfouter1->GetGRF()->SetName(s.Data()); - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- + prfouter2->Read("prf_15006_Gati_047051_d03"); + s=prfouter2->GetGRF()->GetName(); + s+="out2"; + prfouter2->GetGRF()->SetName(s.Data()); - //------------------------------------------------------- - // Get the access to the track hits - //------------------------------------------------------- + f->Close(); + savedir->cd(); - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - TTree *TH = gAlice->TreeH(); // pointer to the hits tree - Stat_t ntracks = TH->GetEntries(); + param->SetInnerPRF(prfinner); + param->SetOuter1PRF(prfouter1); + param->SetOuter2PRF(prfouter2); + param->SetTimeRF(rf); - if( ntracks > 0){ + // set fTPCParam - //------------------------------------------- - // Only if there are any tracks... - //------------------------------------------- + SetParam(param); - - // TObjArrays for three neighouring pad-rows - TObjArray **rowTriplet = new TObjArray* [3]; - - // TObjArray-s for each pad-row + fDefaults = 1; - TObjArray **row; - - for (Int_t trip=0;trip<3;trip++){ - rowTriplet[trip]=new TObjArray; - } +} +//__________________________________________________________________ +void AliTPC::Hits2Digits() +{ + // + // creates digits from hits + // + if (!fTPCParam->IsGeoRead()){ + // + // read transformation matrices for gGeoManager + // + fTPCParam->ReadGeoMatrices(); + } + fLoader->LoadHits("read"); + fLoader->LoadDigits("recreate"); + AliRunLoader* runLoader = fLoader->GetRunLoader(); - - printf("*** Processing sector number %d ***\n",isec); + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + SetActiveSectors(); + Hits2Digits(iEvent); + } - Int_t nrows =fTPCParam->GetNRow(isec); + fLoader->UnloadHits(); + fLoader->UnloadDigits(); +} +//__________________________________________________________________ +void AliTPC::Hits2Digits(Int_t eventnumber) +{ + //---------------------------------------------------- + // Loop over all sectors for a single event + //---------------------------------------------------- + AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName()); + rl->GetEvent(eventnumber); + if (fLoader->TreeH() == 0x0) { + if(fLoader->LoadHits()) { + AliError("Can not load hits."); + } + } + SetTreeAddress(); + + if (fLoader->TreeD() == 0x0 ) { + fLoader->MakeTree("D"); + if (fLoader->TreeD() == 0x0 ) { + AliError("Can not get TreeD"); + return; + } + } - row= new TObjArray* [nrows]; - - MakeSector(isec,nrows,TH,ntracks,row); + if(fDefaults == 0) SetDefaults(); // check if the parameters are set + GenerNoise(500000); //create teble with noise - //-------------------------------------------------------- - // Digitize this sector, row by row - // row[i] is the pointer to the TObjArray of TVectors, - // each one containing electrons accepted on this - // row, assigned into tracks - //-------------------------------------------------------- + //setup TPCDigitsArray - Int_t i; + if(GetDigitsArray()) delete GetDigitsArray(); + + AliTPCDigitsArray *arr = new AliTPCDigitsArray; + arr->SetClass("AliSimDigits"); + arr->Setup(fTPCParam); - for (i=0;iMakeTree(fLoader->TreeD()); + SetDigitsArray(arr); - // Triplets for i = 0 and i=1 are identical! - // The same for i = nrows-1 and nrows! + fDigitsSwitch=0; // standard digits - if(i != 1 && i != nrows-1){ - MakeTriplet(i,rowTriplet,row); - } + for(Int_t isec=0;isecGetNSector();isec++) + if (IsSectorActive(isec)) { + AliDebug(1,Form("Hits2Digits","Sector %d is active.",isec)); + Hits2DigitsSector(isec); + } + else { + AliDebug(1,Form("Hits2Digits","Sector %d is NOT active.",isec)); + } + + fLoader->WriteDigits("OVERWRITE"); + +//this line prevents the crash in the similar one +//on the beginning of this method +//destructor attempts to reset the tree, which is deleted by the loader +//need to be redesign + if(GetDigitsArray()) delete GetDigitsArray(); + SetDigitsArray(0x0); + +} + +//__________________________________________________________________ +void AliTPC::Hits2SDigits2(Int_t eventnumber) +{ + + //----------------------------------------------------------- + // summable digits - 16 bit "ADC", no noise, no saturation + //----------------------------------------------------------- - DigitizeRow(i,isec,rowTriplet); + //---------------------------------------------------- + // Loop over all sectors for a single event + //---------------------------------------------------- - fDigParam->Fill(); + AliRunLoader* rl = fLoader->GetRunLoader(); - Int_t ndig=fDigParam->GetArray()->GetEntriesFast(); + rl->GetEvent(eventnumber); + if (fLoader->TreeH() == 0x0) { + if(fLoader->LoadHits()) { + AliError("Can not load hits."); + return; + } + } + SetTreeAddress(); - printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig); - ResetDigits(); // reset digits for this row after storing them - - } // end of the sector digitization - - // delete the last triplet + if (fLoader->TreeS() == 0x0 ) { + fLoader->MakeTree("S"); + } + + if(fDefaults == 0) SetDefaults(); + + GenerNoise(500000); //create table with noise + //setup TPCDigitsArray - for (i=0;i<3;i++) rowTriplet[i]->Delete(); - - 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 - //------------------------------------------------------------------ + if(GetDigitsArray()) delete GetDigitsArray(); - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- + + AliTPCDigitsArray *arr = new AliTPCDigitsArray; + arr->SetClass("AliSimDigits"); + arr->Setup(fTPCParam); + arr->MakeTree(fLoader->TreeS()); - AliTPCParam * fTPCParam = &(fDigParam->GetParam()); - Float_t gasgain = fTPCParam->GetGasGain(); - Int_t nTracks[3]; + SetDigitsArray(arr); - Int_t nElements,nElectrons; + fDigitsSwitch=1; // summable digits + + // set zero suppression to "0" - TVector *pv; - TVector *tv; + fTPCParam->SetZeroSup(0); - //------------------------------------------------------------------- - // pv is an "old" track, i.e. label + triplets of (x,y,z) - // for each electron - // - //------------------------------------------------------------------- + for(Int_t isec=0;isecGetNSector();isec++) + if (IsSectorActive(isec)) { + Hits2DigitsSector(isec); + } + fLoader->WriteSDigits("OVERWRITE"); - Int_t i1,i2; - Int_t nel,nt; +//this line prevents the crash in the similar one +//on the beginning of this method +//destructor attempts to reset the tree, which is deleted by the loader +//need to be redesign + if(GetDigitsArray()) delete GetDigitsArray(); + SetDigitsArray(0x0); +} +//__________________________________________________________________ - if(row == 0 || row == 1){ +void AliTPC::Hits2SDigits() +{ - // create entire triplet for the first AND the second row + //----------------------------------------------------------- + // summable digits - 16 bit "ADC", no noise, no saturation + //----------------------------------------------------------- - nTracks[0] = prow[0]->GetEntries(); - nTracks[1] = prow[1]->GetEntries(); - nTracks[2] = prow[2]->GetEntries(); + if (!fTPCParam->IsGeoRead()){ + // + // read transformation matrices for gGeoManager + // + fTPCParam->ReadGeoMatrices(); + } + + fLoader->LoadHits("read"); + fLoader->LoadSDigits("recreate"); + AliRunLoader* runLoader = fLoader->GetRunLoader(); - for(i1=0;i1<3;i1++){ - nt = nTracks[i1]; // number of tracks for this row + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + SetTreeAddress(); + SetActiveSectors(); + Hits2SDigits2(iEvent); + } - for(i2=0;i2At(i2); - TVector &v1 = *pv; - nElements = pv->GetNrows(); - nElectrons = (nElements-1)/3; + fLoader->UnloadHits(); + fLoader->UnloadSDigits(); +} +//_____________________________________________________________________________ - tv = new TVector(4*nElectrons+1); // create TVector for a modified track - TVector &v2 = *tv; - v2(0)=v1(0); //track label +void AliTPC::Hits2DigitsSector(Int_t isec) +{ + //------------------------------------------------------------------- + // TPC conversion from hits to digits. + //------------------------------------------------------------------- - 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 - // + //----------------------------------------------------------------- + // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl + //----------------------------------------------------------------- - rowTriplet[i1]->Add(tv); + //------------------------------------------------------- + // Get the access to the track hits + //------------------------------------------------------- + // check if the parameters are set - important if one calls this method + // directly, not from the Hits2Digits - } // end of loop over tracks for this row + if(fDefaults == 0) SetDefaults(); - prow[i1]->Delete(); // remove "old" tracks - delete prow[i1]; // delete TObjArray for this row - prow[i1]=0; // set pointer to NULL + TTree *tH = TreeH(); // pointer to the hits tree + if (tH == 0x0) { + AliFatal("Can not find TreeH in folder"); + return; + } - } // end of loop over row triplets + Stat_t ntracks = tH->GetEntries(); + if( ntracks > 0){ - } - else{ - - rowTriplet[0]->Delete(); // remove old lower row + //------------------------------------------- + // Only if there are any tracks... + //------------------------------------------- - nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row - nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row - nTracks[2]=prow[row+1]->GetEntries(); // next row + TObjArray **row; + Int_t nrows =fTPCParam->GetNRow(isec); - //------------------------------------------- - // shift new tracks downwards - //------------------------------------------- + row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk + + MakeSector(isec,nrows,tH,ntracks,row); - for(i1=0;i1At(i1); - rowTriplet[0]->Add(pv); - } + //-------------------------------------------------------- + // Digitize this sector, row by row + // row[i] is the pointer to the TObjArray of TVectors, + // each one containing electrons accepted on this + // row, assigned into tracks + //-------------------------------------------------------- - rowTriplet[1]->Clear(); // leave tracks on the heap!!! + Int_t i; - for(i1=0;i1At(i1); - rowTriplet[1]->Add(pv); + if (fDigitsArray->GetTree()==0) { + AliFatal("Tree not set in fDigitsArray"); } - rowTriplet[2]->Clear(); // leave tracks on the heap!!! + for (i=0;iCreateRow(isec,i); - //--------------------------------------------- - // Create new upper row - //--------------------------------------------- + DigitizeRow(i,isec,row); - + fDigitsArray->StoreRow(isec,i); - for(i1=0;i1At(i1); - TVector &v1 = *pv; - nElements = pv->GetNrows(); - nElectrons = (nElements-1)/3; + Int_t ndig = dig->GetDigitSize(); + + AliDebug(10, + Form("*** Sector, row, compressed digits %d %d %d ***\n", + isec,i,ndig)); + + fDigitsArray->ClearRow(isec,i); - tv = new TVector(4*nElectrons+1); // create TVector for a modified track - TVector &v2 = *tv; - v2(0)=v1(0); //track label + + } // end of the sector digitization - for(nel=0;nelDelete(); + delete row[i]; + } + + delete [] row; // delete the array of pointers to TObjArray-s - Int_t idx1 = nel*3; - Int_t idx2 = nel*4; - // Avalanche, including fluctuations - Int_t aval = (Int_t) - (-gasgain*TMath::Log(gRandom->Rndm())); - - 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 - - } + } // ntracks >0 -} // 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()); +} // end of Hits2DigitsSector - 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,157 +1464,98 @@ 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); - - 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 + fCurrentIndex[1]= isec; + - - 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 + // + TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated + TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single // - - TMatrix &Total = *m1; + TMatrixF &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); - - } - else{ - - // cross-talk from both inner and outer rows - - GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner - GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer + Int_t i1; + for(lp=0;lpGetEntries(); + 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); + } } - - Total += Cross; // add the cross-talk - - // - // Convert analog signal to ADC counts - // - + Int_t tracks[3]; - Int_t digits[5]; - - - for(Int_t ip=1;ipGetRow(isec,irow); + Int_t gi=-1; + Float_t fzerosup = zerosup+0.5; + for(Int_t it=0;it= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1; // 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 { + if(q <= 0.) continue; // do not fill zeros + if(q>2000.) q=2000.; + q *= 16.; + q = TMath::Nint(q); + } // // "real" signal or electronic noise (list = -1)? // for(Int_t j1=0;j1<3;j1++){ - tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1; + tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2; } - 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 +1564,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 = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2; - 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 + Int_t nElectrons = (tv->GetNrows()-1)/5; + indexRange[0]=9999; // min pad + indexRange[1]=-1; // max pad + indexRange[2]=9999; //min time + indexRange[3]=-1; // max time + TMatrixF &signal = *m1; + TMatrixF &total = *m2; // // Loop over all electrons // - for(Int_t nel=0; nelGetTotalNormFac(); + Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)}; + Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]); + Int_t *index = fTPCParam->GetResBin(0); + Float_t *weight = & (fTPCParam->GetResWeight(0)); - 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 pad-range , lost at the sector edge - - Float_t aval = (absyGetPadPitchWidth(); - for (Int_t i=0;i<7;i++){ - PadSignal[i]=fPRF2D->GetPRF(dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval; - PadSignal[i] *= fTPCParam->GetPadCoupling(); - } + if (n>0) for (Int_t i =0; i=0){ + Int_t time=index[2]; + Float_t qweight = *(weight)*eltoadcfac; + + if (m1!=0) signal(pad,time)+=qweight; + total(pad,time)+=qweight; + if (indexRange[0]>pad) indexRange[0]=pad; + if (indexRange[1]time) indexRange[2]=time; + if (indexRange[3]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 + } + } } // 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, - Float_t **pList) +void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m, + Int_t *indexRange, Float_t **pList) { //---------------------------------------------------------------------- // Updates the list of tracks contributing to digits for a given row @@ -1573,84 +1657,81 @@ void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange, // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl //----------------------------------------------------------------- - TMatrix &signal = *m; + TMatrixF &signal = *m; // lop over nonzero digits - for(Int_t it=IndexRange[2];it 0. - // - - if(signal(ip,it)>0.){ - - pList[GlobalIndex] = new Float_t [6]; - - // set list to -1 + for(Int_t it=indexRange[2];ithighest){ - *(pList[GlobalIndex]+5) = middle; - *(pList[GlobalIndex]+4) = highest; - *(pList[GlobalIndex]+3) = signal(ip,it); + // check the signal magnitude - *(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); + 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 + // + + if(signal(ip,it)highest){ + *(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; + } + else if (signal(ip,it)>middle){ + *(pList[globalIndex]+5) = middle; + *(pList[globalIndex]+4) = signal(ip,it); + + *(pList[globalIndex]+2) = *(pList[globalIndex]+1); + *(pList[globalIndex]+1) = label; + } + else{ + *(pList[globalIndex]+5) = signal(ip,it); + *(pList[globalIndex]+2) = label; + } + } + } // end of loop over pads } // end of loop over time bins - - - }//end of GetList //___________________________________________________________________ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, @@ -1667,302 +1748,202 @@ 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[5]; AliTPChit *tpcHit; // pointer to a sigle TPC hit + //MI change + TBranch * branch=0; + if (fHitType>1) branch = TH->GetBranch("TPC2"); + else branch = TH->GetBranch("TPC"); + //---------------------------------------------- // Create TObjArray-s, one for each row, // each TObjArray will store the TVectors - // of electrons, one TVector per each track. + // of electrons, one TVectors per each track. //---------------------------------------------- - for(i=0; iGetEntry(track); // get next track + + //M.I. changes - TH->GetEvent(track); // get next track - Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track - - 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 - - currentTrack = tpcHit->fTrack; // track number - if(currentTrack != previousTrack){ - - // store already filled fTrack - - for(i=0;i0){ - TVector &v = *tr[i]; - v(0) = previousTrack; - tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary - row[i]->Add(tr[i]); - } - else{ - delete tr[i]; // delete empty TVector - tr[i]=0; - } - } - - n_of_electrons[i]=0; - tr[i] = new TVector(361); // 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) - - //--------------------------------------------------- - // Calculate the electron attachment probability - //--------------------------------------------------- - - Float_t time = 1.e6*(z_end-TMath::Abs(tpcHit->fZ))/fTPCParam->GetDriftV(); - // in microseconds! - 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]++; - //---------------------------------- - // 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 - } - } - - TVector &v = *tr[row_number]; - Int_t idx = 3*n_of_electrons[row_number]-2; - - v(idx)= xyz[0]; // X - v(idx+1)=xyz[1]; // Y (along the pad-row) - v(idx+2)=xyz[2]; // Z - - } // end of loop over electrons - - } // end of loop over hits - } // end of loop over tracks - - // - // store remaining track (the last one) if not empty - // - - for(i=0;i0){ - TVector &v = *tr[i]; - v(0) = previousTrack; - tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary - row[i]->Add(tr[i]); - } - else{ - delete tr[i]; - tr[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()); - - //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 - + Int_t sector=tpcHit->fSector; // sector number + if(sector != isec){ + tpcHit = (AliTPChit*) NextHit(); + continue; + } - 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 + // Remove hits which arrive before the TPC opening gate signal + if(((fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z())) + /fTPCParam->GetDriftV()+tpcHit->Time())GetGateDelay()) { + tpcHit = (AliTPChit*) NextHit(); + continue; + } + currentTrack = tpcHit->Track(); // track number - // loop over pads, from pmin to pmax + if(currentTrack != previousTrack){ + + // store already filled fTrack + + for(i=0;i0){ + TVector &v = *tracks[i]; + v(0) = previousTrack; + tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary + row[i]->Add(tracks[i]); + } + else { + delete tracks[i]; // delete empty TVector + tracks[i]=0; + } + } - for(Int_t i3=pmin;i3 nPadsSignal-nPadsDiff) continue; + } // end of loop over rows + + previousTrack=currentTrack; // update track label + } + + Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons) - TruePad -= nPadsDiff; - signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge! + //--------------------------------------------------- + // Calculate the electron attachment probability + //--------------------------------------------------- - } // end of loop over pads - } // end of loop over time bins - } // end of loop over electrons + Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z())) + /fTPCParam->GetDriftV(); + // in microseconds! + Float_t attProb = fTPCParam->GetAttCoef()* + fTPCParam->GetOxyCont()*time; // fraction! + + //----------------------------------------------- + // Loop over electrons + //----------------------------------------------- + Int_t index[3]; + index[1]=isec; + for(Int_t nel=0;nelRndm(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); + Int_t rowNumber; + fTPCParam->GetPadRow(xyz,index); + // Electron track time (for pileup simulation) + xyz[4] = tpcHit->Time()/fTPCParam->GetTSample(); + // row 0 - cross talk from the innermost row + // row fNRow+1 cross talk from the outermost row + rowNumber = index[2]+1; + //transform position to local digit coordinates + //relative to nearest pad row + if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue; + Float_t x1,y1; + if (isec GetNInnerSector()) { + x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth(); + y1 = fTPCParam->GetYInner(rowNumber); + } + else{ + x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth(); + y1 = fTPCParam->GetYOuter(rowNumber); + } + // gain inefficiency at the wires edges - linear + x1=TMath::Abs(x1); + y1-=1.; + if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); + + nofElectrons[rowNumber]++; + //---------------------------------- + // Expand vector if necessary + //---------------------------------- + if(nofElectrons[rowNumber]>120){ + Int_t range = tracks[rowNumber]->GetNrows(); + if((nofElectrons[rowNumber])>(range-1)/5){ + + tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons + } + } + + TVector &v = *tracks[rowNumber]; + Int_t idx = 5*nofElectrons[rowNumber]-4; + Real_t * position = &(((TVector&)v)(idx)); //make code faster + memcpy(position,xyz,5*sizeof(Float_t)); + + } // end of loop over electrons - } // end of loop over tracks + tpcHit = (AliTPChit*)NextHit(); + + } // end of loop over hits + } // end of loop over tracks -} // end of GetCrossTalk + // + // store remaining track (the last one) if not empty + // + + for(i=0;i0){ + TVector &v = *tracks[i]; + v(0) = previousTrack; + tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary + row[i]->Add(tracks[i]); + } + else{ + delete tracks[i]; + tracks[i]=0; + } + } + + delete [] tracks; + delete [] nofElectrons; +} // end of MakeSector //_____________________________________________________________________________ @@ -1971,16 +1952,7 @@ void AliTPC::Init() // // Initialise TPC detector after definition of geometry // - 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"); + AliDebug(1,"*********************************************"); } //_____________________________________________________________________________ @@ -1989,25 +1961,24 @@ void AliTPC::MakeBranch(Option_t* option) // // Create Tree branches for the TPC. // + AliDebug(1,""); Int_t buffersize = 4000; char branchname[10]; sprintf(branchname,"%s",GetName()); - + + const char *h = strstr(option,"H"); + + if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03 + AliDetector::MakeBranch(option); - - char *D = strstr(option,"D"); - - if (fDigits && gAlice->TreeD() && D) { - gAlice->TreeD()->Branch(branchname,&fDigits, buffersize); - printf("Making Branch %s for digits\n",branchname); + + const char *d = strstr(option,"D"); + + if (fDigits && fLoader->TreeD() && d) { + MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0); } - 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>1) MakeBranch2(option,0); // MI change 14.09.2000 } //_____________________________________________________________________________ @@ -2015,88 +1986,12 @@ 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(); -} - -//_____________________________________________________________________________ -void AliTPC::SetSecAL(Int_t sec) -{ - //--------------------------------------------------- - // Activate/deactivate selection for lower sectors - //--------------------------------------------------- - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - - fSecAL = sec; -} - -//_____________________________________________________________________________ -void AliTPC::SetSecAU(Int_t sec) -{ - //---------------------------------------------------- - // Activate/deactivate selection for upper sectors - //--------------------------------------------------- - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - - fSecAU = sec; + if (fDigits) fDigits->Clear(); } -//_____________________________________________________________________________ -void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6) -{ - //---------------------------------------- - // Select active lower sectors - //---------------------------------------- - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - - fSecLows[0] = s1; - fSecLows[1] = s2; - fSecLows[2] = s3; - fSecLows[3] = s4; - fSecLows[4] = s5; - fSecLows[5] = s6; -} -//_____________________________________________________________________________ -void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6, - Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10, - Int_t s11 , Int_t s12) -{ - //-------------------------------- - // Select active upper sectors - //-------------------------------- - - //----------------------------------------------------------------- - // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl - //----------------------------------------------------------------- - - fSecUps[0] = s1; - fSecUps[1] = s2; - fSecUps[2] = s3; - fSecUps[3] = s4; - fSecUps[4] = s5; - fSecUps[5] = s6; - fSecUps[6] = s7; - fSecUps[7] = s8; - fSecUps[8] = s9; - fSecUps[9] = s10; - fSecUps[10] = s11; - fSecUps[11] = s12; -} //_____________________________________________________________________________ void AliTPC::SetSens(Int_t sens) @@ -2114,99 +2009,48 @@ 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 AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const -{ - // - // Transformation from local to global coordinate system - // - 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) -{ - // - // 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 +void AliTPC::TransportElectron(Float_t *xyz, Int_t *index) { // - //make AliTPCcluster sortabale - return kTRUE; -} - - - -ClassImp(AliTPCdigit) - -//_____________________________________________________________________________ -AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits): - AliDigit(tracks) -{ + // electron transport taking into account: + // 1. diffusion, + // 2.ExB at the wires + // 3. nonisochronity // - // Creates a TPC digit object + // xyz and index must be already transformed to system 1 // - fSector = digits[0]; - fPadRow = digits[1]; - fPad = digits[2]; - fTime = digits[3]; - fSignal = digits[4]; -} - + 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 dx = fTPCParam->Transform2to2NearestWire(xyz,index); + xyz[1]+=dx*(fTPCParam->GetOmegaTau()); + } + //add nonisochronity (not implemented yet) +} + ClassImp(AliTPChit) //_____________________________________________________________________________ @@ -2222,391 +2066,484 @@ AliHit(shunt,track) fY = hits[1]; fZ = hits[2]; fQ = hits[3]; + fTime = hits[4]; } - -ClassImp(AliTPCtrack) - -//_____________________________________________________________________________ -AliTPCtrack::AliTPCtrack(Float_t *hits) +//________________________________________________________________________ +// Additional code because of the AliTPCTrackHitsV2 + +void AliTPC::MakeBranch2(Option_t *option,const char */*file*/) { // - // Default creator for a TPC reconstructed track object + // Create a new branch in the current Root Tree + // The branch of fHits is automatically split + // MI change 14.09.2000 + AliDebug(1,""); + if (fHitType<2) return; + char branchname[10]; + sprintf(branchname,"%s2",GetName()); + // + // Get the pointer to the header + const char *cH = strstr(option,"H"); // - ref=hits[0]; // This is dummy code ! + if (fTrackHits && TreeH() && cH && fHitType&4) { + AliDebug(1,"Making branch for Type 4 Hits"); + TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99); + } + +// if (fTrackHitsOld && TreeH() && cH && fHitType&2) { +// AliDebug(1,"Making branch for Type 2 Hits"); +// AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld, +// TreeH(),fBufferSize,99); +// TreeH()->GetListOfBranches()->Add(branch); +// } } -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 (sec<24) { - fAlpha=(sec%12)*alpha_low; - } else { - fAlpha=(sec%24)*alpha_up; + //Sets tree address for hits + if (fHitType<=1) { + if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03 + AliDetector::SetTreeAddress(); } - clusters.AddLast((AliTPCcluster*)(&c)); + if (fHitType>1) SetTreeAddress2(); } -//_____________________________________________________________________________ -AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C), - clusters(t.clusters.GetEntriesFast()) +void AliTPC::SetTreeAddress2() { // - // 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; iGetBranch(branchname); + if (branch) { + branch->SetAddress(&fTrackHits); + AliDebug(1,"fHitType&4 Setting"); + } + else + AliDebug(1,"fHitType&4 Failed (can not find branch)"); + + } + // if ((treeH)&&(fHitType&2)) { +// branch = treeH->GetBranch(branchname); +// if (branch) { +// branch->SetAddress(&fTrackHitsOld); +// AliDebug(1,"fHitType&2 Setting"); +// } +// else +// AliDebug(1,"fHitType&2 Failed (can not find branch)"); +// } + //set address to TREETR + + TTree *treeTR = TreeTR(); + if (treeTR && fTrackReferences) { + branch = treeTR->GetBranch(GetName()); + if (branch) branch->SetAddress(&fTrackReferences); + } + } -//_____________________________________________________________________________ -Double_t AliTPCtrack::GetY(Double_t xk) const +void AliTPC::FinishPrimary() { + if (fTrackHits &&fHitType&4) fTrackHits->FlushHitStack(); + // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->FlushHitStack(); +} + + +void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits) +{ // - // - // - 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); + // add hit to the list + Int_t rtrack; + if (fIshunt) { + int primary = gAlice->GetMCApp()->GetPrimary(track); + gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit); + rtrack=primary; + } else { + rtrack=track; + gAlice->GetMCApp()->FlagTrack(track); + } + if (fTrackHits && fHitType&4) + fTrackHits->AddHitKartez(vol[0],rtrack, hits[0], + hits[1],hits[2],(Int_t)hits[3],hits[4]); + // if (fTrackHitsOld &&fHitType&2 ) +// fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0], +// hits[1],hits[2],(Int_t)hits[3]); + } -//_____________________________________________________________________________ -int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) +void AliTPC::ResetHits() +{ + if (fHitType&1) AliDetector::ResetHits(); + if (fHitType>1) ResetHits2(); +} + +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 && fHitType&4) fTrackHits->Clear(); + // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->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; -} +} -//_____________________________________________________________________________ -void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) +AliHit* AliTPC::FirstHit(Int_t track) +{ + if (fHitType>1) return FirstHit2(track); + return AliDetector::FirstHit(track); +} +AliHit* AliTPC::NextHit() { // - // Propagate a reconstructed track from the vertex + // gets next hit // - 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>1) 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(); + 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 && fHitType&4) { + fTrackHits->First(); + return fTrackHits->GetHit(); } - - TMatrix saveC=C; - C.Mult(K,tmp); C-=saveC; C*=-1; - - clusters.AddLast((AliTPCcluster*)c); - chi2 += chisq; + // if (fTrackHitsOld && fHitType&2) { +// fTrackHitsOld->First(); +// return fTrackHitsOld->GetHit(); +// } + + 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 (fTrackHitsOld && fHitType&2) { +// fTrackHitsOld->Next(); +// return fTrackHitsOld->GetHit(); +// } + 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_t a = 0; + + if(fHitType==1) AliDetector::LoadPoints(a); + else LoadPoints2(a); +} + + +void AliTPC::RemapTrackHitIDs(Int_t *map) +{ // - // - int num_of_clusters=clusters.GetEntriesFast(); - for (int i=0; iUse(); + // remapping + // + if (!fTrackHits) return; + +// if (fTrackHitsOld && fHitType&2){ +// AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo; +// for (UInt_t i=0;iGetSize();i++){ +// AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i)); +// info->fTrackID = map[info->fTrackID]; +// } +// } +// if (fTrackHitsOld && fHitType&4){ + if (fTrackHits && fHitType&4){ + TClonesArray * arr = fTrackHits->GetArray();; + for (Int_t i=0;iGetEntriesFast();i++){ + AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i)); + info->SetTrackID(map[info->GetTrackID()]); + } } } -//_____________________________________________________________________________ -Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const +Bool_t AliTPC::TrackInVolume(Int_t id,Int_t track) { - // - // 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; + //return bool information - is track in given volume + //load only part of the track information + //return true if current track is in volume + // + // return kTRUE; + // if (fTrackHitsOld && fHitType&2) { +// TBranch * br = TreeH()->GetBranch("fTrackHitsInfo"); +// br->GetEvent(track); +// AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo; +// for (UInt_t j=0;jGetSize();j++){ +// if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE; +// } +// } + + if (fTrackHits && fHitType&4) { + TBranch * br1 = TreeH()->GetBranch("fVolumes"); + TBranch * br2 = TreeH()->GetBranch("fNVolumes"); + br2->GetEvent(track); + br1->GetEvent(track); + Int_t *volumes = fTrackHits->GetVolumes(); + Int_t nvolumes = fTrackHits->GetNVolumes(); + if (!volumes && nvolumes>0) { + AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes)); + return kFALSE; + } + for (Int_t j=0;jGetBranch("fSector"); + br->GetEvent(track); + for (Int_t j=0;jGetEntriesFast();j++){ + if ( ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE; + } + } + return kFALSE; + } //_____________________________________________________________________________ -int AliTPCtrack::GetLab() const +void AliTPC::LoadPoints2(Int_t) { // + // Store x, y, z of all hits in memory // + // if (fTrackHits == 0 && fTrackHitsOld==0) return; + 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; jGetEntriesFast(); + // if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast(); - int max=0; - for (i=0; imax) {max=s[i].max; lab=s[i].lab;} - - for (i=0; ifTracks[1]) == lab || - TMath::Abs(c->fTracks[2]) == lab ) max++; + if (nhits == 0) return; + Int_t tracks = gAlice->GetMCApp()->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); + while (ahit){ + 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->GetMCApp()->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;iGetTrack(); + 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],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 + + +AliLoader* AliTPC::MakeLoader(const char* topfoldername) { - // - // - // - 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; + //Makes TPC loader + fLoader = new AliTPCLoader(GetName(),topfoldername); + return fLoader; +} + +//////////////////////////////////////////////////////////////////////// +AliTPCParam* AliTPC::LoadTPCParam(TFile *file) { +// +// load TPC paarmeters from a given file or create new if the object +// is not found there +// 12/05/2003 This method should be moved to the AliTPCLoader +// and one has to decide where to store the TPC parameters +// M.Kowalski + char paramName[50]; + sprintf(paramName,"75x40_100x60_150x60"); + AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName); + if (paramTPC) { + AliDebugClass(1,Form("TPC parameters %s found.",paramName)); + } else { + AliWarningClass("TPC parameters not found. Create new (they may be incorrect)"); + paramTPC = new AliTPCParamSR; } - return m; + return paramTPC; + +// the older version of parameters can be accessed with this code. +// In some cases, we have old parameters saved in the file but +// digits were created with new parameters, it can be distinguish +// by the name of TPC TreeD. The code here is just for the case +// we would need to compare with old data, uncomment it if needed. +// +// char paramName[50]; +// sprintf(paramName,"75x40_100x60"); +// AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName); +// if (paramTPC) { +// cout<<"TPC parameters "<Get(paramName); +// if (paramTPC) { +// cout<<"TPC parameters "<