X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPC.cxx;h=6020cde38670a2e8879f30244c0f54c3ef5eb149;hb=a149b0816c70e616ddecae4d80f57db2934608a0;hp=326a5e3ac8a280b6818073f10b63a8c1a1af2ed0;hpb=9fe90376639e73d64a7c07e508d68f5980324ba6;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPC.cxx b/TPC/AliTPC.cxx index 326a5e3ac8a..6020cde3867 100644 --- a/TPC/AliTPC.cxx +++ b/TPC/AliTPC.cxx @@ -13,127 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.45 2001/11/03 13:33:48 kowal2 -Updated algorithms in Hits2SDigits, SDigits2Digits, -Hits2ExactClusters. -Added method Merge - -Revision 1.44 2001/08/30 09:28:48 hristov -TTree names are explicitly set via SetName(name) and then Write() is called - -Revision 1.43 2001/07/28 12:02:54 hristov -Branch split level set to 99 - -Revision 1.42 2001/07/28 11:38:52 hristov -Loop variable declared once - -Revision 1.41 2001/07/28 10:53:50 hristov -Digitisation done according to the general scheme (M.Ivanov) - -Revision 1.40 2001/07/27 13:03:14 hristov -Default Branch split level set to 99 - -Revision 1.39 2001/07/26 09:09:34 kowal2 -Hits2Reco method added - -Revision 1.38 2001/07/20 14:32:43 kowal2 -Processing of many events possible now - -Revision 1.37 2001/06/12 07:17:18 kowal2 -Hits2SDigits method implemented (summable digits) - -Revision 1.36 2001/05/16 14:57:25 alibrary -New files for folders and Stack - -Revision 1.35 2001/05/08 16:02:22 kowal2 -Updated material specifications - -Revision 1.34 2001/05/08 15:00:15 hristov -Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov) - -Revision 1.33 2001/04/03 12:40:43 kowal2 -Removed printouts - -Revision 1.32 2001/03/12 17:47:36 hristov -Changes needed on Sun with CC 5.0 - -Revision 1.31 2001/03/12 08:21:50 kowal2 -Corrected C++ bug in the material definitions - -Revision 1.30 2001/03/01 17:34:47 kowal2 -Correction due to the accuracy problem - -Revision 1.29 2001/02/28 16:34:40 kowal2 -Protection against nonphysical values of the avalanche size, -10**6 is the maximum - -Revision 1.28 2001/01/26 19:57:19 hristov -Major upgrade of AliRoot code - -Revision 1.27 2001/01/13 17:29:33 kowal2 -Sun compiler correction - -Revision 1.26 2001/01/10 07:59:43 kowal2 -Corrections to load points from the noncompressed hits. - -Revision 1.25 2000/11/02 07:25:31 kowal2 -Changes due to the new hit structure. -Memory leak removed. - -Revision 1.24 2000/10/05 16:06:09 kowal2 -Forward declarations. Changes due to a new class AliComplexCluster. - -Revision 1.23 2000/10/02 21:28:18 fca -Removal of useless dependecies via forward declarations - -Revision 1.22 2000/07/10 20:57:39 hristov -Update of TPC code and macros by M.Kowalski - -Revision 1.19.2.4 2000/06/26 07:39:42 kowal2 -Changes to obey the coding rules - -Revision 1.19.2.3 2000/06/25 08:38:41 kowal2 -Splitted from AliTPCtracking - -Revision 1.19.2.2 2000/06/16 12:59:28 kowal2 -Changed parameter settings - -Revision 1.19.2.1 2000/06/09 07:15:07 kowal2 - -Defaults loaded automatically (hard-wired) -Optional parameters can be set via macro called in the constructor - -Revision 1.19 2000/04/18 19:00:59 fca -Small bug fixes to TPC files - -Revision 1.18 2000/04/17 09:37:33 kowal2 -removed obsolete AliTPCDigitsDisplay.C - -Revision 1.17.2.2 2000/04/10 08:15:12 kowal2 - -New, experimental data structure from M. Ivanov -New tracking algorithm -Different pad geometry for different sectors -Digitization rewritten - -Revision 1.17.2.1 2000/04/10 07:56:53 kowal2 -Not used anymore - removed - -Revision 1.17 2000/01/19 17:17:30 fca -Introducing a list of lists of hits -- more hits allowed for detector now - -Revision 1.16 1999/11/05 09:29:23 fca -Accept only signals > 0 - -Revision 1.15 1999/10/08 06:26:53 fca -Removed ClustersIndex - not used anymore - -Revision 1.14 1999/09/29 09:24:33 fca -Introduction of the Copyright and cvs Log - -*/ +/* $Id$ */ /////////////////////////////////////////////////////////////////////////////// // // @@ -153,54 +33,72 @@ Introduction of the Copyright and cvs Log // +#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 "AliRun.h" -#include -#include -#include -#include "AliMC.h" -#include "AliMagF.h" - +#include +#include +#include +#include +#include +#include -#include "AliTPCParamSR.h" -#include "AliTPCPRF2D.h" -#include "AliTPCRF1D.h" +#include "AliArrayBranch.h" #include "AliDigits.h" -#include "AliSimDigits.h" -#include "AliTPCTrackHits.h" +#include "AliMagF.h" #include "AliPoints.h" -#include "AliArrayBranch.h" - - +#include "AliRun.h" +#include "AliRunLoader.h" +#include "AliSimDigits.h" +#include "AliTPC.h" +#include "AliTPC.h" #include "AliTPCDigitsArray.h" -#include "AliComplexCluster.h" -#include "AliClusters.h" -#include "AliTPCClustersRow.h" -#include "AliTPCClustersArray.h" - -#include "AliTPCcluster.h" -#include "AliTPCclusterer.h" -#include "AliTPCtracker.h" - -#include -#include - +#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" ClassImp(AliTPC) +//_____________________________________________________________________________ +// helper class for fast matrix and vector manipulation - no range checking +// origin - Marian Ivanov + +class AliTPCFastMatrix : public TMatrix { +public : + AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb); +#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1) + Float_t & UncheckedAt(Int_t rown, Int_t coln) const {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces + Float_t UncheckedAtFast(Int_t rown, Int_t coln) const {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces +#else + Float_t & UncheckedAt(Int_t rown, Int_t coln) const {return (fIndex[coln])[rown];} //fast acces + Float_t UncheckedAtFast(Int_t rown, Int_t coln) const {return (fIndex[coln])[rown];} //fast acces +#endif +}; + +AliTPCFastMatrix::AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb): + TMatrix(rowlwb, rowupb,collwb,colupb) + { + }; //_____________________________________________________________________________ AliTPC::AliTPC() { @@ -211,13 +109,15 @@ AliTPC::AliTPC() fHits = 0; fDigits = 0; fNsectors = 0; - //MI changes fDigitsArray = 0; - fClustersArray = 0; fDefaults = 0; - fTrackHits = 0; - fHitType = 2; - fTPCParam = 0; + fTrackHits = 0; + fTrackHitsOld = 0; + fHitType = 2; //default CONTAINERS - based on ROOT structure + fTPCParam = 0; + fNoiseTable = 0; + fActiveSectors =0; + } //_____________________________________________________________________________ @@ -231,18 +131,24 @@ AliTPC::AliTPC(const char *name, const char *title) // // Initialise arrays of hits and digits fHits = new TClonesArray("AliTPChit", 176); - gAlice->AddHitList(fHits); - //MI change + gAlice->GetMCApp()->AddHitList(fHits); fDigitsArray = 0; - fClustersArray= 0; fDefaults = 0; // - fTrackHits = new AliTPCTrackHits; //MI - 13.09.2000 + fTrackHits = new AliTPCTrackHitsV2; fTrackHits->SetHitPrecision(0.002); fTrackHits->SetStepPrecision(0.003); - fTrackHits->SetMaxDistance(100); + fTrackHits->SetMaxDistance(100); + + fTrackHitsOld = new AliTPCTrackHits; //MI - 13.09.2000 + fTrackHitsOld->SetHitPrecision(0.002); + fTrackHitsOld->SetStepPrecision(0.003); + fTrackHitsOld->SetMaxDistance(100); + + fNoiseTable =0; fHitType = 2; + fActiveSectors = 0; // // Initialise counters fNsectors = 0; @@ -268,6 +174,11 @@ AliTPC::AliTPC(const char *name, const char *title) } //_____________________________________________________________________________ +AliTPC::AliTPC(const AliTPC& t):AliDetector(t){ + // + // dummy copy constructor + // +} AliTPC::~AliTPC() { // @@ -279,6 +190,9 @@ AliTPC::~AliTPC() delete fDigits; delete fTPCParam; delete fTrackHits; //MI 15.09.2000 + delete fTrackHitsOld; //MI 10.12.2001 + if (fNoiseTable) delete [] fNoiseTable; + } //_____________________________________________________________________________ @@ -293,10 +207,10 @@ void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits) TClonesArray &lhits = *fHits; new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits); } - if (fHitType&2) + if (fHitType>1) AddHit2(track,vol,hits); } - + //_____________________________________________________________________________ void AliTPC::BuildGeometry() { @@ -383,7 +297,7 @@ void AliTPC::BuildGeometry() } //_____________________________________________________________________________ -Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) +Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) const { // // Calculate distance from TPC to mouse on the display @@ -392,13 +306,6 @@ Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) return 9999; } -void AliTPC::Clusters2Tracks(TFile *of) { - //----------------------------------------------------------------- - // This is a track finder. - //----------------------------------------------------------------- - AliTPCtracker tracker(fTPCParam); - tracker.Clusters2Tracks(gFile,of); -} //_____________________________________________________________________________ void AliTPC::CreateMaterials() @@ -512,11 +419,11 @@ void AliTPC::CreateMaterials() // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds //---------------------------------------------------------------- - char namate[21]; + char namate[21]=""; density = 0.; Float_t am=0; Int_t nc; - Float_t rho,absl,X0,buf[1]; + Float_t rho,absl,x0,buf[1]; Int_t nbuf; Float_t a,z; @@ -525,7 +432,7 @@ void AliTPC::CreateMaterials() // retrive material constants - gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf); + gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf); amat[nc] = a; zmat[nc] = z; @@ -755,7 +662,7 @@ void AliTPC::CreateMaterials() // get epoxy - gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf); + gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,x0,absl,buf,nbuf); // Carbon fiber @@ -768,7 +675,7 @@ void AliTPC::CreateMaterials() // get SiO2 - gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf); + gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,x0,absl,buf,nbuf); wmat[0]=0.725; // by weight! wmat[1]=0.275; @@ -808,425 +715,249 @@ void AliTPC::CreateMaterials() } - -void AliTPC::Digits2Clusters(TFile *of, Int_t eventnumber) +void AliTPC::GenerNoise(Int_t tablesize) { - //----------------------------------------------------------------- - // This is a simple cluster finder. - //----------------------------------------------------------------- - AliTPCclusterer::Digits2Clusters(fTPCParam,of,eventnumber); -} - -extern Double_t SigmaY2(Double_t, Double_t, Double_t); -extern Double_t SigmaZ2(Double_t, Double_t); -//_____________________________________________________________________________ -void AliTPC::Hits2Clusters(TFile *of, Int_t eventn) -{ - //-------------------------------------------------------- - // 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 - //----------------------------------------------------------------- - // Adopted to Marian's cluster data structure by I.Belikov, CERN, - // Jouri.Belikov@cern.ch - //---------------------------------------------------------------- - - ///////////////////////////////////////////////////////////////////////////// // - //--------------------------------------------------------------------- - // ALICE TPC Cluster Parameters - //-------------------------------------------------------------------- - - - - // Cluster width in rphi - const Float_t kACrphi=0.18322; - const Float_t kBCrphi=0.59551e-3; - const Float_t kCCrphi=0.60952e-1; - // Cluster width in z - const Float_t kACz=0.19081; - const Float_t kBCz=0.55938e-3; - const Float_t kCCz=0.30428; - - TDirectory *savedir=gDirectory; - - if (!of->IsOpen()) { - cerr<<"AliTPC::Hits2Clusters(): output file not open !\n"; - return; - } - - //if(fDefaults == 0) SetDefaults(); - - Float_t sigmaRphi,sigmaZ,clRphi,clZ; + //Generate table with noise // - TParticle *particle; // pointer to a given particle - AliTPChit *tpcHit; // pointer to a sigle TPC hit - Int_t sector; - Int_t ipart; - Float_t xyz[5]; - Float_t pl,pt,tanth,rpad,ratio; - Float_t cph,sph; - - //--------------------------------------------------------------- - // Get the access to the tracks - //--------------------------------------------------------------- + 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 - TTree *tH = gAlice->TreeH(); - Stat_t ntracks = tH->GetEntries(); - - //Switch to the output file - of->cd(); + Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac(); + for (Int_t i=0;iGaus(0,norm); +} - char cname[100]; +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()); +} - sprintf(cname,"TreeC_TPC_%d",eventn); - fTPCParam->Write(fTPCParam->GetTitle()); - AliTPCClustersArray carray; - carray.Setup(fTPCParam); - carray.SetClusterType("AliTPCcluster"); - carray.MakeTree(); +Bool_t AliTPC::IsSectorActive(Int_t sec) const +{ + // + // check if the sector is active + if (!fActiveSectors) return kTRUE; + else return fActiveSectors[sec]; +} - Int_t nclusters=0; //cluster counter - - //------------------------------------------------------------ - // Loop over all sectors (72 sectors for 20 deg - // segmentation for both lower and upper sectors) - // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0 - // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0 - // - // First cluster for sector 0 starts at "0" - //------------------------------------------------------------ - - for(Int_t isec=0;isecGetNSector();isec++){ - //MI change - fTPCParam->AdjustCosSin(isec,cph,sph); - - //------------------------------------------------------------ - // Loop over tracks - //------------------------------------------------------------ +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; - for(Int_t track=0;trackGetEvent(track); - // - // Get number of the TPC hits - // - // nhits=fHits->GetEntriesFast(); - // - - tpcHit = (AliTPChit*)FirstHit(-1); - - // Loop over hits - // - // for(Int_t hit=0;hitUncheckedAt(hit); - - while(tpcHit){ - - if (tpcHit->fQ == 0.) { - tpcHit = (AliTPChit*) NextHit(); - continue; //information about track (I.Belikov) - } - sector=tpcHit->fSector; // sector number - - - // if(sector != isec) continue; //terminate iteration - - if(sector != isec){ - tpcHit = (AliTPChit*) NextHit(); - continue; - } - ipart=tpcHit->Track(); - particle=gAlice->Particle(ipart); - pl=particle->Pz(); - pt=particle->Pt(); - if(pt < 1.e-9) pt=1.e-9; - tanth=pl/pt; - tanth = TMath::Abs(tanth); - rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y()); - ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason - - // space-point resolutions - - sigmaRphi=SigmaY2(rpad,tanth,pt); - sigmaZ =SigmaZ2(rpad,tanth ); - - // cluster widths - - clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio; - clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth; - - // temporary protection - - if(sigmaRphi < 0.) sigmaRphi=0.4e-3; - if(sigmaZ < 0.) sigmaZ=0.4e-3; - if(clRphi < 0.) clRphi=2.5e-3; - if(clZ < 0.) clZ=2.5e-5; - - // - - // - // 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->X()*cph + tpcHit->Y()*sph; - Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph; - xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y - Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ? - fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle(); - Float_t ymax=xprim*TMath::Tan(0.5*alpha); - if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim; - xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z - if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z(); - xyz[2]=sigmaRphi; // fSigmaY2 - xyz[3]=sigmaZ; // fSigmaZ2 - xyz[4]=tpcHit->fQ; // q - - AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow); - if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow); - - Int_t tracks[3]={tpcHit->Track(), -1, -1}; - AliTPCcluster cluster(tracks,xyz); - - clrow->InsertCluster(&cluster); nclusters++; - - tpcHit = (AliTPChit*)NextHit(); - - - } // end of loop over hits - - } // end of loop over tracks +} - Int_t nrows=fTPCParam->GetNRow(isec); - for (Int_t irow=0; irowGetNSector()]; + 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) + { + Fatal("SetActiveSectors","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 + if (GetDebug()) cout<<"\nAliTPC::SetActiveSectors(): Got "<GetBranch("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; + } + } + } +} - } // end of loop over sectors - - cerr<<"Number of made clusters : "<SetName(cname); - carray.GetTree()->Write(); - savedir->cd(); //switch back to the input file - -} // end of function -//_________________________________________________________________ -void AliTPC::Hits2ExactClustersSector(Int_t isec) +//_____________________________________________________________________________ +void AliTPC::Digits2Raw() { - //-------------------------------------------------------- - //calculate exact cross point of track and given pad row - //resulting values are expressed in "digit" coordinata - //-------------------------------------------------------- +// convert digits of the current event to raw data - //----------------------------------------------------------------- - // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de - //----------------------------------------------------------------- - // - if (fClustersArray==0){ + static const Int_t kThreshold = 0; + static const Bool_t kCompress = kTRUE; + + fLoader->LoadDigits(); + TTree* digits = fLoader->TreeD(); + if (!digits) { + Error("Digits2Raw", "no digits tree"); return; } - // - TParticle *particle; // pointer to a given particle - AliTPChit *tpcHit; // pointer to a sigle TPC hit - // Int_t sector,nhits; - Int_t ipart; - const Int_t kcmaxhits=30000; - TVector * xxxx = new TVector(kcmaxhits*4); - TVector & xxx = *xxxx; - Int_t maxhits = kcmaxhits; - //construct array for each padrow - for (Int_t i=0; iGetNRow(isec);i++) - fClustersArray->CreateRow(isec,i); - - //--------------------------------------------------------------- - // Get the access to the tracks - //--------------------------------------------------------------- - - TTree *tH = gAlice->TreeH(); - Stat_t ntracks = tH->GetEntries(); - Int_t npart = gAlice->GetNtrack(); - TBranch * br = tH->GetBranch("fTrackHitsInfo"); - //MI change - TBranch * branch=0; - if (fHitType&2) branch = tH->GetBranch("TPC2"); - else branch = tH->GetBranch("TPC"); - //------------------------------------------------------------ - // Loop over tracks - //------------------------------------------------------------ - - for(Int_t track=0;trackGetEvent(track); - AliObjectArray * ar = fTrackHits->fTrackHitsInfo; - for (UInt_t j=0;jGetSize();j++){ - if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE; - } + 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 (!isInSector) continue; - //MI change - branch->GetEntry(track); // get next track + 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 - // - // Get number of the TPC hits and a pointer - // to the particles - // - //nhits=fHits->GetEntriesFast(); - // - // Loop over hits - // - Int_t currentIndex=0; - Int_t lastrow=-1; //last writen row + delete buffer; + fLoader->UnloadDigits(); - //M.I. changes + AliTPCDDLRawData rawWriter; + rawWriter.SetVerbose(0); - tpcHit = (AliTPChit*)FirstHit(-1); - while(tpcHit){ - - Int_t sector=tpcHit->fSector; // sector number - if(sector != isec){ - tpcHit = (AliTPChit*) NextHit(); - continue; - } + rawWriter.RawData(fileName); + gSystem->Unlink(fileName); - //for(Int_t hit=0;hitUncheckedAt(hit); - //if (tpcHit==0) continue; - //sector=tpcHit->fSector; // sector number - //if(sector != isec) continue; - ipart=tpcHit->Track(); - if (ipartParticle(ipart); - - //find row number - - Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()}; - Int_t index[3]={1,isec,0}; - Int_t currentrow = fTPCParam->GetPadRow(x,index) ; - if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;} - if (lastrow<0) lastrow=currentrow; - if (currentrow==lastrow){ - if ( currentIndex>=maxhits){ - maxhits+=kcmaxhits; - xxx.ResizeTo(4*maxhits); - } - xxx(currentIndex*4)=x[0]; - xxx(currentIndex*4+1)=x[1]; - xxx(currentIndex*4+2)=x[2]; - xxx(currentIndex*4+3)=tpcHit->fQ; - currentIndex++; - } - else - if (currentIndex>2){ - Float_t sumx=0; - Float_t sumx2=0; - Float_t sumx3=0; - Float_t sumx4=0; - Float_t sumy=0; - Float_t sumxy=0; - Float_t sumx2y=0; - Float_t sumz=0; - Float_t sumxz=0; - Float_t sumx2z=0; - Float_t sumq=0; - for (Int_t index=0;indexGetNPads(isec,lastrow)-1)/2; - Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+ - sumx2*(sumx*sumx3-sumx2*sumx2); - - Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+ - sumx2*(sumxy*sumx3-sumx2y*sumx2); - Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+ - sumx2*(sumxz*sumx3-sumx2z*sumx2); - - Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+ - sumx2*(sumx*sumx2y-sumx2*sumxy); - Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+ - sumx2*(sumx*sumx2z-sumx2*sumxz); - - Float_t y=detay/det+centralPad; - Float_t z=detaz/det; - Float_t by=detby/det; //y angle - Float_t bz=detbz/det; //z angle - sumy/=Float_t(currentIndex); - sumz/=Float_t(currentIndex); - - AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow)); - if (row!=0) { - AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ; - // AliComplexCluster cl; - cl->fX=z; - cl->fY=y; - cl->fQ=sumq; - cl->fSigmaX2=bz; - cl->fSigmaY2=by; - cl->fTracks[0]=ipart; - } - currentIndex=0; - lastrow=currentrow; - } //end of calculating cluster for given row - - - tpcHit = (AliTPChit*)NextHit(); - } // end of loop over hits - } // end of loop over tracks - //write padrows to tree - for (Int_t ii=0; iiGetNRow(isec);ii++) { - fClustersArray->StoreRow(isec,ii); - fClustersArray->ClearRow(isec,ii); + if (kCompress) { + Info("Digits2Raw", "compressing raw data"); + rawWriter.RawDataCompDecompress(kTRUE); + gSystem->Unlink("Statistics"); } - xxxx->Delete(); - } +//______________________________________________________________________ +AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const +{ + return new AliTPCDigitizer(manager); +} //__ -void AliTPC::SDigits2Digits2(Int_t eventnumber) +void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/) { //create digits from summable digits - - char sname[100]; - char dname[100]; - sprintf(sname,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber); - sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); + GenerNoise(500000); //create teble with noise //conect tree with sSDigits - TTree *t = (TTree *)gDirectory->Get(sname); + TTree *t = fLoader->TreeS(); + + if (t == 0x0) + { + fLoader->LoadSDigits("READ"); + t = fLoader->TreeS(); + if (t == 0x0) + { + Error("SDigits2Digits2","Can not get input TreeS"); + return; + } + } + + if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D"); + AliSimDigits digarr, *dummy=&digarr; - t->GetBranch("Segment")->SetAddress(&dummy); + TBranch* sdb = t->GetBranch("Segment"); + if (sdb == 0x0) + { + Error("SDigits2Digits2","Can not find branch with segments in TreeS."); + return; + } + + sdb->SetAddress(&dummy); + Stat_t nentries = t->GetEntries(); // set zero suppression @@ -1242,9 +973,9 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber) AliTPCDigitsArray *arr = new AliTPCDigitsArray; arr->SetClass("AliSimDigits"); arr->Setup(fTPCParam); - arr->MakeTree(fDigitsFile); + arr->MakeTree(fLoader->TreeD()); - AliTPCParam * par =fTPCParam; + AliTPCParam * par = fTPCParam; //Loop over segments of the TPC @@ -1255,6 +986,15 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber) cerr<<"AliTPC warning: invalid segment ID ! "<CreateRow(sec,row); Int_t nrows = digrow->GetNRows(); Int_t ncols = digrow->GetNCols(); @@ -1265,23 +1005,38 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber) digarr.ExpandTrackBuffer(); - for (Int_t rows=0;rowsGaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16); - q/=16; //conversion faktor - q=(Short_t)q; // digits are short integers - - if (q > zerosup){ - if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat()); - digrow->SetDigitFast(q,rows,col); - for (Int_t tr=0;tr<3;tr++) - -((AliSimDigits*)digrow)->SetTrackIDFast(digarr.GetTrackIDFast(rows,col,tr)-2,rows,col,tr); - } - } - } + 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(); + //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; + //if (ptracks0[0]==0) + // ptracks1[0]=1; + //else + 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); // cerr<GetTree()->SetName(dname); - arr->GetTree()->Write(); - delete arr; -} -//_________________________________________ -void AliTPC::Merge(TTree * intree, Int_t *mask, Int_t nin, Int_t outid) -{ - - //intree - pointer to array of trees with s digits - //mask - mask for each - //nin - number of inputs - //outid - event number of the output event - - //create digits from summable digits - //conect tree with sSDigits - - - AliSimDigits ** digarr = new AliSimDigits*[nin]; - for (Int_t i1=0;i1SetAddress(&digarr[i1]); - } - Stat_t nentries = intree[0].GetEntries(); - - //make tree with digits - char dname[100]; - sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),outid); - AliTPCDigitsArray *arr = new AliTPCDigitsArray; - arr->SetClass("AliSimDigits"); - arr->Setup(fTPCParam); - arr->MakeTree(fDigitsFile); - - // set zero suppression - - fTPCParam->SetZeroSup(2); - - // get zero suppression - - Int_t zerosup = fTPCParam->GetZeroSup(); - - - AliTPCParam * par =fTPCParam; - - //Loop over segments of the TPC - for (Int_t n=0; nExpandBuffer(); - digarr[i]->ExpandTrackBuffer(); - } - Int_t sec, row; - if (!par->AdjustSectorRow(digarr[0]->GetID(),sec,row)) { - cerr<<"AliTPC warning: invalid segment ID ! "<GetID()<CreateRow(sec,row); - Int_t nrows = digrow->GetNRows(); - Int_t ncols = digrow->GetNCols(); - - digrow->ExpandBuffer(); - digrow->ExpandTrackBuffer(); - - for (Int_t rows=0;rowsGetDigitFast(rows,col); - for (Int_t tr=0;tr<3;tr++) { - Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr); - if ( lab > 1) { - label[labptr]=lab+mask[i]-2; - labptr++; - } - } - } - //add noise - q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16); - - q/=16; //conversion faktor - q=(Short_t)q; - - if (q> zerosup){ - - if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat()); - digrow->SetDigitFast(q,rows,col); - for (Int_t tr=0;tr<3;tr++){ - if (trSetTrackIDFast(label[tr],rows,col,tr); - else - ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr); - } - } - } - } - arr->StoreRow(sec,row); - - arr->ClearRow(sec,row); - // cerr<GetTree()->SetName(dname); - arr->GetTree()->Write(); - - delete arr; - -} - -/*_________________________________________ -void AliTPC::SDigits2Digits(Int_t eventnumber) -{ - - - cerr<<"Digitizing TPC...\n"; - - Hits2Digits(eventnumber); + fLoader->WriteDigits("OVERWRITE"); - - //write results - - // char treeName[100]; - - // sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); - - // GetDigitsArray()->GetTree()->Write(treeName); + delete arr; } -*/ //__________________________________________________________________ void AliTPC::SetDefaults(){ - + // + // setting the defaults + // - cerr<<"Setting default parameters...\n"; + // cerr<<"Setting default parameters...\n"; // Set response functions + // + AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName()); + rl->CdGAFile(); AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60"); + if(param){ + printf("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...\n"); + delete param; + param = new AliTPCParamSR(); + } + else { + param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60"); + } + if(!param){ + printf("No TPC parameters found\n"); + exit(4); + } + + AliTPCPRF2D * prfinner = new AliTPCPRF2D; - AliTPCPRF2D * prfouter = new AliTPCPRF2D; + AliTPCPRF2D * prfouter1 = new AliTPCPRF2D; + AliTPCPRF2D * prfouter2 = new AliTPCPRF2D; AliTPCRF1D * rf = new AliTPCRF1D(kTRUE); rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.); rf->SetOffset(3*param->GetZSigma()); @@ -1446,13 +1090,30 @@ void AliTPC::SetDefaults(){ cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ; exit(3); } + + TString s; prfinner->Read("prf_07504_Gati_056068_d02"); - prfouter->Read("prf_10006_Gati_047051_d03"); + //PH Set different names + s=prfinner->GetGRF()->GetName(); + s+="in"; + prfinner->GetGRF()->SetName(s.Data()); + + prfouter1->Read("prf_10006_Gati_047051_d03"); + s=prfouter1->GetGRF()->GetName(); + s+="out1"; + prfouter1->GetGRF()->SetName(s.Data()); + + prfouter2->Read("prf_15006_Gati_047051_d03"); + s=prfouter2->GetGRF()->GetName(); + s+="out2"; + prfouter2->GetGRF()->SetName(s.Data()); + f->Close(); savedir->cd(); param->SetInnerPRF(prfinner); - param->SetOuterPRF(prfouter); + param->SetOuter1PRF(prfouter1); + param->SetOuter2PRF(prfouter2); param->SetTimeRF(rf); // set fTPCParam @@ -1464,14 +1125,54 @@ void AliTPC::SetDefaults(){ } //__________________________________________________________________ +void AliTPC::Hits2Digits() +{ + // + // creates digits from hits + // + + fLoader->LoadHits("read"); + fLoader->LoadDigits("recreate"); + AliRunLoader* runLoader = fLoader->GetRunLoader(); + + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + SetActiveSectors(); + Hits2Digits(iEvent); + } + + 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()) + { + Error("Hits2Digits","Can not load hits."); + } + } + SetTreeAddress(); + + if (fLoader->TreeD() == 0x0 ) + { + fLoader->MakeTree("D"); + if (fLoader->TreeD() == 0x0 ) + { + Error("Hits2Digits","Can not get TreeD"); + return; + } + } if(fDefaults == 0) SetDefaults(); // check if the parameters are set + GenerNoise(500000); //create teble with noise //setup TPCDigitsArray @@ -1480,29 +1181,36 @@ void AliTPC::Hits2Digits(Int_t eventnumber) AliTPCDigitsArray *arr = new AliTPCDigitsArray; arr->SetClass("AliSimDigits"); arr->Setup(fTPCParam); - arr->MakeTree(fDigitsFile); + + arr->MakeTree(fLoader->TreeD()); SetDigitsArray(arr); fDigitsSwitch=0; // standard digits - cerr<<"Digitizing TPC...\n"; + // cerr<<"Digitizing TPC -- normal digits...\n"; - for(Int_t isec=0;isecGetNSector();isec++) Hits2DigitsSector(isec); + for(Int_t isec=0;isecGetNSector();isec++) + if (IsSectorActive(isec)) + { + if (fDebug) Info("Hits2Digits","Sector %d is active.",isec); + Hits2DigitsSector(isec); + } + else + { + if (fDebug) Info("Hits2Digits","Sector %d is NOT active.",isec); + } - // write results - - char treeName[100]; - - sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); + 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); - GetDigitsArray()->GetTree()->SetName(treeName); - GetDigitsArray()->GetTree()->Write(); - - } - - //__________________________________________________________________ void AliTPC::Hits2SDigits2(Int_t eventnumber) { @@ -1514,21 +1222,43 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber) //---------------------------------------------------- // Loop over all sectors for a single event //---------------------------------------------------- +// AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName); + AliRunLoader* rl = fLoader->GetRunLoader(); + + rl->GetEvent(eventnumber); + if (fLoader->TreeH() == 0x0) + { + if(fLoader->LoadHits()) + { + Error("Hits2Digits","Can not load hits."); + return; + } + } + SetTreeAddress(); - if(fDefaults == 0) SetDefaults(); + if (fLoader->TreeS() == 0x0 ) + { + fLoader->MakeTree("S"); + } + + if(fDefaults == 0) SetDefaults(); + + GenerNoise(500000); //create table with noise //setup TPCDigitsArray if(GetDigitsArray()) delete GetDigitsArray(); + AliTPCDigitsArray *arr = new AliTPCDigitsArray; arr->SetClass("AliSimDigits"); arr->Setup(fTPCParam); - arr->MakeTree(fDigitsFile); + arr->MakeTree(fLoader->TreeS()); + SetDigitsArray(arr); - cerr<<"Digitizing TPC...\n"; + // cerr<<"Digitizing TPC -- summable digits...\n"; fDigitsSwitch=1; // summable digits @@ -1536,23 +1266,24 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber) fTPCParam->SetZeroSup(0); - for(Int_t isec=0;isecGetNSector();isec++) Hits2DigitsSector(isec); - - - // write results - - char treeName[100]; - - sprintf(treeName,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber); - - GetDigitsArray()->GetTree()->SetName(treeName); - GetDigitsArray()->GetTree()->Write(); - + for(Int_t isec=0;isecGetNSector();isec++) + if (IsSectorActive(isec)) + { +// cout<<"Sector "<WriteSDigits("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::Hits2SDigits() { @@ -1560,43 +1291,22 @@ void AliTPC::Hits2SDigits() // summable digits - 16 bit "ADC", no noise, no saturation //----------------------------------------------------------- - //---------------------------------------------------- - // Loop over all sectors for a single event - //---------------------------------------------------- - //MI change - for pp run - Int_t eventnumber = gAlice->GetEvNumber(); - - if(fDefaults == 0) SetDefaults(); - - //setup TPCDigitsArray + fLoader->LoadHits("read"); + fLoader->LoadSDigits("recreate"); + AliRunLoader* runLoader = fLoader->GetRunLoader(); - if(GetDigitsArray()) delete GetDigitsArray(); - - AliTPCDigitsArray *arr = new AliTPCDigitsArray; - arr->SetClass("AliSimDigits"); - arr->Setup(fTPCParam); - arr->MakeTree(fDigitsFile); - SetDigitsArray(arr); - - cerr<<"Digitizing TPC...\n"; - - // fDigitsSwitch=1; // summable digits -for the moment direct - - for(Int_t isec=0;isecGetNSector();isec++) Hits2DigitsSector(isec); - - - // write results - char treeName[100]; - - sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); - - GetDigitsArray()->GetTree()->SetName(treeName); - GetDigitsArray()->GetTree()->Write(); + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + SetTreeAddress(); + SetActiveSectors(); + Hits2SDigits2(iEvent); + } + fLoader->UnloadHits(); + fLoader->UnloadSDigits(); } - - //_____________________________________________________________________________ + void AliTPC::Hits2DigitsSector(Int_t isec) { //------------------------------------------------------------------- @@ -1616,7 +1326,13 @@ void AliTPC::Hits2DigitsSector(Int_t isec) if(fDefaults == 0) SetDefaults(); - TTree *tH = gAlice->TreeH(); // pointer to the hits tree + TTree *tH = TreeH(); // pointer to the hits tree + if (tH == 0x0) + { + Fatal("Hits2DigitsSector","Can not find TreeH in folder"); + return; + } + Stat_t ntracks = tH->GetEntries(); if( ntracks > 0){ @@ -1631,20 +1347,23 @@ void AliTPC::Hits2DigitsSector(Int_t isec) Int_t nrows =fTPCParam->GetNRow(isec); - row= new TObjArray* [nrows]; + row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk MakeSector(isec,nrows,tH,ntracks,row); //-------------------------------------------------------- // Digitize this sector, row by row - // row[i] is the pointer to the TObjArray of TVectors, + // row[i] is the pointer to the TObjArray of AliTPCFastVectors, // each one containing electrons accepted on this // row, assigned into tracks //-------------------------------------------------------- Int_t i; - if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree(fDigitsFile); + if (fDigitsArray->GetTree()==0) + { + Fatal("Hits2DigitsSector","Tree not set in fDigitsArray"); + } for (i=0;iStoreRow(isec,i); - Int_t ndig = dig->GetDigitSize(); - - - //printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig); + Int_t ndig = dig->GetDigitSize(); + if (gDebug > 10) + printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig); + fDigitsArray->ClearRow(isec,i); } // end of the sector digitization - for(i=0;iDelete(); delete row[i]; } @@ -1691,7 +1410,7 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) Float_t zerosup = fTPCParam->GetZeroSup(); - Int_t nrows =fTPCParam->GetNRow(isec); + // Int_t nrows =fTPCParam->GetNRow(isec); fCurrentIndex[1]= isec; @@ -1702,10 +1421,11 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) // Integrated signal for this row // and a single track signal // - TMatrix *m1 = new TMatrix(0,nofPads,0,nofTbins); // integrated - TMatrix *m2 = new TMatrix(0,nofPads,0,nofTbins); // single + + AliTPCFastMatrix *m1 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // integrated + AliTPCFastMatrix *m2 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // single // - TMatrix &total = *m1; + AliTPCFastMatrix &total = *m1; // Array of pointers to the label-signal list @@ -1718,14 +1438,16 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) // //calculate signal // - Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0); - Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1); + //Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0); + //Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1); + Int_t row1=irow; + Int_t row2=irow+2; for (Int_t row= row1;row<=row2;row++){ Int_t nTracks= rows[row]->GetEntries(); for (i1=0;i1Zero(); // clear single track signal matrix Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange); GetList(trackLabel,nofPads,m2,indexRange,pList); @@ -1737,28 +1459,27 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) Int_t tracks[3]; AliDigits *dig = fDigitsArray->GetRow(isec,irow); - for(Int_t ip=0;ipGaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); - - q = (Int_t)q; - - if(q <=zerosup) continue; // do not fill zeros + q+=GetNoise(); + if(q <=fzerosup) continue; // do not fill zeros + q = TMath::Nint(q); if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation } else { if(q <= 0.) continue; // do not fill zeros + if(q>2000.) q=2000.; q *= 16.; - q = (Int_t)q; + q = TMath::Nint(q); } // @@ -1766,7 +1487,7 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) // 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; } //Begin_Html @@ -1805,8 +1526,8 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) //_____________________________________________________________________________ -Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2, - Int_t *indexRange) +Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, + AliTPCFastMatrix *m1, AliTPCFastMatrix *m2,Int_t *indexRange) { //--------------------------------------------------------------- @@ -1820,13 +1541,13 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2, // Modified: Marian Ivanov //----------------------------------------------------------------- - TVector *tv; - - tv = (TVector*)p1->At(ntr); // pointer to a track - TVector &v = *tv; + AliTPCFastVector *tv; + + tv = (AliTPCFastVector*)p1->At(ntr); // pointer to a track + AliTPCFastVector &v = *tv; Float_t label = v(0); - Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2; + Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2; Int_t nElectrons = (tv->GetNrows()-1)/4; indexRange[0]=9999; // min pad @@ -1834,10 +1555,8 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2, indexRange[2]=9999; //min time indexRange[3]=-1; // max time - // Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above - - TMatrix &signal = *m1; - TMatrix &total = *m2; + AliTPCFastMatrix &signal = *m1; + AliTPCFastMatrix &total = *m2; // // Loop over all electrons // @@ -1846,22 +1565,28 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2, Float_t aval = v(idx+4); Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)}; - Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]); - - if (n>0) for (Int_t i =0; iGetResBin(i); + Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]); + + Int_t *index = fTPCParam->GetResBin(0); + Float_t *weight = & (fTPCParam->GetResWeight(0)); + + if (n>0) for (Int_t i =0; iGetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) { + + if (pad>=0){ Int_t time=index[2]; - Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel - weight *= eltoadcfac; + Float_t qweight = *(weight)*eltoadcfac; - if (m1!=0) signal(pad,time)+=weight; - total(pad,time)+=weight; - indexRange[0]=TMath::Min(indexRange[0],pad); - indexRange[1]=TMath::Max(indexRange[1],pad); - indexRange[2]=TMath::Min(indexRange[2],time); - indexRange[3]=TMath::Max(indexRange[3],time); + if (m1!=0) signal.UncheckedAt(pad,time)+=qweight; + total.UncheckedAt(pad,time)+=qweight; + if (indexRange[0]>pad) indexRange[0]=pad; + if (indexRange[1]time) indexRange[2]=time; + if (indexRange[3]GetBranch("TPC2"); + 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. + // each TObjArray will store the AliTPCFastVectors + // of electrons, one AliTPCFastVectors per each track. //---------------------------------------------- - Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row - TVector **tracks = new TVector* [nrows]; //pointers to the track vectors - for(i=0; iGetBranch("fTrackHitsInfo"); - br->GetEvent(track); - AliObjectArray * ar = fTrackHits->fTrackHitsInfo; - for (UInt_t j=0;jGetSize();j++){ - if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE; - } - } + isInSector = TrackInVolume(isec,track); if (!isInSector) continue; //MI change branch->GetEntry(track); // get next track @@ -2039,7 +1756,6 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, while(tpcHit){ Int_t sector=tpcHit->fSector; // sector number - // if(sector != isec) continue; if(sector != isec){ tpcHit = (AliTPChit*) NextHit(); continue; @@ -2052,22 +1768,22 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, // store already filled fTrack - for(i=0;i0){ - TVector &v = *tracks[i]; + AliTPCFastVector &v = *tracks[i]; v(0) = previousTrack; tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary row[i]->Add(tracks[i]); } else{ - delete tracks[i]; // delete empty TVector + delete tracks[i]; // delete empty AliTPCFastVector tracks[i]=0; } } nofElectrons[i]=0; - tracks[i] = new TVector(481); // TVectors for the next fTrack + tracks[i] = new AliTPCFastVector(481); // AliTPCFastVectors for the next fTrack } // end of loop over rows @@ -2105,13 +1821,29 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); index[0]=1; - TransportElectron(xyz,index); //MI change -august + TransportElectron(xyz,index); Int_t rowNumber; - fTPCParam->GetPadRow(xyz,index); //MI change august - rowNumber = index[2]; + fTPCParam->GetPadRow(xyz,index); + // 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)) continue; + 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 @@ -2124,13 +1856,11 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, } } - TVector &v = *tracks[rowNumber]; + AliTPCFastVector &v = *tracks[rowNumber]; Int_t idx = 4*nofElectrons[rowNumber]-3; - - v(idx)= xyz[0]; // X - pad row coordinate - v(idx+1)=xyz[1]; // Y - pad coordinate (along the pad-row) - v(idx+2)=xyz[2]; // Z - time bin coordinate - v(idx+3)=xyz[3]; // avalanche size + Real_t * position = &(((AliTPCFastVector&)v).UncheckedAt(idx)); //make code faster + memcpy(position,xyz,4*sizeof(Float_t)); + } // end of loop over electrons tpcHit = (AliTPChit*)NextHit(); @@ -2142,9 +1872,9 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, // store remaining track (the last one) if not empty // - for(i=0;i0){ - TVector &v = *tracks[i]; + AliTPCFastVector &v = *tracks[i]; v(0) = previousTrack; tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary row[i]->Add(tracks[i]); @@ -2183,25 +1913,30 @@ void AliTPC::Init() } //_____________________________________________________________________________ -void AliTPC::MakeBranch(Option_t* option, const char *file) +void AliTPC::MakeBranch(Option_t* option) { // // Create Tree branches for the TPC. // + if(GetDebug()) Info("MakeBranch",""); Int_t buffersize = 4000; char branchname[10]; sprintf(branchname,"%s",GetName()); + + const char *h = strstr(option,"H"); - AliDetector::MakeBranch(option,file); + if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03 + + AliDetector::MakeBranch(option); const char *d = strstr(option,"D"); + + if (fDigits && fLoader->TreeD() && d) + { + MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0); + } - if (fDigits && gAlice->TreeD() && d) { - MakeBranchInTree(gAlice->TreeD(), - branchname, &fDigits, buffersize, file); - } - - if (fHitType&2) MakeBranch2(option,file); // MI change 14.09.2000 + if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000 } //_____________________________________________________________________________ @@ -2224,7 +1959,6 @@ void AliTPC::SetSecAL(Int_t sec) //----------------------------------------------------------------- // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl //----------------------------------------------------------------- - fSecAL = sec; } @@ -2238,7 +1972,6 @@ void AliTPC::SetSecAU(Int_t sec) //----------------------------------------------------------------- // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl //----------------------------------------------------------------- - fSecAU = sec; } @@ -2359,32 +2092,12 @@ void AliTPC::TransportElectron(Float_t *xyz, Int_t *index) // ExB if (fTPCParam->GetMWPCReadout()==kTRUE){ - Float_t x1=xyz[0]; - fTPCParam->Transform2to2NearestWire(xyz,index); - Float_t dx=xyz[0]-x1; + Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index); xyz[1]+=dx*(fTPCParam->GetOmegaTau()); } - //add nonisochronity (not implemented yet) - + //add nonisochronity (not implemented yet) } -ClassImp(AliTPCdigit) - -//_____________________________________________________________________________ -AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits): - AliDigit(tracks) -{ - // - // Creates a TPC digit object - // - fSector = digits[0]; - fPadRow = digits[1]; - fPad = digits[2]; - fTime = digits[3]; - fSignal = digits[4]; -} - - ClassImp(AliTPChit) //_____________________________________________________________________________ @@ -2402,52 +2115,47 @@ AliHit(shunt,track) fQ = hits[3]; } - //________________________________________________________________________ -// Additional code because of the AliTPCTrackHits +// Additional code because of the AliTPCTrackHitsV2 -void AliTPC::MakeBranch2(Option_t *option,const char *file) +void AliTPC::MakeBranch2(Option_t *option,const char */*file*/) { // // Create a new branch in the current Root Tree // The branch of fHits is automatically split // MI change 14.09.2000 - if (fHitType&2==0) return; + if(GetDebug()) Info("MakeBranch2",""); + if (fHitType<2) return; char branchname[10]; sprintf(branchname,"%s2",GetName()); // // Get the pointer to the header const char *cH = strstr(option,"H"); // - if (fTrackHits && gAlice->TreeH() && cH) { - AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHits, - gAlice->TreeH(),fBufferSize,99); - gAlice->TreeH()->GetListOfBranches()->Add(branch); - if (GetDebug()>1) - printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname); - const char folder [] = "RunMC/Event/Data"; - if (GetDebug()) - printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder); - Publish(folder,&fTrackHits,branchname); - if (file) { - TBranch *b = gAlice->TreeH()->GetBranch(branchname); - TDirectory *wd = gDirectory; - b->SetFile(file); - TIter next( b->GetListOfBranches()); - while ((b=(TBranch*)next())) { - b->SetFile(file); - } - wd->cd(); - if (GetDebug()>1) - cout << "Diverting branch " << branchname << " to file " << file << endl; - } - } + if (fTrackHits && TreeH() && cH && fHitType&4) + { + if(GetDebug()) Info("MakeBranch2","Making branch for Type 4 Hits"); + TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99); + } + + if (fTrackHitsOld && TreeH() && cH && fHitType&2) + { + if(GetDebug()) Info("MakeBranch2","Making branch for Type 2 Hits"); + AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld, + TreeH(),fBufferSize,99); + TreeH()->GetListOfBranches()->Add(branch); + } } void AliTPC::SetTreeAddress() { - if (fHitType&1) AliDetector::SetTreeAddress(); - if (fHitType&2) SetTreeAddress2(); +//Sets tree address for hits + if (fHitType<=1) + { + if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03 + AliDetector::SetTreeAddress(); + } + if (fHitType>1) SetTreeAddress2(); } void AliTPC::SetTreeAddress2() @@ -2455,21 +2163,49 @@ void AliTPC::SetTreeAddress2() // // Set branch address for the TrackHits Tree // + if(GetDebug()) Info("SetTreeAddress2",""); + TBranch *branch; char branchname[20]; sprintf(branchname,"%s2",GetName()); // // Branch address for hit tree - TTree *treeH = gAlice->TreeH(); - if (treeH) { + TTree *treeH = TreeH(); + if ((treeH)&&(fHitType&4)) { branch = treeH->GetBranch(branchname); - if (branch) branch->SetAddress(&fTrackHits); + if (branch) + { + branch->SetAddress(&fTrackHits); + if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Setting"); + } + else + if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Failed (can not find branch)"); + } + if ((treeH)&&(fHitType&2)) { + branch = treeH->GetBranch(branchname); + if (branch) + { + branch->SetAddress(&fTrackHitsOld); + if (GetDebug()) Info("SetTreeAddress2","fHitType&2 Setting"); + } + else if (GetDebug()) + Info("SetTreeAddress2","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); + } + } void AliTPC::FinishPrimary() { - if (fTrackHits) fTrackHits->FlushHitStack(); + if (fTrackHits &&fHitType&4) fTrackHits->FlushHitStack(); + if (fTrackHitsOld && fHitType&2) fTrackHitsOld->FlushHitStack(); } @@ -2479,42 +2215,52 @@ void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits) // add hit to the list Int_t rtrack; if (fIshunt) { - int primary = gAlice->GetPrimary(track); - gAlice->Particle(primary)->SetBit(kKeepBit); + int primary = gAlice->GetMCApp()->GetPrimary(track); + gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit); rtrack=primary; } else { rtrack=track; - gAlice->FlagTrack(track); + gAlice->GetMCApp()->FlagTrack(track); } //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1); //if (hit->fTrack!=rtrack) // cout<<"bad track number\n"; - if (fTrackHits) + if (fTrackHits && fHitType&4) fTrackHits->AddHitKartez(vol[0],rtrack, hits[0], - hits[1],hits[2],(Int_t)hits[3]); + hits[1],hits[2],(Int_t)hits[3]); + if (fTrackHitsOld &&fHitType&2 ) + fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0], + hits[1],hits[2],(Int_t)hits[3]); + } void AliTPC::ResetHits() -{ +{ if (fHitType&1) AliDetector::ResetHits(); - if (fHitType&2) ResetHits2(); + if (fHitType>1) ResetHits2(); } void AliTPC::ResetHits2() { // //reset hits - if (fTrackHits) fTrackHits->Clear(); + if (fTrackHits && fHitType&4) fTrackHits->Clear(); + if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear(); + } AliHit* AliTPC::FirstHit(Int_t track) { - if (fHitType&2) return FirstHit2(track); + if (fHitType>1) return FirstHit2(track); return AliDetector::FirstHit(track); } AliHit* AliTPC::NextHit() { - if (fHitType&2) return NextHit2(); + // + // gets next hit + // + if (fHitType>1) return NextHit2(); + return AliDetector::NextHit(); } @@ -2529,13 +2275,18 @@ AliHit* AliTPC::FirstHit2(Int_t track) // if(track>=0) { gAlice->ResetHits(); - gAlice->TreeH()->GetEvent(track); + TreeH()->GetEvent(track); } // - if (fTrackHits) { + if (fTrackHits && fHitType&4) { fTrackHits->First(); return fTrackHits->GetHit(); } + if (fTrackHitsOld && fHitType&2) { + fTrackHitsOld->First(); + return fTrackHitsOld->GetHit(); + } + else return 0; } @@ -2544,6 +2295,11 @@ AliHit* AliTPC::NextHit2() // //Return the next hit for the current track + + if (fTrackHitsOld && fHitType&2) { + fTrackHitsOld->Next(); + return fTrackHitsOld->GetHit(); + } if (fTrackHits) { fTrackHits->Next(); return fTrackHits->GetHit(); @@ -2569,15 +2325,68 @@ void AliTPC::LoadPoints(Int_t) void AliTPC::RemapTrackHitIDs(Int_t *map) { + // + // remapping + // if (!fTrackHits) return; - AliObjectArray * arr = fTrackHits->fTrackHitsInfo; - for (UInt_t i=0;iGetSize();i++){ - AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i)); - info->fTrackID = map[info->fTrackID]; - } + 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){ + TClonesArray * arr = fTrackHits->GetArray();; + for (Int_t i=0;iGetEntriesFast();i++){ + AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i)); + info->fTrackID = map[info->fTrackID]; + } + } } +Bool_t AliTPC::TrackInVolume(Int_t id,Int_t track) +{ + //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) { + printf("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; + +} //_____________________________________________________________________________ void AliTPC::LoadPoints2(Int_t) @@ -2585,11 +2394,14 @@ void AliTPC::LoadPoints2(Int_t) // // Store x, y, z of all hits in memory // - if (fTrackHits == 0) return; + if (fTrackHits == 0 && fTrackHitsOld==0) return; // - Int_t nhits = fTrackHits->GetEntriesFast(); + Int_t nhits =0; + if (fHitType&4) nhits = fTrackHits->GetEntriesFast(); + if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast(); + if (nhits == 0) return; - Int_t tracks = gAlice->GetNtrack(); + Int_t tracks = gAlice->GetMCApp()->GetNtrack(); if (fPoints == 0) fPoints = new TObjArray(tracks); AliHit *ahit; // @@ -2610,9 +2422,7 @@ void AliTPC::LoadPoints2(Int_t) // Loop over all the hits and store their position // ahit = FirstHit2(-1); - //for (Int_t hit=0;hitUncheckedAt(hit); trk=ahit->GetTrack(); if(ntrk[trk]==limi[trk]) { // @@ -2633,6 +2443,9 @@ void AliTPC::LoadPoints2(Int_t) ntrk[trk]++; ahit = NextHit2(); } + + + // for(trk=0; trkGetEntriesFast(); if (nhits == 0) return; - Int_t tracks = gAlice->GetNtrack(); + Int_t tracks = gAlice->GetMCApp()->GetNtrack(); if (fPoints == 0) fPoints = new TObjArray(2*tracks); fPoints->Expand(2*tracks); AliHit *ahit; @@ -2741,159 +2554,57 @@ void AliTPC::LoadPoints3(Int_t) -void AliTPC::FindTrackHitsIntersection(TClonesArray * arr) +AliLoader* AliTPC::MakeLoader(const char* topfoldername) { - - // - //fill clones array with intersection of current point with the - //middle of the row - Int_t sector; - Int_t ipart; - - const Int_t kcmaxhits=30000; - TVector * xxxx = new TVector(kcmaxhits*4); - TVector & xxx = *xxxx; - Int_t maxhits = kcmaxhits; - - // - AliTPChit * tpcHit=0; - tpcHit = (AliTPChit*)FirstHit2(-1); - Int_t currentIndex=0; - Int_t lastrow=-1; //last writen row - - while (tpcHit){ - if (tpcHit==0) continue; - sector=tpcHit->fSector; // sector number - ipart=tpcHit->Track(); - - //find row number - - Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()}; - Int_t index[3]={1,sector,0}; - Int_t currentrow = fTPCParam->GetPadRow(x,index) ; - if (currentrow<0) continue; - if (lastrow<0) lastrow=currentrow; - if (currentrow==lastrow){ - if ( currentIndex>=maxhits){ - maxhits+=kcmaxhits; - xxx.ResizeTo(4*maxhits); - } - xxx(currentIndex*4)=x[0]; - xxx(currentIndex*4+1)=x[1]; - xxx(currentIndex*4+2)=x[2]; - xxx(currentIndex*4+3)=tpcHit->fQ; - currentIndex++; - } - else - if (currentIndex>2){ - Float_t sumx=0; - Float_t sumx2=0; - Float_t sumx3=0; - Float_t sumx4=0; - Float_t sumy=0; - Float_t sumxy=0; - Float_t sumx2y=0; - Float_t sumz=0; - Float_t sumxz=0; - Float_t sumx2z=0; - Float_t sumq=0; - for (Int_t index=0;indexGetNPads(sector,lastrow)-1)/2; - Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+ - sumx2*(sumx*sumx3-sumx2*sumx2); - - Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+ - sumx2*(sumxy*sumx3-sumx2y*sumx2); - Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+ - sumx2*(sumxz*sumx3-sumx2z*sumx2); - - Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+ - sumx2*(sumx*sumx2y-sumx2*sumxy); - Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+ - sumx2*(sumx*sumx2z-sumx2*sumxz); - - Float_t y=detay/det+centralPad; - Float_t z=detaz/det; - Float_t by=detby/det; //y angle - Float_t bz=detbz/det; //z angle - sumy/=Float_t(currentIndex); - sumz/=Float_t(currentIndex); - - AliComplexCluster cl; - cl.fX=z; - cl.fY=y; - cl.fQ=sumq; - cl.fSigmaX2=bz; - cl.fSigmaY2=by; - cl.fTracks[0]=ipart; - - AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow)); - if (row!=0) row->InsertCluster(&cl); - currentIndex=0; - lastrow=currentrow; - } //end of calculating cluster for given row - - } // end of loop over hits - xxxx->Delete(); - +//Makes TPC loader + fLoader = new AliTPCLoader(GetName(),topfoldername); + return fLoader; } -//_______________________________________________________________________________ -void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent) -{ - // produces rec points from digits and writes them on the root file - // AliTPCclusters.root - - TDirectory *cwd = gDirectory; - - AliTPCParam *dig=(AliTPCParam *)gDirectory->Get("75x40_100x60"); - SetParam(dig); - cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<Getenv("CONFIG_FILE")){ - out=TFile::Open("AliTPCclusters.root","recreate"); - } - else { - const char *tmp1; - const char *tmp2; - char tmp3[80]; - tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX"); - tmp2=gSystem->Getenv("CONFIG_OUTDIR"); - sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2); - out=TFile::Open(tmp3,"recreate"); +//////////////////////////////////////////////////////////////////////// +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) { + cout<<"TPC parameters "<cd(); - for(Int_t iev=firstevent;ievGet(paramName); +// if (paramTPC) { +// cout<<"TPC parameters "<Get(paramName); +// if (paramTPC) { +// cout<<"TPC parameters "<Close(); - cwd->cd(); +} -}