Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /// \class AliTPCSpaceCharge3D
17 : /// \brief The class calculates the space point distortions due to an arbitrary space charge distribution in 3D.
18 : ///
19 : /// The method of calculation is based on the analytical solution for the Poisson
20 : /// problem in 3D (cylindrical coordinates). The solution is used in form of
21 : /// look up tables, where the pre calculated solutions for different voxel
22 : /// positions are stored. These voxel solutions can be summed up according
23 : /// to the weight of the position of the applied space charge distribution.
24 : /// Further details can be found in \cite[chap.5]{PhD-thesis_S.Rossegger}.
25 : ///
26 : /// The class also allows a simple scaling of the resulting distortions
27 : /// via the function SetCorrectionFactor. This speeds up the calculation
28 : /// time under the assumption, that the distortions scales linearly with the
29 : /// magnitude of the space charge distribution $\rho(r,z)$ and the shape stays
30 : /// the same at higher luminosities.
31 : ///
32 : /// In contrast to the implementation in 2D (see the class AliTPCSpaceCharge above),
33 : /// the input charge distribution can be of arbitrary character. An example on how
34 : /// to produce a corresponding charge distribution can be found in the function
35 : /// WriteChargeDistributionToFile. In there, a $\rho(r,z) = (A-B\,z)/r^2$,
36 : /// with slightly different magnitude on the A and C side (due to the muon absorber),
37 : /// is superpositioned with a few leaking wires at arbitrary positions.
38 : ///
39 : /// 
40 : ///
41 : /// \author Stefan Rossegger
42 : /// \date 19/06/2010
43 :
44 :
45 : #include "AliMagF.h"
46 : #include "TGeoGlobalMagField.h"
47 : #include "AliTPCcalibDB.h"
48 : #include "AliTPCParam.h"
49 : #include "AliLog.h"
50 : #include "TH2F.h"
51 : #include "TH3F.h"
52 : #include "TFile.h"
53 : #include "TVector.h"
54 : #include "TMatrix.h"
55 : #include "TMatrixD.h"
56 : #include "AliTPCCorrection.h"
57 : #include "TMath.h"
58 : #include "AliTPCROC.h"
59 : #include "AliTPCSpaceCharge2D2D.h"
60 :
61 : /// \cond CLASSIMP
62 24 : ClassImp(AliTPCSpaceCharge2D2D)
63 : /// \endcond
64 :
65 : AliTPCSpaceCharge2D2D::AliTPCSpaceCharge2D2D()
66 0 : : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
67 0 : fC0(0.),fC1(0.),
68 0 : fCorrectionFactor(1.),
69 0 : fInitLookUp(kFALSE),
70 0 : fSCDataFileName(""),
71 0 : fSCLookUpPOCsFileName3D(""),
72 0 : fSCLookUpPOCsFileNameRZ(""),
73 0 : fSCLookUpPOCsFileNameRPhi(""),
74 0 : fSCdensityInRZ(0),
75 0 : fSCdensityInRPhiA(0),
76 0 : fSCdensityInRPhiC(0)
77 0 : {
78 : //
79 : // default constructor
80 : //
81 :
82 : // Array which will contain the solution according to the setted charge density distribution
83 : // see InitSpaceCharge3DDistortion() function
84 0 : for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
85 0 : fLookUpErOverEz[k] = new TMatrixF(kNR,kNZ);
86 0 : fLookUpEphiOverEz[k] = new TMatrixF(kNR,kNZ);
87 0 : fLookUpDeltaEz[k] = new TMatrixF(kNR,kNZ);
88 0 : fSCdensityDistribution[k] = new TMatrixF(kNR,kNZ);
89 : }
90 0 : fSCdensityInRZ = new TMatrixD(kNR,kNZ);
91 0 : fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
92 0 : fSCdensityInRPhiC = new TMatrixD(kNR,kNPhi);
93 :
94 : // location of the precalculated look up tables
95 :
96 0 : fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
97 0 : fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
98 0 : fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
99 : // fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
100 :
101 :
102 :
103 : // standard location of the space charge distibution ... can be changes
104 0 : fSCDataFileName="$(ALICE_ROOT)/TPC/TPCcalib/maps/sc_3D_distribution_Sim.root";
105 :
106 : // SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
107 :
108 :
109 0 : }
110 :
111 0 : AliTPCSpaceCharge2D2D::~AliTPCSpaceCharge2D2D() {
112 : /// default destructor
113 :
114 0 : for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
115 0 : delete fLookUpErOverEz[k];
116 0 : delete fLookUpEphiOverEz[k];
117 0 : delete fLookUpDeltaEz[k];
118 0 : delete fSCdensityDistribution[k];
119 : }
120 0 : delete fSCdensityInRZ;
121 0 : delete fSCdensityInRPhiA;
122 0 : delete fSCdensityInRPhiC;
123 :
124 0 : }
125 :
126 :
127 : void AliTPCSpaceCharge2D2D::Init() {
128 : /// Initialization funtion
129 :
130 0 : AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
131 0 : if (!magF) AliError("Magneticd field - not initialized");
132 0 : Double_t bzField = magF->SolenoidField()/10.; //field in T
133 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
134 0 : if (!param) AliError("Parameters - not initialized");
135 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
136 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
137 0 : Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
138 : // Correction Terms for effective omegaTau; obtained by a laser calibration run
139 0 : SetOmegaTauT1T2(wt,fT1,fT2);
140 :
141 0 : InitSpaceCharge3DDistortion(); // fill the look up table
142 0 : }
143 :
144 : void AliTPCSpaceCharge2D2D::Update(const TTimeStamp &/*timeStamp*/) {
145 : /// Update function
146 :
147 0 : AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
148 0 : if (!magF) AliError("Magneticd field - not initialized");
149 0 : Double_t bzField = magF->SolenoidField()/10.; //field in T
150 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
151 0 : if (!param) AliError("Parameters - not initialized");
152 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
153 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
154 0 : Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
155 : // Correction Terms for effective omegaTau; obtained by a laser calibration run
156 0 : SetOmegaTauT1T2(wt,fT1,fT2);
157 :
158 : // SetCorrectionFactor(1.); // should come from some database
159 :
160 0 : }
161 :
162 :
163 : void AliTPCSpaceCharge2D2D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
164 : /// Calculates the correction due the Space Charge effect within the TPC drift volume
165 :
166 0 : if (!fInitLookUp) {
167 0 : AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
168 0 : InitSpaceCharge3DDistortion();
169 0 : }
170 :
171 : Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
172 :
173 : Float_t intEr, intEphi, intdEz ;
174 : Double_t r, phi, z ;
175 : Int_t sign;
176 :
177 0 : r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
178 0 : phi = TMath::ATan2(x[1],x[0]) ;
179 0 : if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
180 0 : z = x[2] ; // Create temporary copy of x[2]
181 :
182 0 : if ( (roc%36) < 18 ) {
183 : sign = 1; // (TPC A side)
184 0 : } else {
185 : sign = -1; // (TPC C side)
186 : }
187 :
188 0 : if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
189 0 : if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
190 :
191 :
192 0 : if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
193 0 : AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
194 :
195 : // Get the Er and Ephi field integrals plus the integral over DeltaEz
196 0 : intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
197 0 : fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
198 0 : intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
199 0 : fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
200 0 : intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
201 0 : fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
202 :
203 : // Calculate distorted position
204 0 : if ( r > 0.0 ) {
205 0 : phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
206 0 : r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
207 0 : }
208 0 : Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
209 :
210 : // Calculate correction in cartesian coordinates
211 0 : dx[0] = - (r * TMath::Cos(phi) - x[0]);
212 0 : dx[1] = - (r * TMath::Sin(phi) - x[1]);
213 0 : dx[2] = - dz; // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
214 :
215 0 : }
216 :
217 : void AliTPCSpaceCharge2D2D::InitSpaceCharge3DDistortion() {
218 : /// Initialization of the Lookup table which contains the solutions of the
219 : /// "space charge" (poisson) problem - Faster and more accureate
220 : ///
221 : /// Method: Weighted sum-up of the different fields within the look up table
222 : /// but using two lookup tables with higher granularity in the (r,z) and the (rphi)- plane to emulate
223 : /// more realistic space charges. (r,z) from primary ionisation. (rphi) for possible Gating leaks
224 :
225 0 : if (fInitLookUp) {
226 0 : AliInfo("Lookup table was already initialized! Doing it again anyway ...");
227 : // return;
228 0 : }
229 :
230 : // ------------------------------------------------------------------------------------------------------
231 : // step 1: lookup table in rz, fine grid, radial symetric, to emulate primary ionization
232 :
233 0 : AliInfo("Step 1: Preparation of the weighted look-up tables.");
234 :
235 : // lookup table in rz, fine grid
236 :
237 0 : TFile *fZR = new TFile(fSCLookUpPOCsFileNameRZ.Data(),"READ");
238 0 : if ( !fZR ) {
239 0 : AliError("Precalculated POC-looup-table in ZR could not be found");
240 0 : return;
241 : }
242 :
243 : // units are in [m]
244 0 : TVector *gridf1 = (TVector*) fZR->Get("constants");
245 : TVector &grid1 = *gridf1;
246 0 : TMatrix *coordf1 = (TMatrix*) fZR->Get("coordinates");
247 : TMatrix &coord1 = *coordf1;
248 0 : TMatrix *coordPOCf1 = (TMatrix*) fZR->Get("POCcoord");
249 : TMatrix &coordPOC1 = *coordPOCf1;
250 :
251 0 : Int_t rows = (Int_t)grid1(0); // number of points in r direction - from RZ or RPhi table
252 0 : Int_t phiSlices = (Int_t)grid1(1); // number of points in phi - from RPhi table
253 0 : Int_t columns = (Int_t)grid1(2); // number of points in z direction - from RZ table
254 :
255 0 : Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
256 0 : Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
257 :
258 : // temporary matrices needed for the calculation // for rotational symmetric RZ table, phislices is 1
259 :
260 0 : TMatrixD *arrayofErA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
261 0 : TMatrixD *arrayofErC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
262 :
263 0 : TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
264 0 : TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
265 :
266 :
267 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
268 :
269 0 : arrayofErA[k] = new TMatrixD(rows,columns) ;
270 0 : arrayofdEzA[k] = new TMatrixD(rows,columns) ;
271 0 : arrayofErC[k] = new TMatrixD(rows,columns) ;
272 0 : arrayofdEzC[k] = new TMatrixD(rows,columns) ;
273 :
274 0 : arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
275 0 : arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
276 0 : arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
277 0 : arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
278 :
279 : // zero initialization not necessary, it is done in the constructor of TMatrix
280 :
281 : }
282 :
283 : // list of points as used during sum up
284 0 : Double_t rlist1[kNRows], zedlist1[kNColumns];// , philist1[phiSlices];
285 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
286 0 : rlist1[i] = fgkIFCRadius + i*gridSizeR ;
287 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
288 0 : zedlist1[j] = j * gridSizeZ ;
289 : }
290 : }
291 :
292 0 : TTree *treePOC = (TTree*)fZR->Get("POCall");
293 :
294 0 : TVector *bEr = 0; //TVector *bEphi= 0;
295 0 : TVector *bEz = 0;
296 :
297 0 : treePOC->SetBranchAddress("Er",&bEr);
298 0 : treePOC->SetBranchAddress("Ez",&bEz);
299 :
300 :
301 : // Read the complete tree and do a weighted sum-up over the POC configurations
302 : // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
303 :
304 0 : Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
305 : Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
306 :
307 0 : for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
308 :
309 0 : treePOC->GetEntry(itreepC);
310 :
311 : // center of the POC voxel in [meter]
312 0 : Double_t r0 = coordPOC1(ipC,0);
313 0 : Double_t phi0 = coordPOC1(ipC,1);
314 0 : Double_t z0 = coordPOC1(ipC,2);
315 :
316 0 : ipC++; // POC configuration counter
317 :
318 : // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
319 : // note: coordinates are in [cm]
320 0 : Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100, 1); // partial load in r,z
321 0 : Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100, 1); // partial load in r,z
322 :
323 : // Summing up the vector components according to their weight
324 :
325 : Int_t ip = 0;
326 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
327 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
328 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
329 :
330 : // check wether the coordinates were screwed
331 0 : if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
332 0 : TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
333 0 : AliError("internal error: coordinate system was screwed during the sum-up");
334 0 : printf("sum-up: (r,z)=(%f,%f)\n",rlist1[i],zedlist1[j]);
335 0 : printf("lookup: (r,z)=(%f,%f)\n",coord1(0,ip)*100,coord1(2,ip)*100);
336 0 : AliError("Don't trust the results of the space charge calculation!");
337 0 : }
338 :
339 : // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
340 : // This will be the most frequent usage (hopefully)
341 : // That's why we have to do this here ...
342 :
343 0 : TMatrixD &erA = *arrayofErA[k] ;
344 0 : TMatrixD &dEzA = *arrayofdEzA[k] ;
345 :
346 0 : TMatrixD &erC = *arrayofErC[k] ;
347 0 : TMatrixD &dEzC = *arrayofdEzC[k] ;
348 :
349 : // Sum up - Efield values in [V/m] -> transition to [V/cm]
350 0 : erA(i,j) += ((*bEr)(ip)) * weightA /100;
351 0 : erC(i,j) += ((*bEr)(ip)) * weightC /100;
352 0 : dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
353 0 : dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
354 :
355 : // increase the counter
356 0 : ip++;
357 : }
358 : }
359 : } // end coordinate loop
360 : } // end POC loop
361 :
362 :
363 : // -------------------------------------------------------------------------------
364 : // Division by the Ez (drift) field and integration along z
365 :
366 : // AliInfo("Step 1: Division and integration");
367 :
368 0 : Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
369 :
370 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
371 :
372 : // matrices holding the solution - summation of POC charges // see above
373 0 : TMatrixD &erA = *arrayofErA[k] ;
374 0 : TMatrixD &dezA = *arrayofdEzA[k] ;
375 0 : TMatrixD &erC = *arrayofErC[k] ;
376 0 : TMatrixD &dezC = *arrayofdEzC[k] ;
377 :
378 : // matrices which will contain the integrated fields (divided by the drift field)
379 0 : TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
380 0 : TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
381 0 : TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
382 0 : TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
383 :
384 0 : for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
385 0 : for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
386 : // Count backwards to facilitate integration over Z
387 :
388 : Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
389 : // by trapezoidal rule.
390 :
391 0 : erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
392 0 : deltaEzA(i,j) = 0;
393 0 : erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
394 0 : deltaEzC(i,j) = 0;
395 :
396 0 : for ( Int_t m = j ; m < columns ; m++ ) { // integration
397 :
398 0 : erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
399 0 : erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
400 0 : deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
401 0 : deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
402 :
403 0 : if ( index != 4 ) index = 4; else index = 2 ;
404 :
405 : }
406 :
407 0 : if ( index == 4 ) {
408 0 : erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
409 0 : erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
410 0 : deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
411 0 : deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
412 0 : }
413 0 : if ( index == 2 ) {
414 0 : erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
415 0 : erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
416 0 : deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
417 0 : deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
418 0 : }
419 0 : if ( j == columns-2 ) {
420 0 : erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
421 0 : erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
422 0 : deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
423 0 : deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
424 0 : }
425 0 : if ( j == columns-1 ) {
426 0 : erOverEzA(i,j) = 0;
427 0 : erOverEzC(i,j) = 0;
428 0 : deltaEzA(i,j) = 0;
429 0 : deltaEzC(i,j) = 0;
430 0 : }
431 : }
432 : }
433 :
434 : }
435 :
436 : // AliInfo("Step 1: Interpolation to Standard grid");
437 :
438 : // -------------------------------------------------------------------------------
439 : // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
440 :
441 : const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
442 :
443 : Double_t r, z;//phi, z ;
444 0 : for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
445 : // phi = fgkPhiList[k] ;
446 :
447 : // final lookup table
448 0 : TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
449 0 : TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
450 :
451 : // calculated and integrated tables - just one phi slice
452 0 : TMatrixD &erOverEzA = *arrayofEroverEzA[0] ;
453 0 : TMatrixD &deltaEzA = *arrayofDeltaEzA[0];
454 0 : TMatrixD &erOverEzC = *arrayofEroverEzC[0] ;
455 0 : TMatrixD &deltaEzC = *arrayofDeltaEzC[0];
456 :
457 :
458 0 : for ( Int_t j = 0 ; j < kNZ ; j++ ) {
459 :
460 0 : z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
461 :
462 0 : for ( Int_t i = 0 ; i < kNR ; i++ ) {
463 0 : r = fgkRList[i] ;
464 :
465 : // Interpolate Lookup tables onto standard grid
466 0 : if (fgkZList[j]>0) {
467 0 : erOverEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzA );
468 0 : deltaEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzA );
469 0 : } else {
470 0 : erOverEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzC );
471 0 : deltaEzFinal(i,j) = - Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzC );
472 : // negative coordinate system on C side
473 : }
474 :
475 : } // end r loop
476 : } // end z loop
477 : } // end phi loop
478 :
479 :
480 : // clear the temporary arrays lists
481 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
482 :
483 0 : delete arrayofErA[k];
484 0 : delete arrayofdEzA[k];
485 0 : delete arrayofErC[k];
486 0 : delete arrayofdEzC[k];
487 :
488 0 : delete arrayofEroverEzA[k];
489 0 : delete arrayofDeltaEzA[k];
490 0 : delete arrayofEroverEzC[k];
491 0 : delete arrayofDeltaEzC[k];
492 :
493 : }
494 :
495 0 : fZR->Close();
496 :
497 : // ------------------------------------------------------------------------------------------------------
498 : // Step 2: Load and sum up lookup table in rphi, fine grid, to emulate for example a GG leak
499 :
500 : // AliInfo("Step 2: Preparation of the weighted look-up table");
501 :
502 0 : TFile *fRPhi = new TFile(fSCLookUpPOCsFileNameRPhi.Data(),"READ");
503 0 : if ( !fRPhi ) {
504 0 : AliError("Precalculated POC-looup-table in RPhi could not be found");
505 0 : return;
506 : }
507 :
508 : // units are in [m]
509 0 : TVector *gridf2 = (TVector*) fRPhi->Get("constants");
510 : TVector &grid2 = *gridf2;
511 0 : TMatrix *coordf2 = (TMatrix*) fRPhi->Get("coordinates");
512 : TMatrix &coord2 = *coordf2;
513 0 : TMatrix *coordPOCf2 = (TMatrix*) fRPhi->Get("POCcoord");
514 : TMatrix &coordPOC2 = *coordPOCf2;
515 :
516 0 : rows = (Int_t)grid2(0); // number of points in r direction
517 0 : phiSlices = (Int_t)grid2(1); // number of points in phi
518 0 : columns = (Int_t)grid2(2); // number of points in z direction
519 :
520 0 : gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
521 0 : Float_t gridSizePhi = TMath::TwoPi()/phiSlices; // unit in [rad]
522 0 : gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
523 :
524 : // list of points as used during sum up
525 0 : Double_t rlist2[kNRows], philist2[kNPhiSlices], zedlist2[kNColumns];
526 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
527 0 : philist2[k] = gridSizePhi * k;
528 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
529 0 : rlist2[i] = fgkIFCRadius + i*gridSizeR ;
530 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
531 0 : zedlist2[j] = j * gridSizeZ ;
532 : }
533 : }
534 : } // only done once
535 :
536 : // temporary matrices needed for the calculation
537 :
538 0 : TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
539 0 : TMatrixD *arrayofErC2[kNPhiSlices], *arrayofEphiC2[kNPhiSlices], *arrayofdEzC2[kNPhiSlices];
540 :
541 0 : TMatrixD *arrayofEroverEzA2[kNPhiSlices], *arrayofEphioverEzA2[kNPhiSlices], *arrayofDeltaEzA2[kNPhiSlices];
542 0 : TMatrixD *arrayofEroverEzC2[kNPhiSlices], *arrayofEphioverEzC2[kNPhiSlices], *arrayofDeltaEzC2[kNPhiSlices];
543 :
544 :
545 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
546 :
547 0 : arrayofErA2[k] = new TMatrixD(rows,columns) ;
548 0 : arrayofEphiA2[k] = new TMatrixD(rows,columns) ;
549 0 : arrayofdEzA2[k] = new TMatrixD(rows,columns) ;
550 0 : arrayofErC2[k] = new TMatrixD(rows,columns) ;
551 0 : arrayofEphiC2[k] = new TMatrixD(rows,columns) ;
552 0 : arrayofdEzC2[k] = new TMatrixD(rows,columns) ;
553 :
554 0 : arrayofEroverEzA2[k] = new TMatrixD(rows,columns) ;
555 0 : arrayofEphioverEzA2[k] = new TMatrixD(rows,columns) ;
556 0 : arrayofDeltaEzA2[k] = new TMatrixD(rows,columns) ;
557 0 : arrayofEroverEzC2[k] = new TMatrixD(rows,columns) ;
558 0 : arrayofEphioverEzC2[k] = new TMatrixD(rows,columns) ;
559 0 : arrayofDeltaEzC2[k] = new TMatrixD(rows,columns) ;
560 :
561 : // zero initialization not necessary, it is done in the constructor of TMatrix
562 :
563 : }
564 :
565 :
566 0 : treePOC = (TTree*)fRPhi->Get("POCall");
567 :
568 : // TVector *bEr = 0; // done above
569 0 : TVector *bEphi= 0;
570 : // TVector *bEz = 0; // done above
571 :
572 0 : treePOC->SetBranchAddress("Er",&bEr);
573 0 : treePOC->SetBranchAddress("Ephi",&bEphi);
574 0 : treePOC->SetBranchAddress("Ez",&bEz);
575 :
576 : // Read the complete tree and do a weighted sum-up over the POC configurations
577 : // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
578 :
579 0 : treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
580 : ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
581 :
582 0 : for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
583 :
584 0 : treePOC->GetEntry(itreepC);
585 :
586 : // center of the POC voxel in [meter]
587 0 : Double_t r0 = coordPOC2(ipC,0);
588 0 : Double_t phi0 = coordPOC2(ipC,1);
589 : // Double_t z0 = coordPOC2(ipC,2);
590 :
591 : // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
592 : // note: coordinates are in [cm]
593 0 : Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, 0.499, 2); // partial load in r,phi
594 0 : Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-0.499, 2); // partial load in r,phi
595 :
596 : // printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
597 :
598 : // Summing up the vector components according to their weight
599 :
600 : Int_t ip = 0;
601 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
602 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
603 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
604 :
605 : // check wether the coordinates were screwed
606 0 : if (TMath::Abs((coord2(0,ip)*100-rlist2[i]))>1 ||
607 0 : TMath::Abs((coord2(1,ip)-philist2[k]))>1 ||
608 0 : TMath::Abs((coord2(2,ip)*100-zedlist2[j]))>1) {
609 0 : AliError("internal error: coordinate system was screwed during the sum-up");
610 0 : printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
611 0 : printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist2[i],philist2[k],zedlist2[j]);
612 0 : AliError("Don't trust the results of the space charge calculation!");
613 0 : }
614 :
615 : // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
616 : // This will be the most frequent usage (hopefully)
617 : // That's why we have to do this here ...
618 :
619 0 : TMatrixD &erA = *arrayofErA2[k] ;
620 0 : TMatrixD &ephiA = *arrayofEphiA2[k];
621 0 : TMatrixD &dEzA = *arrayofdEzA2[k] ;
622 :
623 0 : TMatrixD &erC = *arrayofErC2[k] ;
624 0 : TMatrixD &ephiC = *arrayofEphiC2[k];
625 0 : TMatrixD &dEzC = *arrayofdEzC2[k] ;
626 :
627 : // Sum up - Efield values in [V/m] -> transition to [V/cm]
628 0 : erA(i,j) += ((*bEr)(ip)) * weightA /100;
629 0 : erC(i,j) += ((*bEr)(ip)) * weightC /100;
630 0 : ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
631 0 : ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
632 0 : dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
633 0 : dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
634 :
635 : // increase the counter
636 0 : ip++;
637 : }
638 : }
639 : } // end coordinate loop
640 :
641 :
642 : // Rotation and summation in the rest of the dPhiSteps
643 : // which were not stored in the this tree due to storage & symmetry reasons
644 :
645 :
646 0 : Int_t phiPoints = (Int_t) grid2(1);
647 0 : Int_t phiPOC = (Int_t) grid2(4);
648 :
649 : // printf("%d %d\n",phiPOC,flagRadSym);
650 :
651 0 : for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
652 :
653 0 : Double_t phi0R = phiiC*phi0*2 + phi0; // rotate further
654 :
655 : // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
656 : // note: coordinates are in [cm] // ecxept z
657 0 : weightA = GetSpaceChargeDensity(r0*100,phi0R, 0.499, 2); // partial load in r,phi
658 0 : weightC = GetSpaceChargeDensity(r0*100,phi0R,-0.499, 2); // partial load in r,phi
659 :
660 : // printf("%f %f : %e %e\n",r0,phi0R,weightA,weightC);
661 :
662 : // Summing up the vector components according to their weight
663 : ip = 0;
664 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
665 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
666 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
667 :
668 : // Note: rotating the coordinated during the sum up
669 :
670 0 : Int_t rotVal = (phiPoints/phiPOC)*phiiC;
671 : Int_t ipR = -1;
672 :
673 0 : if ((ip%phiPoints)>=rotVal) {
674 0 : ipR = ip-rotVal;
675 0 : } else {
676 0 : ipR = ip+(phiPoints-rotVal);
677 : }
678 :
679 : // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
680 : // This will be the most frequent usage
681 : // That's why we have to do this here and not outside the loop ...
682 :
683 0 : TMatrixD &erA = *arrayofErA2[k] ;
684 0 : TMatrixD &ephiA = *arrayofEphiA2[k];
685 0 : TMatrixD &dEzA = *arrayofdEzA2[k] ;
686 :
687 0 : TMatrixD &erC = *arrayofErC2[k] ;
688 0 : TMatrixD &ephiC = *arrayofEphiC2[k];
689 0 : TMatrixD &dEzC = *arrayofdEzC2[k] ;
690 :
691 : // Sum up - Efield values in [V/m] -> transition to [V/cm]
692 0 : erA(i,j) += ((*bEr)(ipR)) * weightA /100;
693 0 : erC(i,j) += ((*bEr)(ipR)) * weightC /100;
694 0 : ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
695 0 : ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
696 0 : dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
697 0 : dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
698 :
699 : // increase the counter
700 0 : ip++;
701 : }
702 : }
703 : } // end coordinate loop
704 :
705 : } // end phi-POC summation (phiiC)
706 :
707 0 : ipC++; // POC configuration counter
708 :
709 : // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
710 :
711 : }
712 :
713 :
714 :
715 :
716 : // -------------------------------------------------------------------------------
717 : // Division by the Ez (drift) field and integration along z
718 :
719 : // AliInfo("Step 2: Division and integration");
720 :
721 :
722 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
723 :
724 : // matrices holding the solution - summation of POC charges // see above
725 0 : TMatrixD &erA = *arrayofErA2[k] ;
726 0 : TMatrixD &ephiA = *arrayofEphiA2[k];
727 0 : TMatrixD &dezA = *arrayofdEzA2[k] ;
728 0 : TMatrixD &erC = *arrayofErC2[k] ;
729 0 : TMatrixD &ephiC = *arrayofEphiC2[k];
730 0 : TMatrixD &dezC = *arrayofdEzC2[k] ;
731 :
732 : // matrices which will contain the integrated fields (divided by the drift field)
733 0 : TMatrixD &erOverEzA = *arrayofEroverEzA2[k] ;
734 0 : TMatrixD &ephiOverEzA = *arrayofEphioverEzA2[k];
735 0 : TMatrixD &deltaEzA = *arrayofDeltaEzA2[k];
736 0 : TMatrixD &erOverEzC = *arrayofEroverEzC2[k] ;
737 0 : TMatrixD &ephiOverEzC = *arrayofEphioverEzC2[k];
738 0 : TMatrixD &deltaEzC = *arrayofDeltaEzC2[k];
739 :
740 0 : for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
741 0 : for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
742 : // Count backwards to facilitate integration over Z
743 :
744 : Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
745 :
746 0 : erOverEzA(i,j) = 0;
747 0 : ephiOverEzA(i,j) = 0;
748 0 : deltaEzA(i,j) = 0;
749 0 : erOverEzC(i,j) = 0;
750 0 : ephiOverEzC(i,j) = 0;
751 0 : deltaEzC(i,j) = 0;
752 :
753 0 : for ( Int_t m = j ; m < columns ; m++ ) { // integration
754 :
755 0 : erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
756 0 : erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
757 0 : ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
758 0 : ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
759 0 : deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
760 0 : deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
761 :
762 0 : if ( index != 4 ) index = 4; else index = 2 ;
763 :
764 : }
765 :
766 0 : if ( index == 4 ) {
767 0 : erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
768 0 : erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
769 0 : ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
770 0 : ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
771 0 : deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
772 0 : deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
773 0 : }
774 0 : if ( index == 2 ) {
775 0 : erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
776 0 : erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
777 0 : ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
778 0 : ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
779 0 : deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
780 0 : deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
781 0 : }
782 0 : if ( j == columns-2 ) {
783 0 : erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
784 0 : erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
785 0 : ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
786 0 : ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
787 0 : deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
788 0 : deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
789 0 : }
790 0 : if ( j == columns-1 ) {
791 0 : erOverEzA(i,j) = 0;
792 0 : erOverEzC(i,j) = 0;
793 0 : ephiOverEzA(i,j) = 0;
794 0 : ephiOverEzC(i,j) = 0;
795 0 : deltaEzA(i,j) = 0;
796 0 : deltaEzC(i,j) = 0;
797 0 : }
798 : }
799 : }
800 :
801 : }
802 :
803 0 : AliInfo("Step 2: Interpolation to Standard grid");
804 :
805 : // -------------------------------------------------------------------------------
806 : // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
807 :
808 :
809 0 : for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
810 0 : Double_t phi = fgkPhiList[k] ;
811 :
812 : // final lookup table
813 0 : TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
814 0 : TMatrixF &ephiOverEzFinal = *fLookUpEphiOverEz[k];
815 0 : TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
816 :
817 0 : for ( Int_t j = 0 ; j < kNZ ; j++ ) {
818 :
819 0 : z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
820 :
821 0 : for ( Int_t i = 0 ; i < kNR ; i++ ) {
822 0 : r = fgkRList[i] ;
823 :
824 : // Interpolate Lookup tables onto standard grid
825 0 : if (fgkZList[j]>0) {
826 0 : erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
827 0 : rlist2, zedlist2, philist2, arrayofEroverEzA2 );
828 0 : ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
829 0 : rlist2, zedlist2, philist2, arrayofEphioverEzA2);
830 0 : deltaEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
831 0 : rlist2, zedlist2, philist2, arrayofDeltaEzA2 );
832 0 : } else {
833 0 : erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
834 0 : rlist2, zedlist2, philist2, arrayofEroverEzC2 );
835 0 : ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
836 0 : rlist2, zedlist2, philist2, arrayofEphioverEzC2);
837 0 : deltaEzFinal(i,j) -= Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
838 0 : rlist2, zedlist2, philist2, arrayofDeltaEzC2 );
839 : }
840 :
841 : } // end r loop
842 : } // end z loop
843 : } // end phi loop
844 :
845 :
846 : // clear the temporary arrays lists
847 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
848 :
849 0 : delete arrayofErA2[k];
850 0 : delete arrayofEphiA2[k];
851 0 : delete arrayofdEzA2[k];
852 0 : delete arrayofErC2[k];
853 0 : delete arrayofEphiC2[k];
854 0 : delete arrayofdEzC2[k];
855 :
856 0 : delete arrayofEroverEzA2[k];
857 0 : delete arrayofEphioverEzA2[k];
858 0 : delete arrayofDeltaEzA2[k];
859 0 : delete arrayofEroverEzC2[k];
860 0 : delete arrayofEphioverEzC2[k];
861 0 : delete arrayofDeltaEzC2[k];
862 :
863 : }
864 :
865 0 : fRPhi->Close();
866 :
867 : // FINISHED
868 :
869 0 : fInitLookUp = kTRUE;
870 :
871 0 : }
872 :
873 : void AliTPCSpaceCharge2D2D::InitSpaceCharge3DDistortionCourse() {
874 : /// Initialization of the Lookup table which contains the solutions of the
875 : /// "space charge" (poisson) problem
876 : ///
877 : /// The sum-up uses a look-up table which contains different discretized Space charge fields
878 : /// in order to calculate the corresponding field deviations due to a given (discretized)
879 : /// space charge distribution ....
880 : ///
881 : /// Method of calculation: Weighted sum-up of the different fields within the look up table
882 : /// Note: Full 3d version: Course and slow ...
883 :
884 0 : if (fInitLookUp) {
885 0 : AliInfo("Lookup table was already initialized!");
886 : // return;
887 0 : }
888 :
889 0 : AliInfo("Preparation of the weighted look-up table");
890 :
891 0 : TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
892 0 : if ( !f ) {
893 0 : AliError("Precalculated POC-looup-table could not be found");
894 0 : return;
895 : }
896 :
897 : // units are in [m]
898 0 : TVector *gridf = (TVector*) f->Get("constants");
899 : TVector &grid = *gridf;
900 0 : TMatrix *coordf = (TMatrix*) f->Get("coordinates");
901 : TMatrix &coord = *coordf;
902 0 : TMatrix *coordPOCf = (TMatrix*) f->Get("POCcoord");
903 : TMatrix &coordPOC = *coordPOCf;
904 :
905 : Bool_t flagRadSym = 0;
906 0 : if (grid(1)==1 && grid(4)==1) {
907 : // AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
908 : flagRadSym=1;
909 0 : }
910 :
911 0 : Int_t rows = (Int_t)grid(0); // number of points in r direction
912 0 : Int_t phiSlices = (Int_t)grid(1); // number of points in phi
913 0 : Int_t columns = (Int_t)grid(2); // number of points in z direction
914 :
915 0 : const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
916 0 : const Float_t gridSizePhi = TMath::TwoPi()/phiSlices; // unit in [rad]
917 0 : const Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
918 :
919 : // temporary matrices needed for the calculation
920 0 : TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
921 0 : TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
922 :
923 0 : TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
924 0 : TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
925 :
926 :
927 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
928 :
929 0 : arrayofErA[k] = new TMatrixD(rows,columns) ;
930 0 : arrayofEphiA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
931 0 : arrayofdEzA[k] = new TMatrixD(rows,columns) ;
932 0 : arrayofErC[k] = new TMatrixD(rows,columns) ;
933 0 : arrayofEphiC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
934 0 : arrayofdEzC[k] = new TMatrixD(rows,columns) ;
935 :
936 0 : arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
937 0 : arrayofEphioverEzA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
938 0 : arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
939 0 : arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
940 0 : arrayofEphioverEzC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
941 0 : arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
942 :
943 : // Set the values to zero the lookup tables
944 : // not necessary, it is done in the constructor of TMatrix - code deleted
945 :
946 : }
947 :
948 : // list of points as used in the interpolation (during sum up)
949 0 : Double_t rlist[kNRows], zedlist[kNColumns] , philist[kNPhiSlices];
950 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
951 0 : philist[k] = gridSizePhi * k;
952 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
953 0 : rlist[i] = fgkIFCRadius + i*gridSizeR ;
954 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
955 0 : zedlist[j] = j * gridSizeZ ;
956 : }
957 : }
958 : } // only done once
959 :
960 :
961 0 : TTree *treePOC = (TTree*)f->Get("POCall");
962 :
963 0 : TVector *bEr = 0; TVector *bEphi= 0; TVector *bEz = 0;
964 :
965 0 : treePOC->SetBranchAddress("Er",&bEr);
966 0 : if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
967 0 : treePOC->SetBranchAddress("Ez",&bEz);
968 :
969 :
970 : // Read the complete tree and do a weighted sum-up over the POC configurations
971 : // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
972 :
973 0 : Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
974 : Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
975 :
976 0 : for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
977 :
978 0 : treePOC->GetEntry(itreepC);
979 :
980 : // center of the POC voxel in [meter]
981 0 : Double_t r0 = coordPOC(ipC,0);
982 0 : Double_t phi0 = coordPOC(ipC,1);
983 0 : Double_t z0 = coordPOC(ipC,2);
984 :
985 0 : ipC++; // POC configuration counter
986 :
987 : // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
988 : // note: coordinates are in [cm]
989 0 : Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
990 0 : Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
991 :
992 : // Summing up the vector components according to their weight
993 :
994 : Int_t ip = 0;
995 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
996 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
997 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
998 :
999 : // check wether the coordinates were screwed
1000 0 : if (TMath::Abs((coord(0,ip)*100-rlist[i]))>1 ||
1001 0 : TMath::Abs((coord(1,ip)-philist[k]))>1 ||
1002 0 : TMath::Abs((coord(2,ip)*100-zedlist[j]))>1) {
1003 0 : AliError("internal error: coordinate system was screwed during the sum-up");
1004 0 : printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
1005 0 : printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist[i],philist[k],zedlist[j]);
1006 0 : AliError("Don't trust the results of the space charge calculation!");
1007 0 : }
1008 :
1009 : // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1010 : // This will be the most frequent usage (hopefully)
1011 : // That's why we have to do this here ...
1012 :
1013 0 : TMatrixD &erA = *arrayofErA[k] ;
1014 0 : TMatrixD &ephiA = *arrayofEphiA[k];
1015 0 : TMatrixD &dEzA = *arrayofdEzA[k] ;
1016 :
1017 0 : TMatrixD &erC = *arrayofErC[k] ;
1018 0 : TMatrixD &ephiC = *arrayofEphiC[k];
1019 0 : TMatrixD &dEzC = *arrayofdEzC[k] ;
1020 :
1021 : // Sum up - Efield values in [V/m] -> transition to [V/cm]
1022 0 : erA(i,j) += ((*bEr)(ip)) * weightA /100;
1023 0 : erC(i,j) += ((*bEr)(ip)) * weightC /100;
1024 0 : if (!flagRadSym) {
1025 0 : ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
1026 0 : ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
1027 0 : }
1028 0 : dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
1029 0 : dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
1030 :
1031 : // increase the counter
1032 0 : ip++;
1033 : }
1034 : }
1035 : } // end coordinate loop
1036 :
1037 :
1038 : // Rotation and summation in the rest of the dPhiSteps
1039 : // which were not stored in the this tree due to storage & symmetry reasons
1040 :
1041 0 : Int_t phiPoints = (Int_t) grid(1);
1042 0 : Int_t phiPOC = (Int_t) grid(4);
1043 :
1044 : // printf("%d %d\n",phiPOC,flagRadSym);
1045 :
1046 0 : for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
1047 :
1048 0 : r0 = coordPOC(ipC,0);
1049 0 : phi0 = coordPOC(ipC,1);
1050 0 : z0 = coordPOC(ipC,2);
1051 :
1052 0 : ipC++; // POC conf. counter
1053 :
1054 : // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1055 : // note: coordinates are in [cm]
1056 0 : weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
1057 0 : weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
1058 :
1059 : // printf("%f %f %f: %e %e\n",r0,phi0,z0,weightA,weightC);
1060 :
1061 : // Summing up the vector components according to their weight
1062 : ip = 0;
1063 0 : for ( Int_t j = 0 ; j < columns ; j++ ) {
1064 0 : for ( Int_t i = 0 ; i < rows ; i++ ) {
1065 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1066 :
1067 : // Note: rotating the coordinated during the sum up
1068 :
1069 0 : Int_t rotVal = (phiPoints/phiPOC)*phiiC;
1070 : Int_t ipR = -1;
1071 :
1072 0 : if ((ip%phiPoints)>=rotVal) {
1073 0 : ipR = ip-rotVal;
1074 0 : } else {
1075 0 : ipR = ip+(phiPoints-rotVal);
1076 : }
1077 :
1078 : // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1079 : // This will be the most frequent usage
1080 : // That's why we have to do this here and not outside the loop ...
1081 :
1082 0 : TMatrixD &erA = *arrayofErA[k] ;
1083 0 : TMatrixD &ephiA = *arrayofEphiA[k];
1084 0 : TMatrixD &dEzA = *arrayofdEzA[k] ;
1085 :
1086 0 : TMatrixD &erC = *arrayofErC[k] ;
1087 0 : TMatrixD &ephiC = *arrayofEphiC[k];
1088 0 : TMatrixD &dEzC = *arrayofdEzC[k] ;
1089 :
1090 : // Sum up - Efield values in [V/m] -> transition to [V/cm]
1091 0 : erA(i,j) += ((*bEr)(ipR)) * weightA /100;
1092 0 : erC(i,j) += ((*bEr)(ipR)) * weightC /100;
1093 0 : if (!flagRadSym) {
1094 0 : ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
1095 0 : ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
1096 0 : }
1097 0 : dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
1098 0 : dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
1099 :
1100 : // increase the counter
1101 0 : ip++;
1102 : }
1103 : }
1104 : } // end coordinate loop
1105 :
1106 : } // end phi-POC summation (phiiC)
1107 :
1108 :
1109 : // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
1110 :
1111 : }
1112 :
1113 :
1114 :
1115 : // -------------------------------------------------------------------------------
1116 : // Division by the Ez (drift) field and integration along z
1117 :
1118 0 : AliInfo("Division and integration");
1119 :
1120 0 : Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
1121 :
1122 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
1123 :
1124 : // matrices holding the solution - summation of POC charges // see above
1125 0 : TMatrixD &erA = *arrayofErA[k] ;
1126 0 : TMatrixD &ephiA = *arrayofEphiA[k];
1127 0 : TMatrixD &dezA = *arrayofdEzA[k] ;
1128 0 : TMatrixD &erC = *arrayofErC[k] ;
1129 0 : TMatrixD &ephiC = *arrayofEphiC[k];
1130 0 : TMatrixD &dezC = *arrayofdEzC[k] ;
1131 :
1132 : // matrices which will contain the integrated fields (divided by the drift field)
1133 0 : TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
1134 0 : TMatrixD &ephiOverEzA = *arrayofEphioverEzA[k];
1135 0 : TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
1136 0 : TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
1137 0 : TMatrixD &ephiOverEzC = *arrayofEphioverEzC[k];
1138 0 : TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
1139 :
1140 0 : for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
1141 0 : for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
1142 : // Count backwards to facilitate integration over Z
1143 :
1144 : Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
1145 :
1146 0 : erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0; deltaEzA(i,j) = 0;
1147 0 : erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0; deltaEzC(i,j) = 0;
1148 :
1149 0 : for ( Int_t m = j ; m < columns ; m++ ) { // integration
1150 :
1151 0 : erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
1152 0 : erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
1153 0 : if (!flagRadSym) {
1154 0 : ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
1155 0 : ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
1156 0 : }
1157 0 : deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
1158 0 : deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
1159 :
1160 0 : if ( index != 4 ) index = 4; else index = 2 ;
1161 :
1162 : }
1163 :
1164 0 : if ( index == 4 ) {
1165 0 : erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
1166 0 : erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
1167 0 : if (!flagRadSym) {
1168 0 : ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
1169 0 : ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
1170 0 : }
1171 0 : deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
1172 0 : deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
1173 0 : }
1174 0 : if ( index == 2 ) {
1175 0 : erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
1176 0 : erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
1177 0 : if (!flagRadSym) {
1178 0 : ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
1179 0 : ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
1180 0 : }
1181 0 : deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
1182 0 : deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
1183 0 : }
1184 0 : if ( j == columns-2 ) {
1185 0 : erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
1186 0 : erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
1187 0 : if (!flagRadSym) {
1188 0 : ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
1189 0 : ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
1190 0 : }
1191 0 : deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
1192 0 : deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
1193 0 : }
1194 0 : if ( j == columns-1 ) {
1195 0 : erOverEzA(i,j) = 0;
1196 0 : erOverEzC(i,j) = 0;
1197 0 : if (!flagRadSym) {
1198 0 : ephiOverEzA(i,j) = 0;
1199 0 : ephiOverEzC(i,j) = 0;
1200 0 : }
1201 0 : deltaEzA(i,j) = 0;
1202 0 : deltaEzC(i,j) = 0;
1203 0 : }
1204 : }
1205 : }
1206 :
1207 : }
1208 :
1209 :
1210 :
1211 0 : AliInfo("Interpolation to Standard grid");
1212 :
1213 : // -------------------------------------------------------------------------------
1214 : // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
1215 :
1216 : const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
1217 :
1218 : Double_t r, phi, z ;
1219 0 : for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1220 0 : phi = fgkPhiList[k] ;
1221 :
1222 0 : TMatrixF &erOverEz = *fLookUpErOverEz[k] ;
1223 0 : TMatrixF &ephiOverEz = *fLookUpEphiOverEz[k];
1224 0 : TMatrixF &deltaEz = *fLookUpDeltaEz[k] ;
1225 :
1226 0 : for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1227 :
1228 0 : z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
1229 :
1230 0 : for ( Int_t i = 0 ; i < kNR ; i++ ) {
1231 0 : r = fgkRList[i] ;
1232 :
1233 : // Interpolate Lookup tables onto standard grid
1234 0 : if (fgkZList[j]>0) {
1235 0 : erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1236 0 : rlist, zedlist, philist, arrayofEroverEzA );
1237 0 : ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1238 0 : rlist, zedlist, philist, arrayofEphioverEzA);
1239 0 : deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1240 0 : rlist, zedlist, philist, arrayofDeltaEzA );
1241 0 : } else {
1242 0 : erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1243 0 : rlist, zedlist, philist, arrayofEroverEzC );
1244 0 : ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1245 0 : rlist, zedlist, philist, arrayofEphioverEzC);
1246 0 : deltaEz(i,j) = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1247 0 : rlist, zedlist, philist, arrayofDeltaEzC );
1248 : // negative coordinate system on C side
1249 : }
1250 :
1251 : } // end r loop
1252 : } // end z loop
1253 : } // end phi loop
1254 :
1255 :
1256 : // clear the temporary arrays lists
1257 0 : for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1258 :
1259 0 : delete arrayofErA[k];
1260 0 : delete arrayofEphiA[k];
1261 0 : delete arrayofdEzA[k];
1262 0 : delete arrayofErC[k];
1263 0 : delete arrayofEphiC[k];
1264 0 : delete arrayofdEzC[k];
1265 :
1266 0 : delete arrayofEroverEzA[k];
1267 0 : delete arrayofEphioverEzA[k];
1268 0 : delete arrayofDeltaEzA[k];
1269 0 : delete arrayofEroverEzC[k];
1270 0 : delete arrayofEphioverEzC[k];
1271 0 : delete arrayofDeltaEzC[k];
1272 :
1273 : }
1274 :
1275 0 : fInitLookUp = kTRUE;
1276 :
1277 0 : }
1278 :
1279 :
1280 : void AliTPCSpaceCharge2D2D::SetSCDataFileName(TString fname) {
1281 : /// Set & load the Space charge density distribution from a file
1282 : /// (linear interpolation onto a standard grid)
1283 :
1284 :
1285 0 : fSCDataFileName = fname;
1286 :
1287 0 : TFile *f = new TFile(fSCDataFileName.Data(),"READ");
1288 0 : if (!f) {
1289 0 : AliError(Form("File %s, which should contain the space charge distribution, could not be found",
1290 : fSCDataFileName.Data()));
1291 0 : return;
1292 : }
1293 :
1294 0 : TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
1295 0 : if (!densityRZ) {
1296 0 : AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1297 : fSCDataFileName.Data(),"SpaceChargeInRZ"));
1298 0 : return;
1299 : }
1300 :
1301 0 : TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
1302 0 : if (!densityRPhi) {
1303 0 : AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1304 : fSCDataFileName.Data(),"SpaceChargeInRPhi"));
1305 0 : return;
1306 : }
1307 :
1308 :
1309 : Double_t r, phi, z ;
1310 :
1311 0 : TMatrixD &scDensityInRZ = *fSCdensityInRZ;
1312 0 : TMatrixD &scDensityInRPhiA = *fSCdensityInRPhiA;
1313 0 : TMatrixD &scDensityInRPhiC = *fSCdensityInRPhiC;
1314 0 : for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1315 0 : phi = fgkPhiList[k] ;
1316 0 : TMatrixF &scDensity = *fSCdensityDistribution[k] ;
1317 0 : for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1318 0 : z = fgkZList[j] ;
1319 0 : for ( Int_t i = 0 ; i < kNR ; i++ ) {
1320 0 : r = fgkRList[i] ;
1321 :
1322 : // partial load in (r,z)
1323 0 : if (k==0) // do just once
1324 0 : scDensityInRZ(i,j) = densityRZ->Interpolate(r,z);
1325 :
1326 : // partial load in (r,phi)
1327 0 : if ( j==0 || j == kNZ/2 ) {
1328 0 : if (z>0)
1329 0 : scDensityInRPhiA(i,k) = densityRPhi->Interpolate(r,phi,0.499); // A side
1330 : else
1331 0 : scDensityInRPhiC(i,k) = densityRPhi->Interpolate(r,phi,-0.499); // C side
1332 : }
1333 :
1334 : // Full 3D configuration ...
1335 0 : if (z>0)
1336 0 : scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
1337 : else
1338 0 : scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
1339 : }
1340 : }
1341 : }
1342 :
1343 0 : f->Close();
1344 :
1345 0 : fInitLookUp = kFALSE;
1346 :
1347 :
1348 0 : }
1349 :
1350 :
1351 : Float_t AliTPCSpaceCharge2D2D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
1352 : /// returns the (input) space charge density at a given point according
1353 : /// Note: input in [cm], output in [C/m^3/e0] !!
1354 :
1355 0 : while (phi<0) phi += TMath::TwoPi();
1356 0 : while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
1357 :
1358 :
1359 : // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
1360 : Int_t order = 1; //
1361 : Float_t sc = 0;
1362 :
1363 0 : if (mode == 0) { // return full load
1364 0 : sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
1365 0 : fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
1366 :
1367 0 : } else if (mode == 1) { // return partial load in (r,z)
1368 0 : TMatrixD &scDensityInRZ = *fSCdensityInRZ;
1369 0 : sc = Interpolate2DTable(order, r, z, kNR, kNZ, fgkRList, fgkZList, scDensityInRZ );
1370 :
1371 0 : } else if (mode == 2) { // return partial load in (r,phi)
1372 :
1373 0 : if (z>0) {
1374 0 : TMatrixD &scDensityInRPhi = *fSCdensityInRPhiA;
1375 0 : sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1376 0 : } else {
1377 0 : TMatrixD &scDensityInRPhi = *fSCdensityInRPhiC;
1378 0 : sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1379 : }
1380 :
1381 : } else {
1382 : // should i give a warning?
1383 : sc = 0;
1384 : }
1385 :
1386 : // printf("%f %f %f: %f\n",r,phi,z,sc);
1387 :
1388 0 : return sc;
1389 : }
1390 :
1391 :
1392 : TH2F * AliTPCSpaceCharge2D2D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
1393 : /// return a simple histogramm containing the space charge distribution (input for the calculation)
1394 :
1395 0 : TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1396 : nx,-250.,250.,ny,-250.,250.);
1397 :
1398 0 : for (Int_t iy=1;iy<=ny;++iy) {
1399 0 : Double_t yp = h->GetYaxis()->GetBinCenter(iy);
1400 0 : for (Int_t ix=1;ix<=nx;++ix) {
1401 0 : Double_t xp = h->GetXaxis()->GetBinCenter(ix);
1402 :
1403 0 : Float_t r = TMath::Sqrt(xp*xp+yp*yp);
1404 0 : Float_t phi = TMath::ATan2(yp,xp);
1405 :
1406 0 : if (85.<=r && r<=250.) {
1407 0 : Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1408 0 : h->SetBinContent(ix,iy,sc);
1409 0 : } else {
1410 0 : h->SetBinContent(ix,iy,0.);
1411 : }
1412 : }
1413 : }
1414 :
1415 0 : return h;
1416 : }
1417 :
1418 : TH2F * AliTPCSpaceCharge2D2D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
1419 : /// return a simple histogramm containing the space charge distribution (input for the calculation)
1420 :
1421 0 : TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1422 : nz,-250.,250.,nr,85.,250.);
1423 :
1424 0 : for (Int_t ir=1;ir<=nr;++ir) {
1425 0 : Float_t r = h->GetYaxis()->GetBinCenter(ir);
1426 0 : for (Int_t iz=1;iz<=nz;++iz) {
1427 0 : Float_t z = h->GetXaxis()->GetBinCenter(iz);
1428 0 : Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1429 0 : h->SetBinContent(iz,ir,sc);
1430 : }
1431 : }
1432 :
1433 0 : return h;
1434 : }
1435 :
1436 : void AliTPCSpaceCharge2D2D::WriteChargeDistributionToFile(const char* fname) {
1437 : /// Example on how to write a Space charge distribution into a File
1438 : /// (see below: estimate from scaling STAR measurements to Alice)
1439 : /// Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1440 : /// the needed calculation time
1441 :
1442 0 : TFile *f = new TFile(fname,"RECREATE");
1443 :
1444 : // some grid, not too course
1445 : Int_t nr = 350;
1446 : Int_t nphi = 180;
1447 : Int_t nz = 500;
1448 :
1449 0 : Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1450 0 : Double_t dphi = TMath::TwoPi()/(nphi+1);
1451 0 : Double_t dz = 500./(nz+1);
1452 : Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
1453 :
1454 :
1455 : // Charge distribution in ZR (rotational symmetric) ------------------
1456 :
1457 0 : TH2F *histoZR = new TH2F("chargeZR","chargeZR",
1458 0 : nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1459 0 : nz,-250-dz-safty,250+dz+safty);
1460 :
1461 0 : for (Int_t ir=1;ir<=nr;++ir) {
1462 0 : Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1463 0 : for (Int_t iz=1;iz<=nz;++iz) {
1464 0 : Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
1465 :
1466 : // recalculation to meter
1467 : Double_t lZ = 2.5; // approx. TPC drift length
1468 0 : Double_t rpM = rp/100.; // in [m]
1469 0 : Double_t zpM = TMath::Abs(zp/100.); // in [m]
1470 :
1471 : // setting of mb multiplicity and Interaction rate
1472 : Double_t multiplicity = 950;
1473 : Double_t intRate = 7800;
1474 :
1475 : // calculation of "scaled" parameters
1476 0 : Double_t a = multiplicity*intRate/79175;
1477 0 : Double_t b = a/lZ;
1478 :
1479 0 : Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
1480 :
1481 0 : charge = charge*fgke0; // [C/m^3]
1482 :
1483 0 : if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
1484 :
1485 : // charge = 0; // for tests
1486 0 : histoZR->SetBinContent(ir,iz,charge);
1487 : }
1488 : }
1489 :
1490 0 : histoZR->Write("SpaceChargeInRZ");
1491 :
1492 :
1493 : // Charge distribution in RPhi (e.g. Floating GG wire) ------------
1494 :
1495 0 : TH3F *histoRPhi = new TH3F("chargeRPhi","chargeRPhi",
1496 0 : nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1497 0 : nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
1498 : 2,-1,1); // z part - to allow A and C side differences
1499 :
1500 : // some 'arbitrary' GG leaks
1501 : Int_t nGGleaks = 5;
1502 0 : Double_t secPosA[5] = {3,6,6,11,13}; // sector
1503 0 : Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm
1504 0 : Double_t secPosC[5] = {1,8,12,15,15}; // sector
1505 0 : Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm
1506 :
1507 0 : for (Int_t ir=1;ir<=nr;++ir) {
1508 0 : Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir);
1509 0 : for (Int_t iphi=1;iphi<=nphi;++iphi) {
1510 0 : Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi);
1511 0 : for (Int_t iz=1;iz<=2;++iz) {
1512 0 : Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz);
1513 :
1514 : Double_t charge = 0;
1515 :
1516 0 : for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks
1517 :
1518 : // A side
1519 0 : Double_t secPos = secPosA[igg];
1520 0 : Double_t radialPos = radialPosA[igg];
1521 :
1522 0 : if (zp<0) { // C side
1523 0 : secPos = secPosC[igg];
1524 0 : radialPos = radialPosC[igg];
1525 0 : }
1526 :
1527 : // some 'arbitrary' GG leaks
1528 0 : if ( (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice
1529 0 : if ( rp>(radialPos-2.5) && rp<(radialPos+2.5)) // 5 cm slice
1530 0 : charge = 300;
1531 : }
1532 :
1533 : }
1534 :
1535 0 : charge = charge*fgke0; // [C/m^3]
1536 :
1537 0 : histoRPhi->SetBinContent(ir,iphi,iz,charge);
1538 : }
1539 : }
1540 : }
1541 :
1542 0 : histoRPhi->Write("SpaceChargeInRPhi");
1543 :
1544 0 : f->Close();
1545 :
1546 0 : }
1547 :
1548 :
1549 : void AliTPCSpaceCharge2D2D::Print(const Option_t* option) const {
1550 : /// Print function to check the settings of the boundary vectors
1551 : /// option=="a" prints the C0 and C1 coefficents for calibration purposes
1552 :
1553 0 : TString opt = option; opt.ToLower();
1554 0 : printf("%s\n",GetTitle());
1555 0 : printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
1556 0 : printf(" SC correction factor: %f \n",fCorrectionFactor);
1557 :
1558 0 : if (opt.Contains("a")) { // Print all details
1559 0 : printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
1560 0 : printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
1561 : }
1562 :
1563 0 : if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
1564 :
1565 0 : }
|
