X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPC.cxx;h=1820b73555fe7f6b99796f4781d6da56491d8adc;hb=8f4188a7c02eb02e669af8d546ab1efbe895fd6f;hp=9c5d75c535ceff01bfa261638ac69e531f5fd177;hpb=f03e342345a8fa1862aa7bee50b2e9af9aabda45;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPC.cxx b/TPC/AliTPC.cxx index 9c5d75c535c..1820b73555f 100644 --- a/TPC/AliTPC.cxx +++ b/TPC/AliTPC.cxx @@ -13,151 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.53 2002/02/25 11:02:56 kowal2 -Changes towards speeding up the code. Thanks to Marian Ivanov. - -Revision 1.52 2002/02/18 09:26:09 kowal2 -Removed compiler warning - -Revision 1.51 2002/01/21 17:13:21 kowal2 -New track hits using root containers. Setting active sectors added. - -Revision 1.50 2001/12/06 14:16:19 kowal2 -meaningfull printouts - -Revision 1.49 2001/11/30 11:55:37 hristov -Noise table created in Hits2SDigits (M.Ivanov) - -Revision 1.48 2001/11/24 16:10:21 kowal2 -Faster algorithms. - -Revision 1.47 2001/11/19 10:25:34 kowal2 -Nearest integer instead of integer when converting to ADC counts - -Revision 1.46 2001/11/07 06:47:12 kowal2 -Removed printouts - -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$ */ /////////////////////////////////////////////////////////////////////////////// // // @@ -177,105 +33,108 @@ Introduction of the Copyright and cvs Log // +#include +#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 "AliRun.h" -#include -#include -#include -#include "AliMC.h" -#include "AliMagF.h" - +#include +#include +#include -#include "AliTPCParamSR.h" -#include "AliTPCPRF2D.h" -#include "AliTPCRF1D.h" #include "AliDigits.h" +#include "AliMagF.h" +#include "AliRun.h" +#include "AliRunLoader.h" #include "AliSimDigits.h" -#include "AliTPCTrackHits.h" -#include "AliTPCTrackHitsV2.h" -#include "AliPoints.h" -#include "AliArrayBranch.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 "AliTPCTrackHitsV2.h" +#include "AliTrackReference.h" +#include "AliMC.h" +#include "AliStack.h" +#include "AliTPCDigitizer.h" +#include "AliTPCBuffer.h" +#include "AliTPCDDLRawData.h" +#include "AliLog.h" +#include "AliTPCcalibDB.h" +#include "AliTPCCalPad.h" +#include "AliTPCCalROC.h" +#include "AliTPCExB.h" +#include "AliRawReader.h" +#include "AliTPCRawStream.h" +#include "TTreeStream.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 row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb); - inline Float_t & UncheckedAt(Int_t rown, Int_t coln) const {return (fIndex[coln])[rown];} //fast acces - inline Float_t UncheckedAtFast(Int_t rown, Int_t coln) const {return (fIndex[coln])[rown];} //fast acces -}; - -AliTPCFastMatrix::AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb): - TMatrix(row_lwb, row_upb,col_lwb,col_upb) - { - }; //_____________________________________________________________________________ -class AliTPCFastVector : public TVector { -public : - AliTPCFastVector(Int_t size); - inline Float_t & UncheckedAt(Int_t index) const {return fElements[index];} //fast acces -}; - -AliTPCFastVector::AliTPCFastVector(Int_t size):TVector(size){ -}; + AliTPC::AliTPC():AliDetector(), + fDefaults(0), + fSens(0), + fNsectors(0), + fDigitsArray(0), + fTPCParam(0), + fTrackHits(0), + fHitType(0), + fDigitsSwitch(0), + fSide(0), + fPrimaryIonisation(0), + fNoiseDepth(0), + fNoiseTable(0), + fCurrentNoise(0), + fActiveSectors(0), + fGainFactor(1.), + fDebugStreamer(0) -//_____________________________________________________________________________ -AliTPC::AliTPC() { // // Default constructor // fIshunt = 0; - fHits = 0; - fDigits = 0; - fNsectors = 0; - fDigitsArray = 0; - fClustersArray = 0; - fDefaults = 0; - fTrackHits = 0; - fTrackHitsOld = 0; - fHitType = 4; //default CONTAINERS - based on ROOT structure - fTPCParam = 0; - fNoiseTable = 0; - fActiveSectors =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 } //_____________________________________________________________________________ AliTPC::AliTPC(const char *name, const char *title) - : AliDetector(name,title) + : AliDetector(name,title), + fDefaults(0), + fSens(0), + fNsectors(0), + fDigitsArray(0), + fTPCParam(0), + fTrackHits(0), + fHitType(0), + fDigitsSwitch(0), + fSide(0), + fPrimaryIonisation(0), + fNoiseDepth(0), + fNoiseTable(0), + fCurrentNoise(0), + fActiveSectors(0), + fGainFactor(1.), + fDebugStreamer(0) + { // // Standard constructor @@ -284,34 +143,32 @@ AliTPC::AliTPC(const char *name, const char *title) // // Initialise arrays of hits and digits fHits = new TClonesArray("AliTPChit", 176); - gAlice->AddHitList(fHits); - fDigitsArray = 0; - fClustersArray= 0; - fDefaults = 0; + gAlice->GetMCApp()->AddHitList(fHits); // 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); + //fTrackHitsOld = new AliTPCTrackHits; //MI - 13.09.2000 + //fTrackHitsOld->SetHitPrecision(0.002); + //fTrackHitsOld->SetStepPrecision(0.003); + //fTrackHitsOld->SetMaxDistance(100); + + +#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1) + fHitType = 4; // ROOT containers +#else + fHitType = 2; +#endif - fNoiseTable =0; - fHitType = 4; - fActiveSectors = 0; - // - // Initialise counters - fNsectors = 0; // fIshunt = 0; // // Initialise color attributes - SetMarkerColor(kYellow); + //PH SetMarkerColor(kYellow); // // Set TPC parameters @@ -319,9 +176,10 @@ AliTPC::AliTPC(const char *name, const char *title) if (!strcmp(title,"Default")) { - fTPCParam = new AliTPCParamSR; + //fTPCParam = new AliTPCParamSR; + fTPCParam = AliTPCcalibDB::Instance()->GetParameters(); } else { - cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n"; + AliWarning("In Config.C you must set non-default parameters."); fTPCParam=0; } @@ -337,11 +195,14 @@ AliTPC::~AliTPC() fIshunt = 0; delete fHits; delete fDigits; - delete fTPCParam; + //delete fTPCParam; delete fTrackHits; //MI 15.09.2000 - delete fTrackHitsOld; //MI 10.12.2001 - if (fNoiseTable) delete [] fNoiseTable; - + // delete fTrackHitsOld; //MI 10.12.2001 + + fDigitsArray = 0x0; + delete [] fNoiseTable; + delete [] fActiveSectors; + if (fDebugStreamer) delete fDebugStreamer; } //_____________________________________________________________________________ @@ -350,117 +211,12 @@ 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); 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 *nNode, *nTop; - TTUBS *tubs; - Int_t i; - const int kColorTPC=19; - char name[5], title[25]; - const Double_t kDegrad=TMath::Pi()/180; - 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 - // - - nTop=gAlice->GetGeometry()->GetNode("alice"); - - // inner sectors - - rl = fTPCParam->GetInnerRadiusLow(); - ru = fTPCParam->GetInnerRadiusUp(); - - - for(i=0;iSetNumberOfDivisions(1); - 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); - 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 ) -{ - // - // Calculate distance from TPC to mouse on the display - // Dummy procedure - // - return 9999; -} - -void AliTPC::Clusters2Tracks(TFile *of) { - //----------------------------------------------------------------- - // This is a track finder. - //----------------------------------------------------------------- - AliTPCtracker tracker(fTPCParam); - tracker.Clusters2Tracks(gFile,of); + AddHit2(track,vol,hits); } //_____________________________________________________________________________ @@ -474,51 +230,24 @@ void AliTPC::CreateMaterials() // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl //----------------------------------------------------------------- - Int_t iSXFLD=gAlice->Field()->Integ(); - Float_t sXMGMX=gAlice->Field()->Max(); + Int_t iSXFLD=((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ(); + Float_t sXMGMX=((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max(); Float_t amat[5]; // atomic numbers Float_t zmat[5]; // z Float_t wmat[5]; // proportions Float_t density; - Float_t apure[2]; - - - //***************** Gases ************************* - - //------------------------------------------------- - // pure gases - //------------------------------------------------- - - // Neon - - - amat[0]= 20.18; - zmat[0]= 10.; - density = 0.0009; - apure[0]=amat[0]; - AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.); - // Argon - - amat[0]= 39.948; - zmat[0]= 18.; - density = 0.001782; - - apure[1]=amat[0]; + //***************** Gases ************************* - AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.); - //-------------------------------------------------------------- - // gases - compounds + // gases - air and CO2 //-------------------------------------------------------------- - Float_t amol[3]; - // CO2 amat[0]=12.011; @@ -527,106 +256,57 @@ void AliTPC::CreateMaterials() zmat[0]=6.; zmat[1]=8.; - wmat[0]=1.; - wmat[1]=2.; + wmat[0]=0.2729; + wmat[1]=0.7271; density=0.001977; - amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1]; - - AliMixture(10,"CO2",amat,zmat,density,-2,wmat); - - // CF4 - - amat[0]=12.011; - amat[1]=18.998; - - zmat[0]=6.; - zmat[1]=9.; - - wmat[0]=1.; - wmat[1]=4.; - - density=0.003034; - - amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1]; - - AliMixture(11,"CF4",amat,zmat,density,-2,wmat); - - - // CH4 - - amat[0]=12.011; - amat[1]=1.; - - zmat[0]=6.; - zmat[1]=1.; - - wmat[0]=1.; - wmat[1]=4.; - density=0.000717; - - amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1]; - - AliMixture(12,"CH4",amat,zmat,density,-2,wmat); + AliMixture(10,"CO2",amat,zmat,density,2,wmat); + // + // Air + // + amat[0]=15.9994; + amat[1]=14.007; + // + zmat[0]=8.; + zmat[1]=7.; + // + wmat[0]=0.233; + wmat[1]=0.767; + // + density=0.001205; + AliMixture(11,"Air",amat,zmat,density,2,wmat); + //---------------------------------------------------------------- - // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds + // drift gases //---------------------------------------------------------------- - char namate[21]; - density = 0.; - Float_t am=0; - Int_t nc; - Float_t rho,absl,X0,buf[1]; - Int_t nbuf; - Float_t a,z; - - for(nc = 0;ncGfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf); - - amat[nc] = a; - zmat[nc] = z; - - Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10; - - am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]); - density += fMixtProp[nc]*rho; // density of the mixture - - } - - // mixture proportions by weight! - - for(nc = 0;nc=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10; - - wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? - apure[nnc] : amol[nnc])/am; - - } // Drift gases 1 - nonsensitive, 2 - sensitive + // Ne-CO2-N (85-10-5) - AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat); - AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat); - - - // Air - - amat[0] = 14.61; - zmat[0] = 7.3; - density = 0.001205; + amat[0]= 20.18; + amat[1]=12.011; + amat[2]=15.9994; + amat[3]=14.007; - AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.); + 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; + AliMixture(12,"Ne-CO2-N-1",amat,zmat,density,4,wmat); + AliMixture(13,"Ne-CO2-N-2",amat,zmat,density,4,wmat); + AliMixture(30,"Ne-CO2-N-3",amat,zmat,density,4,wmat); //---------------------------------------------------------------------- // solid materials //---------------------------------------------------------------------- @@ -651,7 +331,7 @@ void AliTPC::CreateMaterials() density = 1.45; - AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat); + AliMixture(14,"Kevlar",amat,zmat,density,-4,wmat); // NOMEX @@ -670,10 +350,9 @@ void AliTPC::CreateMaterials() wmat[2] = 2.; wmat[3] = 2.; - density = 0.03; - - - AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat); + density = 0.029; + + AliMixture(15,"NOMEX",amat,zmat,density,-4,wmat); // Makrolon C16H18O3 @@ -691,7 +370,25 @@ void AliTPC::CreateMaterials() density = 1.2; - AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat); + AliMixture(16,"Makrolon",amat,zmat,density,-3,wmat); + + // Tedlar C2H3F + + amat[0] = 12.011; + amat[1] = 1.; + amat[2] = 18.998; + + zmat[0] = 6.; + zmat[1] = 1.; + zmat[2] = 9.; + + wmat[0] = 2.; + wmat[1] = 3.; + wmat[2] = 1.; + + density = 1.71; + + AliMixture(17, "Tedlar",amat,zmat,density,-3,wmat); // Mylar C5H4O2 @@ -709,22 +406,91 @@ void AliTPC::CreateMaterials() density = 1.39; - AliMixture(37, "Mylar",amat,zmat,density,-3,wmat); + AliMixture(18, "Mylar",amat,zmat,density,-3,wmat); + // material for "prepregs" + // Epoxy - C14 H20 O3 + // Quartz SiO2 + // Carbon C + // prepreg1 60% C-fiber, 40% epoxy (vol) + amat[0]=12.011; + amat[1]=1.; + amat[2]=15.994; - // SiO2 - used later for the glass fiber + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; - amat[0]=28.086; - amat[1]=15.9994; + wmat[0]=0.923; + wmat[1]=0.023; + wmat[2]=0.054; - zmat[0]=14.; - zmat[1]=8.; + density=1.859; + + AliMixture(19, "Prepreg1",amat,zmat,density,3,wmat); - wmat[0]=1.; - wmat[1]=2.; + //prepreg2 60% glass-fiber, 40% epoxy + + amat[0]=12.01; + amat[1]=1.; + amat[2]=15.994; + amat[3]=28.086; + + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; + zmat[3]=14.; + + wmat[0]=0.194; + wmat[1]=0.023; + wmat[2]=0.443; + wmat[3]=0.340; + + density=1.82; + + AliMixture(20, "Prepreg2",amat,zmat,density,4,wmat); + + //prepreg3 50% glass-fiber, 50% epoxy + + amat[0]=12.01; + amat[1]=1.; + amat[2]=15.994; + amat[3]=28.086; + + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; + zmat[3]=14.; + + wmat[0]=0.225; + wmat[1]=0.03; + wmat[2]=0.443; + wmat[3]=0.3; + + density=1.163; + + AliMixture(21, "Prepreg3",amat,zmat,density,4,wmat); + + // G10 60% SiO2 40% epoxy + + amat[0]=12.01; + amat[1]=1.; + amat[2]=15.994; + amat[3]=28.086; + + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; + zmat[3]=14.; + wmat[0]=0.194; + wmat[1]=0.023; + wmat[2]=0.443; + wmat[3]=0.340; - AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2) + density=1.7; + AliMixture(22, "G10",amat,zmat,density,4,wmat); + // Al amat[0] = 26.98; @@ -732,16 +498,16 @@ void AliTPC::CreateMaterials() density = 2.7; - AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.); + AliMaterial(23,"Al",amat[0],zmat[0],density,999.,999.); - // Si + // Si (for electronics amat[0] = 28.086; zmat[0] = 14.; density = 2.33; - AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.); + AliMaterial(24,"Si",amat[0],zmat[0],density,999.,999.); // Cu @@ -750,29 +516,10 @@ void AliTPC::CreateMaterials() density = 8.96; - AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.); - - // Tedlar C2H3F - - amat[0] = 12.011; - amat[1] = 1.; - amat[2] = 18.998; - - zmat[0] = 6.; - zmat[1] = 1.; - zmat[2] = 9.; - - wmat[0] = 2.; - wmat[1] = 3.; - wmat[2] = 1.; - - density = 1.71; - - AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat); - - - // Plexiglas C5H8O2 + AliMaterial(25,"Cu",amat[0],zmat[0],density,999.,999.); + // Epoxy - C14 H20 O3 + amat[0]=12.011; amat[1]=1.; amat[2]=15.9994; @@ -781,17 +528,16 @@ void AliTPC::CreateMaterials() zmat[1]=1.; zmat[2]=8.; - wmat[0]=5.; - wmat[1]=8.; - wmat[2]=2.; + wmat[0]=14.; + wmat[1]=20.; + wmat[2]=3.; - density=1.18; + density=1.25; - AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat); + AliMixture(26,"Epoxy",amat,zmat,density,-3,wmat); - // Epoxy - C14 H20 O3 + // Plexiglas C5H8O2 - amat[0]=12.011; amat[1]=1.; amat[2]=15.9994; @@ -800,13 +546,13 @@ void AliTPC::CreateMaterials() zmat[1]=1.; zmat[2]=8.; - wmat[0]=14.; - wmat[1]=20.; - wmat[2]=3.; + wmat[0]=5.; + wmat[1]=8.; + wmat[2]=2.; - density=1.25; + density=1.18; - AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat); + AliMixture(27,"Plexiglas",amat,zmat,density,-3,wmat); // Carbon @@ -814,61 +560,102 @@ void AliTPC::CreateMaterials() zmat[0]=6.; density= 2.265; - AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.); + AliMaterial(28,"C",amat[0],zmat[0],density,999.,999.); + + // Fe (steel for the inner heat screen) + + amat[0]=55.845; - // get epoxy + zmat[0]=26.; - gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf); + density=7.87; - // Carbon fiber + AliMaterial(29,"Fe",amat[0],zmat[0],density,999.,999.); + // + // Peek - (C6H4-O-OC6H4-O-C6H4-CO)n + amat[0]=12.011; + amat[1]=1.; + amat[2]=15.9994; - wmat[0]=0.644; // by weight! - wmat[1]=0.356; + zmat[0]=6.; + zmat[1]=1.; + zmat[2]=8.; - density=0.5*(1.25+2.265); + wmat[0]=19.; + wmat[1]=12.; + wmat[2]=3.; + // + density=1.3; + // + AliMixture(30,"Peek",amat,zmat,density,-3,wmat); + // + // Ceramics - Al2O3 + // + amat[0] = 26.98; + amat[1]= 15.9994; - AliMixture(47,"Cfiber",amat,zmat,density,2,wmat); + zmat[0] = 13.; + zmat[1]=8.; + + wmat[0]=2.; + wmat[1]=3.; + + density = 3.97; - // get SiO2 + AliMixture(31,"Alumina",amat,zmat,density,-2,wmat); - gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf); + // + // liquids + // - wmat[0]=0.725; // by weight! - wmat[1]=0.275; + // water - density=1.7; + amat[0]=1.; + amat[1]=15.9994; - AliMixture(39,"G10",amat,zmat,density,2,wmat); + zmat[0]=1.; + zmat[1]=8.; - + wmat[0]=2.; + wmat[1]=1.; + density=1.; + AliMixture(32,"Water",amat,zmat,density,-2,wmat); + //---------------------------------------------------------- // tracking media for gases //---------------------------------------------------------- - AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1); - AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); - AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); + AliMedium(0, "Air", 11, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1); + AliMedium(1, "Ne-CO2-N-1", 12, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); + AliMedium(2, "Ne-CO2-N-2", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001); - + AliMedium(20, "Ne-CO2-N-3", 30, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001); //----------------------------------------------------------- // tracking media for solids //----------------------------------------------------------- - AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); - AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); - AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); - AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); - AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001); - AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); - + 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); + AliMedium(19,"Peek",30,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(21,"Alumina",31,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); + AliMedium(22,"Water",32,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); } void AliTPC::GenerNoise(Int_t tablesize) @@ -901,7 +688,7 @@ Float_t AliTPC::GetNoise() } -Bool_t AliTPC::IsSectorActive(Int_t sec) +Bool_t AliTPC::IsSectorActive(Int_t sec) const { // // check if the sector is active @@ -912,446 +699,341 @@ Bool_t AliTPC::IsSectorActive(Int_t sec) void AliTPC::SetActiveSectors(Int_t * sectors, Int_t n) { // activate interesting sectors - if (fActiveSectors) delete [] fActiveSectors; - fActiveSectors = new Bool_t[fTPCParam->GetNSector()]; + SetTreeAddress();//just for security + if (!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; } -void AliTPC::SetActiveSectors() +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 (!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 (fHitType>1) branch = gAlice->TreeH()->GetBranch("TPC2"); - else branch = gAlice->TreeH()->GetBranch("TPC"); - Stat_t ntracks = gAlice->TreeH()->GetEntries(); + if (fLoader->TreeH() == 0x0) + { + AliFatal("Can not find TreeH in folder"); + return; + } + if (fHitType>1) branch = fLoader->TreeH()->GetBranch("TPC2"); + else branch = fLoader->TreeH()->GetBranch("TPC"); + Stat_t ntracks = fLoader->TreeH()->GetEntries(); // loop over all hits - for(Int_t track=0;trackTreeH()->GetBranch("fVolumes"); - TBranch * br2 = gAlice->TreeH()->GetBranch("fNVolumes"); + TBranch * br1 = fLoader->TreeH()->GetBranch("fVolumes"); + TBranch * br2 = fLoader->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; + for (Int_t j=0;jGetNVolumes(); j++) { + if (volumes[j]>-1 && volumes[j]GetNSector()) { + fActiveSectors[volumes[j]]=kTRUE; + } + else { + AliError(Form("Volume %d -> sector number %d is outside (0..%d)", + j, + volumes[j], + fTPCParam->GetNSector())); + } + } } // - if (fTrackHitsOld && fHitType&2) { - TBranch * br = gAlice->TreeH()->GetBranch("fTrackHitsInfo"); - br->GetEvent(track); - AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo; - for (UInt_t j=0;jGetSize();j++){ - fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE; - } - } +// if (fTrackHitsOld && fHitType&2) { +// TBranch * br = fLoader->TreeH()->GetBranch("fTrackHitsInfo"); +// br->GetEvent(track); +// AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo; +// for (UInt_t j=0;jGetSize();j++){ +// fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE; +// } +// } } - } -void AliTPC::Digits2Clusters(TFile *of, Int_t eventnumber) -{ - //----------------------------------------------------------------- - // 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) +void AliTPC::Digits2Raw() { - //-------------------------------------------------------- - // TPC simple cluster generator from hits - // obtained from the TPC Fast Simulator - // The point errors are taken from the parametrization - //-------------------------------------------------------- +// convert digits of the current event to raw data - //----------------------------------------------------------------- - // 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; + static const Int_t kThreshold = 0; - TDirectory *savedir=gDirectory; - - if (!of->IsOpen()) { - cerr<<"AliTPC::Hits2Clusters(): output file not open !\n"; - return; + fLoader->LoadDigits(); + TTree* digits = fLoader->TreeD(); + if (!digits) { + AliError("No digits tree"); + return; } - //if(fDefaults == 0) SetDefaults(); - - Float_t sigmaRphi,sigmaZ,clRphi,clZ; - // - 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 - //--------------------------------------------------------------- - - TTree *tH = gAlice->TreeH(); - Stat_t ntracks = tH->GetEntries(); - - //Switch to the output file - of->cd(); + 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; + } - char cname[100]; + 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 - sprintf(cname,"TreeC_TPC_%d",eventn); + delete buffer; + fLoader->UnloadDigits(); - fTPCParam->Write(fTPCParam->GetTitle()); - AliTPCClustersArray carray; - carray.Setup(fTPCParam); - carray.SetClusterType("AliTPCcluster"); - carray.MakeTree(); + AliTPCDDLRawData rawWriter; + rawWriter.SetVerbose(0); - 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 - //------------------------------------------------------------ - - for(Int_t track=0;trackGetEvent(track); - // - // Get number of the TPC hits - // - tpcHit = (AliTPChit*)FirstHit(-1); + rawWriter.RawData(fileName); + gSystem->Unlink(fileName); - // Loop over hits - // - while(tpcHit){ - - if (tpcHit->fQ == 0.) { - tpcHit = (AliTPChit*) NextHit(); - continue; //information about track (I.Belikov) - } - sector=tpcHit->fSector; // sector number - - 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; irowTreeS() == 0x0 ) { + fLoader->MakeTree("S"); + } - } // end of loop over sectors + if(fDefaults == 0) SetDefaults(); // check if the parameters are set - cerr<<"Number of made clusters : "<SetName(cname); - carray.GetTree()->Write(); + //setup TPCDigitsArray + if(GetDigitsArray()) delete GetDigitsArray(); - savedir->cd(); //switch back to the input file - -} // end of function + AliTPCDigitsArray *arr = new AliTPCDigitsArray; + arr->SetClass("AliSimDigits"); + arr->Setup(fTPCParam); + arr->MakeTree(fLoader->TreeS()); -//_________________________________________________________________ -void AliTPC::Hits2ExactClustersSector(Int_t isec) -{ - //-------------------------------------------------------- - //calculate exact cross point of track and given pad row - //resulting values are expressed in "digit" coordinata - //-------------------------------------------------------- + SetDigitsArray(arr); - //----------------------------------------------------------------- - // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de - //----------------------------------------------------------------- - // - if (fClustersArray==0){ - return; - } - // - TParticle *particle; // pointer to a given particle - AliTPChit *tpcHit; // pointer to a sigle TPC hit - // Int_t sector,nhits; - Int_t ipart; - const Int_t kcmaxhits=30000; - AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4); - AliTPCFastVector & 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 - //--------------------------------------------------------------- + // 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 - TTree *tH = gAlice->TreeH(); - Stat_t ntracks = tH->GetEntries(); - Int_t npart = gAlice->GetNtrack(); - //MI change - TBranch * branch=0; - if (fHitType>1) branch = tH->GetBranch("TPC2"); - else branch = tH->GetBranch("TPC"); + for(Int_t iSector = 0; iSector < kNS; iSector++) { + + Int_t nRows = fTPCParam->GetNRow(iSector); + Int_t nDDLs = 0, indexDDL = 0; + if (iSector < kNIS) { + nDDLs = 2; + indexDDL = iSector * 2; + } + else { + nDDLs = 4; + indexDDL = (iSector-kNIS) * 4 + kNIS * 2; + } - //------------------------------------------------------------ - // Loop over tracks - //------------------------------------------------------------ + // 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); + } - for(Int_t track=0;trackGetEntry(track); // get next track - // - // Get number of the TPC hits and a pointer - // to the particles - // Loop over hits - // - Int_t currentIndex=0; - Int_t lastrow=-1; //last writen row + // Begin loop over altro data + while (input.Next()) { - //M.I. changes + if (input.GetSector() != iSector) + AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSector,input.GetSector())); - tpcHit = (AliTPChit*)FirstHit(-1); - while(tpcHit){ - - Int_t sector=tpcHit->fSector; // sector number - if(sector != isec){ - tpcHit = (AliTPChit*) NextHit(); - continue; - } + 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); - ipart=tpcHit->Track(); - if (ipartParticle(ipart); + 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)); - //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++; + 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 + + // Now fill the digits array + if (fDigitsArray->GetTree()==0) { + AliFatal("Tree not set in fDigitsArray"); + } + + 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); + } } - 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 + fDigitsArray->StoreRow(iSector,iRow); + Int_t ndig = dig->GetDigitSize(); - - 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); + AliDebug(10, + Form("*** Sector, row, compressed digits %d %d %d ***\n", + iSector,iRow,ndig)); + + fDigitsArray->ClearRow(iSector,iRow); + + } // end of the sector digitization + + for (Int_t iRow = 0; iRow < nRows; iRow++) + delete [] allBins[iRow]; + + delete [] allBins; + } - xxxx->Delete(); - -} + fLoader->WriteSDigits("OVERWRITE"); + + if(GetDigitsArray()) delete GetDigitsArray(); + SetDigitsArray(0x0); + return kTRUE; +} +//______________________________________________________________________ +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 GenerNoise(500000); //create teble with noise - char sname[100]; - char dname[100]; - sprintf(sname,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber); - sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); //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) { + AliError("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) { + AliError("Can not find branch with segments in TreeS."); + return; + } + + sdb->SetAddress(&dummy); + Stat_t nentries = t->GetEntries(); // set zero suppression @@ -1367,9 +1049,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 @@ -1377,10 +1059,11 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber) t->GetEvent(n); Int_t sec, row; if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) { - cerr<<"AliTPC warning: invalid segment ID ! "<CreateRow(sec,row); Int_t nrows = digrow->GetNRows(); Int_t ncols = digrow->GetNCols(); @@ -1396,23 +1079,16 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber) Short_t * pamp1 = digrow->GetDigits(); Int_t * ptracks1 = digrow->GetTracks(); Int_t nelems =nrows*ncols; - Int_t saturation = fTPCParam->GetADCSat(); + 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; - //if (ptracks0[0]==0) - // ptracks1[0]=1; - //else + ptracks1[0]=ptracks0[0]; ptracks1[nelems]=ptracks0[nelems]; ptracks1[2*nelems]=ptracks0[2*nelems]; @@ -1425,172 +1101,89 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber) 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 = fTPCParam->GetADCSat(); - digrow->SetDigitFast((Short_t)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); - - } - - delete digarr; - arr->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"; - // 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"); +// } + param = (AliTPCParamSR*)AliTPCcalibDB::Instance()->GetParameters(); + if (!param->IsGeoRead()){ + // + // read transformation matrices for gGeoManager + // + param->ReadGeoMatrices(); + } + if(!param){ + AliFatal("No TPC parameters found"); + } + + AliTPCPRF2D * prfinner = 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()); - rf->Update(); + + + //AliTPCRF1D * rf = new AliTPCRF1D(kTRUE); + //rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.); + //rf->SetOffset(3*param->GetZSigma()); + //rf->Update(); + // + // Use gamma 4 + // + char strgamma4[1000]; + sprintf(strgamma4,"AliTPCRF1D::Gamma4((x-0.135+%f)*%f,55,160)",3*param->GetZSigma(), 1000000000*param->GetTSample()/param->GetZWidth()); + TF1 * fgamma4 = new TF1("fgamma4",strgamma4, -1,1); + AliTPCRF1D * rf = new AliTPCRF1D(kTRUE,1000); + rf->SetParam(fgamma4,param->GetZWidth(), 1,0.2); + rf->SetOffset(3*param->GetZSigma()); + rf->Update(); + TDirectory *savedir=gDirectory; TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root"); - if (!f->IsOpen()) { - cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ; - exit(3); - } + 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()); + 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(); @@ -1608,12 +1201,53 @@ void AliTPC::SetDefaults(){ } //__________________________________________________________________ +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(); + + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + //PH runLoader->GetEvent(iEvent); + 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); + SetActiveSectors(); + 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; + } + } if(fDefaults == 0) SetDefaults(); // check if the parameters are set GenerNoise(500000); //create teble with noise @@ -1621,33 +1255,36 @@ void AliTPC::Hits2Digits(Int_t eventnumber) //setup TPCDigitsArray if(GetDigitsArray()) delete GetDigitsArray(); - + 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 -- normal digits...\n"; - - for(Int_t isec=0;isecGetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec); - - // write results - - char treeName[100]; - - sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); + 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); - GetDigitsArray()->GetTree()->SetName(treeName); - GetDigitsArray()->GetTree()->Write(); - - } - - //__________________________________________________________________ void AliTPC::Hits2SDigits2(Int_t eventnumber) { @@ -1656,24 +1293,40 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber) // summable digits - 16 bit "ADC", no noise, no saturation //----------------------------------------------------------- - //---------------------------------------------------- - // Loop over all sectors for a single event - //---------------------------------------------------- + //---------------------------------------------------- + // Loop over all sectors for a single event + //---------------------------------------------------- + + AliRunLoader* rl = fLoader->GetRunLoader(); + + rl->GetEvent(eventnumber); + if (fLoader->TreeH() == 0x0) { + if(fLoader->LoadHits()) { + AliError("Can not load hits."); + return; + } + } + SetTreeAddress(); + 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); - SetDigitsArray(arr); + arr->MakeTree(fLoader->TreeS()); - cerr<<"Digitizing TPC -- summable digits...\n"; + SetDigitsArray(arr); fDigitsSwitch=1; // summable digits @@ -1681,68 +1334,57 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber) fTPCParam->SetZeroSup(0); - for(Int_t isec=0;isecGetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec); - - - // write results - - char treeName[100]; + for(Int_t isec=0;isecGetNSector();isec++) + if (IsSectorActive(isec)) { + Hits2DigitsSector(isec); + } - sprintf(treeName,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber); - - GetDigitsArray()->GetTree()->SetName(treeName); - GetDigitsArray()->GetTree()->Write(); + fLoader->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() { //----------------------------------------------------------- // summable digits - 16 bit "ADC", no noise, no saturation //----------------------------------------------------------- + if (0) fDebugStreamer = new TTreeSRedirector("TPCSimdebug.root"); - //---------------------------------------------------- - // Loop over all sectors for a single event - //---------------------------------------------------- - //MI change - for pp run - Int_t eventnumber = gAlice->GetEvNumber(); - - 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); - SetDigitsArray(arr); - - cerr<<"Digitizing TPC -- summable digits...\n"; - - // fDigitsSwitch=1; // summable digits -for the moment direct - - for(Int_t isec=0;isecGetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec); - - - // write results - char treeName[100]; - - sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber); + if (!fTPCParam->IsGeoRead()){ + // + // read transformation matrices for gGeoManager + // + fTPCParam->ReadGeoMatrices(); + } - GetDigitsArray()->GetTree()->SetName(treeName); - GetDigitsArray()->GetTree()->Write(); + fLoader->LoadHits("read"); + fLoader->LoadSDigits("recreate"); + AliRunLoader* runLoader = fLoader->GetRunLoader(); + + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + SetTreeAddress(); + SetActiveSectors(); + Hits2SDigits2(iEvent); + } + fLoader->UnloadHits(); + fLoader->UnloadSDigits(); + if (fDebugStreamer) { + delete fDebugStreamer; + fDebugStreamer=0; + } } - - //_____________________________________________________________________________ + void AliTPC::Hits2DigitsSector(Int_t isec) { //------------------------------------------------------------------- @@ -1762,60 +1404,63 @@ void AliTPC::Hits2DigitsSector(Int_t isec) if(fDefaults == 0) SetDefaults(); - TTree *tH = gAlice->TreeH(); // pointer to the hits tree + TTree *tH = fLoader->TreeH(); // pointer to the hits tree + if (tH == 0x0) { + AliFatal("Can not find TreeH in folder"); + return; + } + Stat_t ntracks = tH->GetEntries(); - if( ntracks > 0){ - //------------------------------------------- - // Only if there are any tracks... - //------------------------------------------- TObjArray **row; - //printf("*** Processing sector number %d ***\n",isec); + Int_t nrows =fTPCParam->GetNRow(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 AliTPCFastVectors, - // each one containing electrons accepted on this - // row, assigned into tracks - //-------------------------------------------------------- + MakeSector(isec,nrows,tH,ntracks,row); - Int_t i; + //-------------------------------------------------------- + // 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 + //-------------------------------------------------------- - if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree(fDigitsFile); + Int_t i; - for (i=0;iGetTree()==0) { + AliFatal("Tree not set in fDigitsArray"); + } - AliDigits * dig = fDigitsArray->CreateRow(isec,i); + for (i=0;iCreateRow(isec,i); - DigitizeRow(i,isec,row); + DigitizeRow(i,isec,row); - fDigitsArray->StoreRow(isec,i); + fDigitsArray->StoreRow(isec,i); - Int_t ndig = dig->GetDigitSize(); - if (gDebug > 10) printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig); + Int_t ndig = dig->GetDigitSize(); + + AliDebug(10, + Form("*** Sector, row, compressed digits %d %d %d ***\n", + isec,i,ndig)); - fDigitsArray->ClearRow(isec,i); + fDigitsArray->ClearRow(isec,i); - } // end of the sector digitization + } // end of the sector digitization - for(i=0;iDelete(); - delete row[i]; - } + for(i=0;iDelete(); + delete row[i]; + } - delete [] row; // delete the array of pointers to TObjArray-s + delete [] row; // delete the array of pointers to TObjArray-s - } // ntracks >0 } // end of Hits2DigitsSector @@ -1833,9 +1478,13 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de //----------------------------------------------------------------- - Float_t zerosup = fTPCParam->GetZeroSup(); - Int_t nrows =fTPCParam->GetNRow(isec); + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); + AliTPCCalROC * gainROC = gainTPC->GetCalROC(isec); // pad gains per given sector + AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(isec); // noise per given sector + + fCurrentIndex[1]= isec; @@ -1847,10 +1496,10 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) // and a single track signal // - AliTPCFastMatrix *m1 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // integrated - AliTPCFastMatrix *m2 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // single + TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated + TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single // - AliTPCFastMatrix &total = *m1; + TMatrixF &total = *m1; // Array of pointers to the label-signal list @@ -1863,14 +1512,14 @@ 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=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); @@ -1885,24 +1534,26 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) Int_t gi=-1; Float_t fzerosup = zerosup+0.5; for(Int_t it=0;itGetValue(irow,ip); // get gain for given - pad-row pad + Float_t noisePad = noiseROC->GetValue(irow,ip); + // + q*=gain; + q+=GetNoise()*noisePad; if(q <=fzerosup) continue; // do not fill zeros q = TMath::Nint(q); - if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation + if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1; // saturation } else { - if(q <= 0.) continue; // do not fill zeros - if(q>2000.) q=2000.; - q *= 16.; - q = TMath::Nint(q); + if(q <= 0.) continue; // do not fill zeros + if(q>2000.) q=2000.; + q *= 16.; + q = TMath::Nint(q); } // @@ -1920,17 +1571,22 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) */ //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); - } - + 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 + // + // test + // + // + // glitch filters if normal simulated digits + // + if(!fDigitsSwitch) ((AliSimDigits*)dig)->GlitchFilter(); // // This row has been digitized, delete nonused stuff // @@ -1943,14 +1599,13 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows) delete m1; delete m2; - // delete m3; } // end of DigitizeRow //_____________________________________________________________________________ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, - AliTPCFastMatrix *m1, AliTPCFastMatrix *m2,Int_t *indexRange) + TMatrixF *m1, TMatrixF *m2,Int_t *indexRange) { //--------------------------------------------------------------- @@ -1964,53 +1619,53 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, // Modified: Marian Ivanov //----------------------------------------------------------------- - AliTPCFastVector *tv; + TVector *tv; - tv = (AliTPCFastVector*)p1->At(ntr); // pointer to a track - AliTPCFastVector &v = *tv; + 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 = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2; - Int_t nElectrons = (tv->GetNrows()-1)/4; + 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 - AliTPCFastMatrix &signal = *m1; - AliTPCFastMatrix &total = *m2; + TMatrixF &signal = *m1; + TMatrixF &total = *m2; // // Loop over all electrons // for(Int_t nel=0; nelGetTotalNormFac(); - Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)}; + 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)); - if (n>0) for (Int_t i =0; i=0){ - Int_t time=index[2]; - Float_t qweight = *(weight)*eltoadcfac; - - 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]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]highest){ - *(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, @@ -2120,13 +1771,46 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, // contains the track label and the position of electrons. //----------------------------------------------------------------- + // + // The trasport of the electrons through TPC drift volume + // Drift (drift velocity + velocity map(not yet implemented))) + // Application of the random processes (diffusion, gas gain) + // Systematic effects (ExB effect in drift volume + ROCs) + // + // Algorithm: + // Loop over primary electrons: + // Creation of the secondary electrons + // Loop over electrons (primary+ secondaries) + // Global coordinate frame: + // 1. Skip electrons if attached + // 2. ExB effect in drift volume + // a.) Simulation calib->GetExB()->CorrectInverse(dxyz0,dxyz1); + // b.) Reconstruction - calib->GetExB()->CorrectInverse(dxyz0,dxyz1); + // 3. Generation of gas gain (Random - Exponential distribution) + // 4. TransportElectron function (diffusion) + // + // 5. Conversion to the local coordinate frame pad-row, pad, timebin + // 6. Apply Time0 shift - AliTPCCalPad class + // a.) Plus sign in simulation + // b.) Minus sign in reconstruction + // end of loop + // //----------------------------------------------------------------- // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl + // Origin: Marian Ivanov, marian.ivanov@cern.ch //----------------------------------------------------------------- + AliTPCcalibDB* const calib=AliTPCcalibDB::Instance(); + if (gAlice){ // Set correctly the magnetic field in the ExB calculation + AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField()); + if (field) { + calib->SetExBField(field->SolenoidField()); + } + } Float_t gasgain = fTPCParam->GetGasGain(); + gasgain = gasgain/fGainFactor; Int_t i; - Float_t xyz[4]; + Float_t xyz[5]; AliTPChit *tpcHit; // pointer to a sigle TPC hit //MI change @@ -2137,14 +1821,14 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, //---------------------------------------------- // Create TObjArray-s, one for each row, - // each TObjArray will store the AliTPCFastVectors - // of electrons, one AliTPCFastVectors per each track. + // each TObjArray will store the TVectors + // of electrons, one TVectors per each track. //---------------------------------------------- - Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row - AliTPCFastVector **tracks = new AliTPCFastVector* [nrows]; //pointers to the track vectors + Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row + TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors - for(i=0; iGetEntry(track); // get next track - + //M.I. changes tpcHit = (AliTPChit*)FirstHit(-1); @@ -2184,134 +1868,193 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, continue; } - currentTrack = tpcHit->Track(); // track number + // Remove hits which arrive before the TPC opening gate signal + if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z())) + /fTPCParam->GetDriftV()+tpcHit->Time())GetGateDelay()) { + tpcHit = (AliTPChit*) NextHit(); + continue; + } + currentTrack = tpcHit->Track(); // track number - if(currentTrack != previousTrack){ + if(currentTrack != previousTrack){ - // store already filled fTrack + // store already filled fTrack - for(i=0;i0){ - 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 AliTPCFastVector - tracks[i]=0; - } - } - - nofElectrons[i]=0; - tracks[i] = new AliTPCFastVector(481); // AliTPCFastVectors for the next fTrack - - } // end of loop over rows + 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; + } + } + + nofElectrons[i]=0; + tracks[i] = new TVector(601); // TVectors for the next fTrack + + } // end of loop over rows - previousTrack=currentTrack; // update track label - } + previousTrack=currentTrack; // update track label + } - Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons) + Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons) - //--------------------------------------------------- - // Calculate the electron attachment probability - //--------------------------------------------------- + //--------------------------------------------------- + // Calculate the electron attachment probability + //--------------------------------------------------- - Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z())) - /fTPCParam->GetDriftV(); - // in microseconds! - Float_t attProb = fTPCParam->GetAttCoef()* - fTPCParam->GetOxyCont()*time; // fraction! + Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-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; + //----------------------------------------------- + // Loop over electrons + //----------------------------------------------- + Int_t index[3]; + index[1]=isec; + for(Int_t nel=0;nelRndm(0)) < attProb) continue; // electron lost! + + // + // ExB effect + // + Double_t dxyz0[3],dxyz1[3]; + dxyz0[0]=tpcHit->X(); + dxyz0[1]=tpcHit->Y(); + dxyz0[2]=tpcHit->Z(); + if (calib->GetExB()){ + calib->GetExB()->CorrectInverse(dxyz0,dxyz1); + }else{ + AliError("Not valid ExB calibration"); + dxyz1[0]=tpcHit->X(); + dxyz1[1]=tpcHit->Y(); + dxyz1[2]=tpcHit->Z(); + } + xyz[0]=dxyz1[0]; + xyz[1]=dxyz1[1]; + xyz[2]=dxyz1[2]; + // + // + // + // protection for the nonphysical avalanche size (10**6 maximum) + // + Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22); + xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); + index[0]=1; - TransportElectron(xyz,index); //MI change -august - Int_t rowNumber; - fTPCParam->GetPadRow(xyz,index); //MI change august - rowNumber = index[2]; - //transform position to local digit coordinates - //relative to nearest pad row - if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue; - 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); + TransportElectron(xyz,index); + Int_t rowNumber; + Int_t padrow = fTPCParam->GetPadRow(xyz,index); + // + // Add Time0 correction due unisochronity + // xyz[0] - pad row coordinate + // xyz[1] - pad coordinate + // xyz[2] - is in now time bin coordinate system + Float_t correction =0; + if (calib->GetPadTime0()){ + if (!calib->GetPadTime0()->GetCalROC(isec)) continue; + Int_t npads = fTPCParam->GetNPads(isec,padrow); + // Int_t pad = TMath::Nint(xyz[1]+fTPCParam->GetNPads(isec,TMath::Nint(xyz[0]))*0.5); + // pad numbering from -npads/2 .. npads/2-1 + Int_t pad = TMath::Nint(xyz[1]+npads/2); + if (pad<0) pad=0; + if (pad>=npads) pad=npads-1; + correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(padrow,pad); + // printf("%d\t%d\t%d\t%f\n",isec,padrow,pad,correction); + if (fDebugStreamer){ + (*fDebugStreamer)<<"Time0"<< + "isec="<GetNTBinsL1(); // adding Level 1 time bin offset + // + // Electron track time (for pileup simulation) + xyz[2]+=tpcHit->Time()/fTPCParam->GetTSample(); // adding time of flight + xyz[4] =0; - x1=TMath::Abs(x1); - if (y1-0.5 120){ - Int_t range = tracks[rowNumber]->GetNrows(); - if((nofElectrons[rowNumber])>(range-1)/4){ - - tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons - } + // + // 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 } - - AliTPCFastVector &v = *tracks[rowNumber]; - Int_t idx = 4*nofElectrons[rowNumber]-3; - Real_t * position = &(((AliTPCFastVector&)v).UncheckedAt(idx)); //make code faster - memcpy(position,xyz,4*sizeof(Float_t)); - - } // end of loop over 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 - tpcHit = (AliTPChit*)NextHit(); - - } // end of loop over hits - } // end of loop over tracks + tpcHit = (AliTPChit*)NextHit(); + + } // end of loop over hits + } // end of loop over tracks // // store remaining track (the last one) if not empty // - - for(i=0;i0){ - 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]; - tracks[i]=0; - } - } - - delete [] tracks; - delete [] nofElectrons; - + + 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 @@ -2322,42 +2065,9 @@ void AliTPC::Init() // // Initialise TPC detector after definition of geometry // - Int_t i; - // - if(fDebug) { - printf("\n%s: ",ClassName()); - for(i=0;i<35;i++) printf("*"); - printf(" TPC_INIT "); - for(i=0;i<35;i++) printf("*"); - printf("\n%s: ",ClassName()); - // - for(i=0;i<80;i++) printf("*"); - printf("\n"); - } + AliDebug(1,"*********************************************"); } -//_____________________________________________________________________________ -void AliTPC::MakeBranch(Option_t* option, const char *file) -{ - // - // Create Tree branches for the TPC. - // - Int_t buffersize = 4000; - char branchname[10]; - sprintf(branchname,"%s",GetName()); - - AliDetector::MakeBranch(option,file); - - const char *d = strstr(option,"D"); - - if (fDigits && gAlice->TreeD() && d) { - MakeBranchInTree(gAlice->TreeD(), - branchname, &fDigits, buffersize, file); - } - - if (fHitType>1) MakeBranch2(option,file); // MI change 14.09.2000 -} - //_____________________________________________________________________________ void AliTPC::ResetDigits() { @@ -2368,79 +2078,7 @@ void AliTPC::ResetDigits() if (fDigits) fDigits->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; -} - -//_____________________________________________________________________________ -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) @@ -2465,25 +2103,6 @@ void AliTPC::SetSide(Float_t side=0.) fSide = side; -} -//____________________________________________________________________________ -void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1, - Float_t p2,Float_t p3) -{ - - // gax mixture definition - - fNoComp = nc; - - fMixtComp[0]=c1; - fMixtComp[1]=c2; - fMixtComp[2]=c3; - - fMixtProp[0]=p1; - fMixtProp[1]=p2; - fMixtProp[2]=p3; - - } //_____________________________________________________________________________ @@ -2498,7 +2117,7 @@ void AliTPC::TransportElectron(Float_t *xyz, Int_t *index) // xyz and index must be already transformed to system 1 // - fTPCParam->Transform1to2(xyz,index); + fTPCParam->Transform1to2(xyz,index); // mis-alignment applied in this step //add diffusion Float_t driftl=xyz[2]; @@ -2513,37 +2132,35 @@ 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) +ClassImp(AliTPChit) + //______________________________________________________________________ + AliTPChit::AliTPChit() + :AliHit(), + fSector(0), + fPadRow(0), + fQ(0), + fTime(0) { // - // Creates a TPC digit object + // default // - fSector = digits[0]; - fPadRow = digits[1]; - fPad = digits[2]; - fTime = digits[3]; - fSignal = digits[4]; -} - -ClassImp(AliTPChit) - +} //_____________________________________________________________________________ -AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits): -AliHit(shunt,track) +AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits) + :AliHit(shunt,track), + fSector(0), + fPadRow(0), + fQ(0), + fTime(0) { // // Creates a TPC hit object @@ -2554,18 +2171,19 @@ AliHit(shunt,track) fY = hits[1]; fZ = hits[2]; fQ = hits[3]; + fTime = hits[4]; } - //________________________________________________________________________ // Additional code because of the AliTPCTrackHitsV2 -void AliTPC::MakeBranch2(Option_t *option,const char *file) +void AliTPC::MakeBranch(Option_t *option) { // // 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()); @@ -2573,65 +2191,26 @@ void AliTPC::MakeBranch2(Option_t *option,const char *file) // Get the pointer to the header const char *cH = strstr(option,"H"); // - if (fTrackHits && gAlice->TreeH() && cH && fHitType&4) { - // AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHitsV2",&fTrackHits, - // gAlice->TreeH(),fBufferSize,99); - //TBranch * branch = - gAlice->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits, - 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 (fTrackHitsOld && gAlice->TreeH() && cH && fHitType&2) { - AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld, - gAlice->TreeH(),fBufferSize,99); - //TBranch * branch = gAlice->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits, - // 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,&fTrackHitsOld,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 && fLoader->TreeH() && cH && fHitType&4) { + AliDebug(1,"Making branch for Type 4 Hits"); + fLoader->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99); + } + +// if (fTrackHitsOld && fLoader->TreeH() && cH && fHitType&2) { +// AliDebug(1,"Making branch for Type 2 Hits"); +// AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld, +// fLoader->TreeH(),fBufferSize,99); +// fLoader->TreeH()->GetListOfBranches()->Add(branch); +// } } void AliTPC::SetTreeAddress() { - if (fHitType<=1) AliDetector::SetTreeAddress(); + //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(); } @@ -2640,57 +2219,67 @@ void AliTPC::SetTreeAddress2() // // Set branch address for the TrackHits Tree // + AliDebug(1,""); + TBranch *branch; char branchname[20]; sprintf(branchname,"%s2",GetName()); // // Branch address for hit tree - TTree *treeH = gAlice->TreeH(); + TTree *treeH = fLoader->TreeH(); if ((treeH)&&(fHitType&4)) { branch = treeH->GetBranch(branchname); - if (branch) branch->SetAddress(&fTrackHits); - } - if ((treeH)&&(fHitType&2)) { - branch = treeH->GetBranch(branchname); - if (branch) branch->SetAddress(&fTrackHitsOld); + 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)"); +// } } void AliTPC::FinishPrimary() { if (fTrackHits &&fHitType&4) fTrackHits->FlushHitStack(); - if (fTrackHitsOld && fHitType&2) fTrackHitsOld->FlushHitStack(); + // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->FlushHitStack(); } void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits) { // - // add hit to the list + // 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 && fHitType&4) fTrackHits->AddHitKartez(vol[0],rtrack, hits[0], - 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]); + 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]); } void AliTPC::ResetHits() -{ +{ if (fHitType&1) AliDetector::ResetHits(); if (fHitType>1) ResetHits2(); } @@ -2700,7 +2289,7 @@ void AliTPC::ResetHits2() // //reset hits if (fTrackHits && fHitType&4) fTrackHits->Clear(); - if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear(); + // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear(); } @@ -2711,6 +2300,9 @@ AliHit* AliTPC::FirstHit(Int_t track) } AliHit* AliTPC::NextHit() { + // + // gets next hit + // if (fHitType>1) return NextHit2(); return AliDetector::NextHit(); @@ -2726,18 +2318,18 @@ AliHit* AliTPC::FirstHit2(Int_t track) // track is returned // if(track>=0) { - gAlice->ResetHits(); - gAlice->TreeH()->GetEvent(track); + gAlice->GetMCApp()->ResetHits(); + fLoader->TreeH()->GetEvent(track); } // if (fTrackHits && fHitType&4) { fTrackHits->First(); return fTrackHits->GetHit(); } - if (fTrackHitsOld && fHitType&2) { - fTrackHitsOld->First(); - return fTrackHitsOld->GetHit(); - } + // if (fTrackHitsOld && fHitType&2) { +// fTrackHitsOld->First(); +// return fTrackHitsOld->GetHit(); +// } else return 0; } @@ -2748,10 +2340,10 @@ AliHit* AliTPC::NextHit2() //Return the next hit for the current track - if (fTrackHitsOld && fHitType&2) { - fTrackHitsOld->Next(); - return fTrackHitsOld->GetHit(); - } +// if (fTrackHitsOld && fHitType&2) { +// fTrackHitsOld->Next(); +// return fTrackHitsOld->GetHit(); +// } if (fTrackHits) { fTrackHits->Next(); return fTrackHits->GetHit(); @@ -2760,37 +2352,26 @@ AliHit* AliTPC::NextHit2() return 0; } -void AliTPC::LoadPoints(Int_t) -{ - // - Int_t a = 0; - /* if(fHitType==1) return AliDetector::LoadPoints(a); - LoadPoints2(a); - */ - if(fHitType==1) AliDetector::LoadPoints(a); - else LoadPoints2(a); - - // LoadPoints3(a); - -} - - void AliTPC::RemapTrackHitIDs(Int_t *map) { + // + // 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 (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->fTrackID = map[info->fTrackID]; + info->SetTrackID(map[info->GetTrackID()]); } } } @@ -2802,24 +2383,24 @@ Bool_t AliTPC::TrackInVolume(Int_t id,Int_t track) //return true if current track is in volume // // return kTRUE; - if (fTrackHitsOld && fHitType&2) { - TBranch * br = gAlice->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 (fTrackHitsOld && fHitType&2) { +// TBranch * br = fLoader->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 = gAlice->TreeH()->GetBranch("fVolumes"); - TBranch * br2 = gAlice->TreeH()->GetBranch("fNVolumes"); + TBranch * br1 = fLoader->TreeH()->GetBranch("fVolumes"); + TBranch * br2 = fLoader->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); + AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes)); return kFALSE; } for (Int_t j=0;jTreeH()->GetBranch("fSector"); + TBranch * br = fLoader->TreeH()->GetBranch("fSector"); br->GetEvent(track); for (Int_t j=0;jGetEntriesFast();j++){ if ( ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE; @@ -2837,325 +2418,65 @@ Bool_t AliTPC::TrackInVolume(Int_t id,Int_t track) } -//_____________________________________________________________________________ -void AliTPC::LoadPoints2(Int_t) -{ - // - // Store x, y, z of all hits in memory - // - if (fTrackHits == 0 && fTrackHitsOld==0) return; - // - 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(); - 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;iGetTrack(); - 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(); - } - - - // - 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; +AliLoader* AliTPC::MakeLoader(const char* topfoldername) +{ + //Makes TPC loader + fLoader = new AliTPCLoader(GetName(),topfoldername); + return fLoader; } - -//_____________________________________________________________________________ -void AliTPC::LoadPoints3(Int_t) -{ - // - // Store x, y, z of all hits in memory - // - only intersection point with pad row - if (fTrackHits == 0) return; - // - Int_t nhits = fTrackHits->GetEntriesFast(); - if (nhits == 0) return; - Int_t tracks = gAlice->GetNtrack(); - if (fPoints == 0) fPoints = new TObjArray(2*tracks); - fPoints->Expand(2*tracks); - AliHit *ahit; - // - Int_t *ntrk=new Int_t[tracks]; - Int_t *limi=new Int_t[tracks]; - Float_t **coor=new Float_t*[tracks]; - for(Int_t i=0;iUncheckedAt(hit); - trk=ahit->GetTrack(); - Float_t x[3]={ahit->X(),ahit->Y(),ahit->Z()}; - Int_t index[3]={1,((AliTPChit*)ahit)->fSector,0}; - Int_t currentrow = fTPCParam->GetPadRow(x,index) ; - if (currentrow!=lastrow){ - lastrow = currentrow; - //later calculate intersection point - if(ntrk[trk]==limi[trk]) { - // - // Initialise a new track - fp=new Float_t[3*(limi[trk]+chunk)]; - if(coor[trk]) { - memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]); - delete [] coor[trk]; - } - limi[trk]+=chunk; - coor[trk] = fp; - } else { - fp = coor[trk]; - } - fp[3*ntrk[trk] ] = ahit->X(); - fp[3*ntrk[trk]+1] = ahit->Y(); - fp[3*ntrk[trk]+2] = ahit->Z(); - ntrk[trk]++; +//////////////////////////////////////////////////////////////////////// +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; + paramTPC = AliTPCcalibDB::Instance()->GetParameters(); + if (!paramTPC->IsGeoRead()){ + // + // read transformation matrices for gGeoManager + // + paramTPC->ReadGeoMatrices(); } - ahit = NextHit2(); - } - // - for(trk=0; trkSetMarkerColor(GetMarkerColor()+1); - points->SetMarkerStyle(5); - points->SetMarkerSize(0.2); - points->SetDetector(this); - points->SetParticle(trk); - // points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20); - points->SetPolyMarker(ntrk[trk],coor[trk],30); - fPoints->AddAt(points,tracks+trk); - delete [] coor[trk]; - coor[trk]=0; - } } - delete [] coor; - delete [] ntrk; - delete [] limi; -} - - - -void AliTPC::FindTrackHitsIntersection(TClonesArray * arr) -{ - - // - //fill clones array with intersection of current point with the - //middle of the row - Int_t sector; - Int_t ipart; - - const Int_t kcmaxhits=30000; - AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4); - AliTPCFastVector & xxx = *xxxx; - Int_t maxhits = kcmaxhits; - - // - AliTPChit * tpcHit=0; - tpcHit = (AliTPChit*)FirstHit2(-1); - Int_t currentIndex=0; - Int_t lastrow=-1; //last writen row + return paramTPC; - 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(); +// 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 "<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"); - } - - TStopwatch timer; - cout<<"AliTPC::Digits2Reco - determination of rec points begins"<cd(); - for(Int_t iev=firstevent;ievClose(); - cwd->cd(); - - -}