]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/TPCbase/AliTPCSpaceCharge3D.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
doxy: TPC/TPCbase converted
[u/mrichter/AliRoot.git] / TPC / TPCbase / AliTPCSpaceCharge3D.cxx
CommitLineData
cc3e558a 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
7d855b04 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 AliTPCSpaceChargeabove),
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/// Marian Ivanov change: 26.06.2013
40/// Usage of the realy 3D space charge map as an optional input
41/// SetInputSpaceCharge map.
42/// In case given map is used 2 2D maps are ignored and scaling functions $\rho(r,z) = (A-B\,z)/r^2$,
43/// will not work
92a85338 44
45
b4caed64 46// End_Html
47//
6a1caa6b 48// Begin_Macro(source)
b4caed64 49// {
50// gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
7d855b04 51// TCanvas *c2 = new TCanvas("cAliTPCSpaceCharge3D","cAliTPCSpaceCharge3D",500,400);
b4caed64 52// AliTPCSpaceCharge3D sc;
53// sc.WriteChargeDistributionToFile("SC_zr2_GGleaks.root");
54// sc.SetSCDataFileName("SC_zr2_GGleaks.root");
55// sc.SetOmegaTauT1T2(0,1,1); // B=0
56// sc.InitSpaceCharge3DDistortion();
57// sc.CreateHistoDRinXY(15,300,300)->Draw("colz");
58// return c2;
7d855b04 59// }
b4caed64 60// End_Macro
61//
62// Begin_Html
63// <p>
7d855b04 64// Date: 19/06/2010 <br>
65// Authors: Stefan Rossegger
66// End_Html
b4caed64 67// _________________________________________________________________
68
cc3e558a 69
70#include "AliMagF.h"
71#include "TGeoGlobalMagField.h"
72#include "AliTPCcalibDB.h"
73#include "AliTPCParam.h"
74#include "AliLog.h"
75#include "TH2F.h"
76#include "TH3F.h"
77#include "TFile.h"
78#include "TVector.h"
79#include "TMatrix.h"
80#include "TMatrixD.h"
81
82#include "TMath.h"
83#include "AliTPCROC.h"
84#include "AliTPCSpaceCharge3D.h"
92a85338 85#include "AliSysInfo.h"
cc3e558a 86
7d855b04 87/// \cond CLASSIMP
cc3e558a 88ClassImp(AliTPCSpaceCharge3D)
7d855b04 89/// \endcond
cc3e558a 90
91AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
92 : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
93 fC0(0.),fC1(0.),
94 fCorrectionFactor(1.),
95 fInitLookUp(kFALSE),
15687d71 96 fSCDataFileName(""),
97 fSCLookUpPOCsFileName3D(""),
98 fSCLookUpPOCsFileNameRZ(""),
99 fSCLookUpPOCsFileNameRPhi(""),
100 fSCdensityInRZ(0),
7d855b04 101 fSCdensityInRPhiA(0),
92a85338 102 fSCdensityInRPhiC(0),
103 fSpaceChargeHistogram3D(0),
104 fSpaceChargeHistogramRPhi(0),
105 fSpaceChargeHistogramRZ(0)
cc3e558a 106{
107 //
108 // default constructor
109 //
110
111 // Array which will contain the solution according to the setted charge density distribution
112 // see InitSpaceCharge3DDistortion() function
113 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
7d855b04 114 fLookUpErOverEz[k] = new TMatrixF(kNR,kNZ);
2bf29b72 115 fLookUpEphiOverEz[k] = new TMatrixF(kNR,kNZ);
7d855b04 116 fLookUpDeltaEz[k] = new TMatrixF(kNR,kNZ);
2bf29b72 117 fSCdensityDistribution[k] = new TMatrixF(kNR,kNZ);
cc3e558a 118 }
15687d71 119 fSCdensityInRZ = new TMatrixD(kNR,kNZ);
120 fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
121 fSCdensityInRPhiC = new TMatrixD(kNR,kNPhi);
122
123 // location of the precalculated look up tables
124
125 fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
126 fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/Calib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
752b0cc7 127 fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
128 // fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
7d855b04 129
15687d71 130
131
132 // standard location of the space charge distibution ... can be changes
133 fSCDataFileName="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root";
134
135 // SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
cc3e558a 136
cc3e558a 137
138}
139
140AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
7d855b04 141 /// default destructor
142
cc3e558a 143 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
144 delete fLookUpErOverEz[k];
145 delete fLookUpEphiOverEz[k];
146 delete fLookUpDeltaEz[k];
147 delete fSCdensityDistribution[k];
148 }
15687d71 149 delete fSCdensityInRZ;
150 delete fSCdensityInRPhiA;
151 delete fSCdensityInRPhiC;
92a85338 152 delete fSpaceChargeHistogram3D;
153 delete fSpaceChargeHistogramRPhi;
154 delete fSpaceChargeHistogramRZ;
cc3e558a 155}
156
157
cc3e558a 158void AliTPCSpaceCharge3D::Init() {
7d855b04 159 /// Initialization funtion
160
cc3e558a 161 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
162 if (!magF) AliError("Magneticd field - not initialized");
163 Double_t bzField = magF->SolenoidField()/10.; //field in T
164 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
165 if (!param) AliError("Parameters - not initialized");
166 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
167 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
7d855b04 168 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
cc3e558a 169 // Correction Terms for effective omegaTau; obtained by a laser calibration run
170 SetOmegaTauT1T2(wt,fT1,fT2);
171
172 InitSpaceCharge3DDistortion(); // fill the look up table
173}
174
175void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
7d855b04 176 /// Update function
177
cc3e558a 178 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
179 if (!magF) AliError("Magneticd field - not initialized");
180 Double_t bzField = magF->SolenoidField()/10.; //field in T
181 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
182 if (!param) AliError("Parameters - not initialized");
183 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
184 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
7d855b04 185 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
cc3e558a 186 // Correction Terms for effective omegaTau; obtained by a laser calibration run
187 SetOmegaTauT1T2(wt,fT1,fT2);
188
189 // SetCorrectionFactor(1.); // should come from some database
190
191}
192
193
cc3e558a 194void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
7d855b04 195 /// Calculates the correction due the Space Charge effect within the TPC drift volume
cc3e558a 196
197 if (!fInitLookUp) {
198 AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
199 InitSpaceCharge3DDistortion();
200 }
201
7d855b04 202 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
203
2bf29b72 204 Float_t intEr, intEphi, intdEz ;
cc3e558a 205 Double_t r, phi, z ;
206 Int_t sign;
207
208 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
209 phi = TMath::ATan2(x[1],x[0]) ;
210 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
211 z = x[2] ; // Create temporary copy of x[2]
212
213 if ( (roc%36) < 18 ) {
214 sign = 1; // (TPC A side)
215 } else {
216 sign = -1; // (TPC C side)
217 }
7d855b04 218
cc3e558a 219 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
220 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
7d855b04 221
cc3e558a 222
223 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
224 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
225
7d855b04 226 // Get the Er and Ephi field integrals plus the integral over DeltaEz
227 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
cc3e558a 228 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
7d855b04 229 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
cc3e558a 230 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
7d855b04 231 intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
cc3e558a 232 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
233
234 // Calculate distorted position
235 if ( r > 0.0 ) {
7d855b04 236 phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
237 r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
cc3e558a 238 }
239 Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
7d855b04 240
cc3e558a 241 // Calculate correction in cartesian coordinates
752b0cc7 242 dx[0] = - (r * TMath::Cos(phi) - x[0]);
7d855b04 243 dx[1] = - (r * TMath::Sin(phi) - x[1]);
752b0cc7 244 dx[2] = - dz; // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
cc3e558a 245
246}
247
cc3e558a 248void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
7d855b04 249 /// Initialization of the Lookup table which contains the solutions of the
250 /// "space charge" (poisson) problem - Faster and more accureate
251 ///
252 /// Method: Weighted sum-up of the different fields within the look up table
253 /// but using two lookup tables with higher granularity in the (r,z) and the (rphi)- plane to emulate
254 /// more realistic space charges. (r,z) from primary ionisation. (rphi) for possible Gating leaks
15687d71 255
256 if (fInitLookUp) {
257 AliInfo("Lookup table was already initialized! Doing it again anyway ...");
db59c7fb 258 return;
15687d71 259 }
7d855b04 260
15687d71 261 // ------------------------------------------------------------------------------------------------------
262 // step 1: lookup table in rz, fine grid, radial symetric, to emulate primary ionization
263
264 AliInfo("Step 1: Preparation of the weighted look-up tables.");
7d855b04 265
15687d71 266 // lookup table in rz, fine grid
267
268 TFile *fZR = new TFile(fSCLookUpPOCsFileNameRZ.Data(),"READ");
269 if ( !fZR ) {
270 AliError("Precalculated POC-looup-table in ZR could not be found");
271 return;
7d855b04 272 }
15687d71 273
274 // units are in [m]
7d855b04 275 TVector *gridf1 = (TVector*) fZR->Get("constants");
15687d71 276 TVector &grid1 = *gridf1;
277 TMatrix *coordf1 = (TMatrix*) fZR->Get("coordinates");
278 TMatrix &coord1 = *coordf1;
279 TMatrix *coordPOCf1 = (TMatrix*) fZR->Get("POCcoord");
280 TMatrix &coordPOC1 = *coordPOCf1;
7d855b04 281
15687d71 282 Int_t rows = (Int_t)grid1(0); // number of points in r direction - from RZ or RPhi table
283 Int_t phiSlices = (Int_t)grid1(1); // number of points in phi - from RPhi table
284 Int_t columns = (Int_t)grid1(2); // number of points in z direction - from RZ table
285
286 Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
287 Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
288
289 // temporary matrices needed for the calculation // for rotational symmetric RZ table, phislices is 1
7d855b04 290
291 TMatrixD *arrayofErA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
292 TMatrixD *arrayofErC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
15687d71 293
9f3b99e2 294 TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
295 TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
15687d71 296
7d855b04 297
15687d71 298 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 299
15687d71 300 arrayofErA[k] = new TMatrixD(rows,columns) ;
301 arrayofdEzA[k] = new TMatrixD(rows,columns) ;
302 arrayofErC[k] = new TMatrixD(rows,columns) ;
303 arrayofdEzC[k] = new TMatrixD(rows,columns) ;
304
305 arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
306 arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
307 arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
308 arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
309
7d855b04 310 // zero initialization not necessary, it is done in the constructor of TMatrix
15687d71 311
312 }
7d855b04 313
15687d71 314 // list of points as used during sum up
9f3b99e2 315 Double_t rlist1[kNRows], zedlist1[kNColumns];// , philist1[phiSlices];
15687d71 316 for ( Int_t i = 0 ; i < rows ; i++ ) {
317 rlist1[i] = fgkIFCRadius + i*gridSizeR ;
7d855b04 318 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 319 zedlist1[j] = j * gridSizeZ ;
320 }
321 }
7d855b04 322
15687d71 323 TTree *treePOC = (TTree*)fZR->Get("POCall");
324
325 TVector *bEr = 0; //TVector *bEphi= 0;
326 TVector *bEz = 0;
7d855b04 327
15687d71 328 treePOC->SetBranchAddress("Er",&bEr);
329 treePOC->SetBranchAddress("Ez",&bEz);
330
331
332 // Read the complete tree and do a weighted sum-up over the POC configurations
333 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7d855b04 334
15687d71 335 Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
336 Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
337
338 for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
7d855b04 339
15687d71 340 treePOC->GetEntry(itreepC);
341
342 // center of the POC voxel in [meter]
343 Double_t r0 = coordPOC1(ipC,0);
344 Double_t phi0 = coordPOC1(ipC,1);
345 Double_t z0 = coordPOC1(ipC,2);
346
347 ipC++; // POC configuration counter
348
349 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
350 // note: coordinates are in [cm]
351 Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100, 1); // partial load in r,z
352 Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100, 1); // partial load in r,z
7d855b04 353
15687d71 354 // Summing up the vector components according to their weight
355
356 Int_t ip = 0;
7d855b04 357 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 358 for ( Int_t i = 0 ; i < rows ; i++ ) {
359 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 360
15687d71 361 // check wether the coordinates were screwed
7d855b04 362 if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
363 TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
15687d71 364 AliError("internal error: coordinate system was screwed during the sum-up");
9f3b99e2 365 printf("sum-up: (r,z)=(%f,%f)\n",rlist1[i],zedlist1[j]);
366 printf("lookup: (r,z)=(%f,%f)\n",coord1(0,ip)*100,coord1(2,ip)*100);
15687d71 367 AliError("Don't trust the results of the space charge calculation!");
368 }
7d855b04 369
15687d71 370 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
371 // This will be the most frequent usage (hopefully)
372 // That's why we have to do this here ...
373
374 TMatrixD &erA = *arrayofErA[k] ;
375 TMatrixD &dEzA = *arrayofdEzA[k] ;
7d855b04 376
15687d71 377 TMatrixD &erC = *arrayofErC[k] ;
378 TMatrixD &dEzC = *arrayofdEzC[k] ;
7d855b04 379
15687d71 380 // Sum up - Efield values in [V/m] -> transition to [V/cm]
381 erA(i,j) += ((*bEr)(ip)) * weightA /100;
382 erC(i,j) += ((*bEr)(ip)) * weightC /100;
383 dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
384 dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
385
386 // increase the counter
387 ip++;
388 }
389 }
390 } // end coordinate loop
391 } // end POC loop
392
393
394 // -------------------------------------------------------------------------------
395 // Division by the Ez (drift) field and integration along z
396
397 // AliInfo("Step 1: Division and integration");
398
399 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
400
401 for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
402
403 // matrices holding the solution - summation of POC charges // see above
404 TMatrixD &erA = *arrayofErA[k] ;
405 TMatrixD &dezA = *arrayofdEzA[k] ;
406 TMatrixD &erC = *arrayofErC[k] ;
407 TMatrixD &dezC = *arrayofdEzC[k] ;
408
409 // matrices which will contain the integrated fields (divided by the drift field)
410 TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
411 TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
412 TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
7d855b04 413 TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
414
15687d71 415 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
7d855b04 416 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
15687d71 417 // Count backwards to facilitate integration over Z
418
7d855b04 419 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
420 // by trapezoidal rule.
15687d71 421
422 erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
423 deltaEzA(i,j) = 0;
7d855b04 424 erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
15687d71 425 deltaEzC(i,j) = 0;
426
427 for ( Int_t m = j ; m < columns ; m++ ) { // integration
428
429 erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
430 erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
431 deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
432 deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
433
434 if ( index != 4 ) index = 4; else index = 2 ;
435
436 }
437
438 if ( index == 4 ) {
439 erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
440 erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
441 deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
442 deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
443 }
444 if ( index == 2 ) {
445 erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
446 erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
447 deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
448 deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
449 }
450 if ( j == columns-2 ) {
451 erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
452 erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
453 deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
454 deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
455 }
456 if ( j == columns-1 ) {
7d855b04 457 erOverEzA(i,j) = 0;
15687d71 458 erOverEzC(i,j) = 0;
7d855b04 459 deltaEzA(i,j) = 0;
15687d71 460 deltaEzC(i,j) = 0;
461 }
462 }
463 }
464
465 }
7d855b04 466
15687d71 467 // AliInfo("Step 1: Interpolation to Standard grid");
468
469 // -------------------------------------------------------------------------------
470 // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
471
7d855b04 472 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
15687d71 473
474 Double_t r, z;//phi, z ;
475 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
476 // phi = fgkPhiList[k] ;
7d855b04 477
15687d71 478 // final lookup table
2bf29b72 479 TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
480 TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
7d855b04 481
15687d71 482 // calculated and integrated tables - just one phi slice
483 TMatrixD &erOverEzA = *arrayofEroverEzA[0] ;
484 TMatrixD &deltaEzA = *arrayofDeltaEzA[0];
485 TMatrixD &erOverEzC = *arrayofEroverEzC[0] ;
7d855b04 486 TMatrixD &deltaEzC = *arrayofDeltaEzC[0];
487
488
15687d71 489 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
490
491 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
7d855b04 492
493 for ( Int_t i = 0 ; i < kNR ; i++ ) {
15687d71 494 r = fgkRList[i] ;
495
496 // Interpolate Lookup tables onto standard grid
497 if (fgkZList[j]>0) {
498 erOverEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzA );
499 deltaEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzA );
500 } else {
501 erOverEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzC );
502 deltaEzFinal(i,j) = - Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzC );
503 // negative coordinate system on C side
504 }
505
506 } // end r loop
507 } // end z loop
508 } // end phi loop
509
7d855b04 510
15687d71 511 // clear the temporary arrays lists
512 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
513
7d855b04 514 delete arrayofErA[k];
15687d71 515 delete arrayofdEzA[k];
7d855b04 516 delete arrayofErC[k];
15687d71 517 delete arrayofdEzC[k];
518
7d855b04 519 delete arrayofEroverEzA[k];
15687d71 520 delete arrayofDeltaEzA[k];
7d855b04 521 delete arrayofEroverEzC[k];
15687d71 522 delete arrayofDeltaEzC[k];
523
524 }
525
526 fZR->Close();
527
528 // ------------------------------------------------------------------------------------------------------
529 // Step 2: Load and sum up lookup table in rphi, fine grid, to emulate for example a GG leak
7d855b04 530
15687d71 531 // AliInfo("Step 2: Preparation of the weighted look-up table");
7d855b04 532
15687d71 533 TFile *fRPhi = new TFile(fSCLookUpPOCsFileNameRPhi.Data(),"READ");
534 if ( !fRPhi ) {
535 AliError("Precalculated POC-looup-table in RPhi could not be found");
536 return;
7d855b04 537 }
15687d71 538
539 // units are in [m]
7d855b04 540 TVector *gridf2 = (TVector*) fRPhi->Get("constants");
15687d71 541 TVector &grid2 = *gridf2;
542 TMatrix *coordf2 = (TMatrix*) fRPhi->Get("coordinates");
543 TMatrix &coord2 = *coordf2;
544 TMatrix *coordPOCf2 = (TMatrix*) fRPhi->Get("POCcoord");
545 TMatrix &coordPOC2 = *coordPOCf2;
546
7d855b04 547 rows = (Int_t)grid2(0); // number of points in r direction
548 phiSlices = (Int_t)grid2(1); // number of points in phi
549 columns = (Int_t)grid2(2); // number of points in z direction
15687d71 550
551 gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
552 Float_t gridSizePhi = TMath::TwoPi()/phiSlices; // unit in [rad]
553 gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
7d855b04 554
15687d71 555 // list of points as used during sum up
7d855b04 556 Double_t rlist2[kNRows], philist2[kNPhiSlices], zedlist2[kNColumns];
15687d71 557 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
558 philist2[k] = gridSizePhi * k;
559 for ( Int_t i = 0 ; i < rows ; i++ ) {
560 rlist2[i] = fgkIFCRadius + i*gridSizeR ;
7d855b04 561 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 562 zedlist2[j] = j * gridSizeZ ;
563 }
564 }
565 } // only done once
7d855b04 566
567 // temporary matrices needed for the calculation
15687d71 568
9f3b99e2 569 TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
7d855b04 570 TMatrixD *arrayofErC2[kNPhiSlices], *arrayofEphiC2[kNPhiSlices], *arrayofdEzC2[kNPhiSlices];
571
572 TMatrixD *arrayofEroverEzA2[kNPhiSlices], *arrayofEphioverEzA2[kNPhiSlices], *arrayofDeltaEzA2[kNPhiSlices];
573 TMatrixD *arrayofEroverEzC2[kNPhiSlices], *arrayofEphioverEzC2[kNPhiSlices], *arrayofDeltaEzC2[kNPhiSlices];
15687d71 574
15687d71 575
15687d71 576 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 577
15687d71 578 arrayofErA2[k] = new TMatrixD(rows,columns) ;
579 arrayofEphiA2[k] = new TMatrixD(rows,columns) ;
580 arrayofdEzA2[k] = new TMatrixD(rows,columns) ;
581 arrayofErC2[k] = new TMatrixD(rows,columns) ;
582 arrayofEphiC2[k] = new TMatrixD(rows,columns) ;
583 arrayofdEzC2[k] = new TMatrixD(rows,columns) ;
584
585 arrayofEroverEzA2[k] = new TMatrixD(rows,columns) ;
7d855b04 586 arrayofEphioverEzA2[k] = new TMatrixD(rows,columns) ;
15687d71 587 arrayofDeltaEzA2[k] = new TMatrixD(rows,columns) ;
588 arrayofEroverEzC2[k] = new TMatrixD(rows,columns) ;
7d855b04 589 arrayofEphioverEzC2[k] = new TMatrixD(rows,columns) ;
15687d71 590 arrayofDeltaEzC2[k] = new TMatrixD(rows,columns) ;
591
7d855b04 592 // zero initialization not necessary, it is done in the constructor of TMatrix
15687d71 593
594 }
7d855b04 595
596
15687d71 597 treePOC = (TTree*)fRPhi->Get("POCall");
598
599 // TVector *bEr = 0; // done above
600 TVector *bEphi= 0;
601 // TVector *bEz = 0; // done above
602
603 treePOC->SetBranchAddress("Er",&bEr);
604 treePOC->SetBranchAddress("Ephi",&bEphi);
605 treePOC->SetBranchAddress("Ez",&bEz);
606
607 // Read the complete tree and do a weighted sum-up over the POC configurations
608 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7d855b04 609
15687d71 610 treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
611 ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
612
613 for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
7d855b04 614
15687d71 615 treePOC->GetEntry(itreepC);
616
617 // center of the POC voxel in [meter]
618 Double_t r0 = coordPOC2(ipC,0);
619 Double_t phi0 = coordPOC2(ipC,1);
620 // Double_t z0 = coordPOC2(ipC,2);
621
622 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
623 // note: coordinates are in [cm]
624 Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, 0.499, 2); // partial load in r,phi
625 Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-0.499, 2); // partial load in r,phi
7d855b04 626
9f3b99e2 627 // printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
15687d71 628
629 // Summing up the vector components according to their weight
630
631 Int_t ip = 0;
7d855b04 632 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 633 for ( Int_t i = 0 ; i < rows ; i++ ) {
634 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 635
15687d71 636 // check wether the coordinates were screwed
7d855b04 637 if (TMath::Abs((coord2(0,ip)*100-rlist2[i]))>1 ||
638 TMath::Abs((coord2(1,ip)-philist2[k]))>1 ||
639 TMath::Abs((coord2(2,ip)*100-zedlist2[j]))>1) {
15687d71 640 AliError("internal error: coordinate system was screwed during the sum-up");
9f3b99e2 641 printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
642 printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist2[i],philist2[k],zedlist2[j]);
15687d71 643 AliError("Don't trust the results of the space charge calculation!");
644 }
7d855b04 645
15687d71 646 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
647 // This will be the most frequent usage (hopefully)
648 // That's why we have to do this here ...
649
650 TMatrixD &erA = *arrayofErA2[k] ;
651 TMatrixD &ephiA = *arrayofEphiA2[k];
652 TMatrixD &dEzA = *arrayofdEzA2[k] ;
7d855b04 653
15687d71 654 TMatrixD &erC = *arrayofErC2[k] ;
655 TMatrixD &ephiC = *arrayofEphiC2[k];
656 TMatrixD &dEzC = *arrayofdEzC2[k] ;
7d855b04 657
15687d71 658 // Sum up - Efield values in [V/m] -> transition to [V/cm]
659 erA(i,j) += ((*bEr)(ip)) * weightA /100;
660 erC(i,j) += ((*bEr)(ip)) * weightC /100;
661 ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
662 ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
663 dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
664 dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
665
666 // increase the counter
667 ip++;
668 }
669 }
670 } // end coordinate loop
7d855b04 671
672
15687d71 673 // Rotation and summation in the rest of the dPhiSteps
674 // which were not stored in the this tree due to storage & symmetry reasons
675
7d855b04 676
15687d71 677 Int_t phiPoints = (Int_t) grid2(1);
678 Int_t phiPOC = (Int_t) grid2(4);
7d855b04 679
15687d71 680 // printf("%d %d\n",phiPOC,flagRadSym);
7d855b04 681
682 for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
683
15687d71 684 Double_t phi0R = phiiC*phi0*2 + phi0; // rotate further
685
686 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
687 // note: coordinates are in [cm] // ecxept z
688 weightA = GetSpaceChargeDensity(r0*100,phi0R, 0.499, 2); // partial load in r,phi
689 weightC = GetSpaceChargeDensity(r0*100,phi0R,-0.499, 2); // partial load in r,phi
7d855b04 690
9f3b99e2 691 // printf("%f %f : %e %e\n",r0,phi0R,weightA,weightC);
7d855b04 692
15687d71 693 // Summing up the vector components according to their weight
694 ip = 0;
7d855b04 695 for ( Int_t j = 0 ; j < columns ; j++ ) {
15687d71 696 for ( Int_t i = 0 ; i < rows ; i++ ) {
697 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 698
15687d71 699 // Note: rotating the coordinated during the sum up
7d855b04 700
15687d71 701 Int_t rotVal = (phiPoints/phiPOC)*phiiC;
702 Int_t ipR = -1;
7d855b04 703
15687d71 704 if ((ip%phiPoints)>=rotVal) {
705 ipR = ip-rotVal;
706 } else {
707 ipR = ip+(phiPoints-rotVal);
708 }
7d855b04 709
15687d71 710 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
7d855b04 711 // This will be the most frequent usage
15687d71 712 // That's why we have to do this here and not outside the loop ...
7d855b04 713
15687d71 714 TMatrixD &erA = *arrayofErA2[k] ;
715 TMatrixD &ephiA = *arrayofEphiA2[k];
716 TMatrixD &dEzA = *arrayofdEzA2[k] ;
7d855b04 717
15687d71 718 TMatrixD &erC = *arrayofErC2[k] ;
719 TMatrixD &ephiC = *arrayofEphiC2[k];
720 TMatrixD &dEzC = *arrayofdEzC2[k] ;
7d855b04 721
15687d71 722 // Sum up - Efield values in [V/m] -> transition to [V/cm]
723 erA(i,j) += ((*bEr)(ipR)) * weightA /100;
724 erC(i,j) += ((*bEr)(ipR)) * weightC /100;
725 ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
726 ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
727 dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
728 dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
729
730 // increase the counter
731 ip++;
732 }
733 }
734 } // end coordinate loop
735
736 } // end phi-POC summation (phiiC)
737
738 ipC++; // POC configuration counter
739
9f3b99e2 740 // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
7d855b04 741
15687d71 742 }
743
744
745
746
747 // -------------------------------------------------------------------------------
748 // Division by the Ez (drift) field and integration along z
749
750 // AliInfo("Step 2: Division and integration");
751
752
753 for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
754
755 // matrices holding the solution - summation of POC charges // see above
756 TMatrixD &erA = *arrayofErA2[k] ;
757 TMatrixD &ephiA = *arrayofEphiA2[k];
758 TMatrixD &dezA = *arrayofdEzA2[k] ;
759 TMatrixD &erC = *arrayofErC2[k] ;
760 TMatrixD &ephiC = *arrayofEphiC2[k];
761 TMatrixD &dezC = *arrayofdEzC2[k] ;
762
763 // matrices which will contain the integrated fields (divided by the drift field)
764 TMatrixD &erOverEzA = *arrayofEroverEzA2[k] ;
765 TMatrixD &ephiOverEzA = *arrayofEphioverEzA2[k];
766 TMatrixD &deltaEzA = *arrayofDeltaEzA2[k];
767 TMatrixD &erOverEzC = *arrayofEroverEzC2[k] ;
768 TMatrixD &ephiOverEzC = *arrayofEphioverEzC2[k];
7d855b04 769 TMatrixD &deltaEzC = *arrayofDeltaEzC2[k];
770
15687d71 771 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
7d855b04 772 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
15687d71 773 // Count backwards to facilitate integration over Z
774
7d855b04 775 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
15687d71 776
7d855b04 777 erOverEzA(i,j) = 0;
15687d71 778 ephiOverEzA(i,j) = 0;
779 deltaEzA(i,j) = 0;
7d855b04 780 erOverEzC(i,j) = 0;
781 ephiOverEzC(i,j) = 0;
15687d71 782 deltaEzC(i,j) = 0;
783
784 for ( Int_t m = j ; m < columns ; m++ ) { // integration
785
786 erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
787 erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
788 ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
789 ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
790 deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
791 deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
792
793 if ( index != 4 ) index = 4; else index = 2 ;
794
795 }
796
797 if ( index == 4 ) {
798 erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
799 erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
800 ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
801 ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
802 deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
803 deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
804 }
805 if ( index == 2 ) {
806 erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
807 erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
808 ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
809 ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
810 deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
811 deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
812 }
813 if ( j == columns-2 ) {
814 erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
815 erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
816 ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
817 ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
818 deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
819 deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
820 }
821 if ( j == columns-1 ) {
7d855b04 822 erOverEzA(i,j) = 0;
15687d71 823 erOverEzC(i,j) = 0;
7d855b04 824 ephiOverEzA(i,j) = 0;
15687d71 825 ephiOverEzC(i,j) = 0;
7d855b04 826 deltaEzA(i,j) = 0;
15687d71 827 deltaEzC(i,j) = 0;
828 }
829 }
830 }
831
832 }
7d855b04 833
15687d71 834 AliInfo("Step 2: Interpolation to Standard grid");
835
836 // -------------------------------------------------------------------------------
837 // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
838
839
840 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
841 Double_t phi = fgkPhiList[k] ;
7d855b04 842
15687d71 843 // final lookup table
2bf29b72 844 TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
845 TMatrixF &ephiOverEzFinal = *fLookUpEphiOverEz[k];
846 TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
7d855b04 847
15687d71 848 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
849
850 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
7d855b04 851
852 for ( Int_t i = 0 ; i < kNR ; i++ ) {
15687d71 853 r = fgkRList[i] ;
854
855 // Interpolate Lookup tables onto standard grid
856 if (fgkZList[j]>0) {
7d855b04 857 erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
15687d71 858 rlist2, zedlist2, philist2, arrayofEroverEzA2 );
859 ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
860 rlist2, zedlist2, philist2, arrayofEphioverEzA2);
861 deltaEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
862 rlist2, zedlist2, philist2, arrayofDeltaEzA2 );
863 } else {
7d855b04 864 erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
15687d71 865 rlist2, zedlist2, philist2, arrayofEroverEzC2 );
866 ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
867 rlist2, zedlist2, philist2, arrayofEphioverEzC2);
868 deltaEzFinal(i,j) -= Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
869 rlist2, zedlist2, philist2, arrayofDeltaEzC2 );
870 }
871
872 } // end r loop
873 } // end z loop
874 } // end phi loop
7d855b04 875
876
15687d71 877 // clear the temporary arrays lists
878 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
879
7d855b04 880 delete arrayofErA2[k];
15687d71 881 delete arrayofEphiA2[k];
882 delete arrayofdEzA2[k];
7d855b04 883 delete arrayofErC2[k];
15687d71 884 delete arrayofEphiC2[k];
885 delete arrayofdEzC2[k];
886
7d855b04 887 delete arrayofEroverEzA2[k];
15687d71 888 delete arrayofEphioverEzA2[k];
889 delete arrayofDeltaEzA2[k];
7d855b04 890 delete arrayofEroverEzC2[k];
15687d71 891 delete arrayofEphioverEzC2[k];
892 delete arrayofDeltaEzC2[k];
893
894 }
7d855b04 895
15687d71 896 fRPhi->Close();
7d855b04 897
15687d71 898 // FINISHED
899
900 fInitLookUp = kTRUE;
901
902}
903
904void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortionCourse() {
7d855b04 905 /// Initialization of the Lookup table which contains the solutions of the
906 /// "space charge" (poisson) problem
907 ///
908 /// The sum-up uses a look-up table which contains different discretized Space charge fields
909 /// in order to calculate the corresponding field deviations due to a given (discretized)
910 /// space charge distribution ....
911 ///
912 /// Method of calculation: Weighted sum-up of the different fields within the look up table
913 /// Note: Full 3d version: Course and slow ...
cc3e558a 914
915 if (fInitLookUp) {
916 AliInfo("Lookup table was already initialized!");
917 // return;
918 }
919
920 AliInfo("Preparation of the weighted look-up table");
7d855b04 921
15687d71 922 TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
923 if ( !f ) {
cc3e558a 924 AliError("Precalculated POC-looup-table could not be found");
925 return;
926 }
927
928 // units are in [m]
7d855b04 929 TVector *gridf = (TVector*) f->Get("constants");
cc3e558a 930 TVector &grid = *gridf;
931 TMatrix *coordf = (TMatrix*) f->Get("coordinates");
932 TMatrix &coord = *coordf;
933 TMatrix *coordPOCf = (TMatrix*) f->Get("POCcoord");
934 TMatrix &coordPOC = *coordPOCf;
7d855b04 935
cc3e558a 936 Bool_t flagRadSym = 0;
937 if (grid(1)==1 && grid(4)==1) {
15687d71 938 // AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
cc3e558a 939 flagRadSym=1;
940 }
941
7d855b04 942 Int_t rows = (Int_t)grid(0); // number of points in r direction
943 Int_t phiSlices = (Int_t)grid(1); // number of points in phi
944 Int_t columns = (Int_t)grid(2); // number of points in z direction
cc3e558a 945
15687d71 946 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
947 const Float_t gridSizePhi = TMath::TwoPi()/phiSlices; // unit in [rad]
948 const Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
7d855b04 949
cc3e558a 950 // temporary matrices needed for the calculation
9f3b99e2 951 TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
952 TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
cc3e558a 953
9f3b99e2 954 TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
955 TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
cc3e558a 956
7d855b04 957
cc3e558a 958 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 959
cc3e558a 960 arrayofErA[k] = new TMatrixD(rows,columns) ;
961 arrayofEphiA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
962 arrayofdEzA[k] = new TMatrixD(rows,columns) ;
963 arrayofErC[k] = new TMatrixD(rows,columns) ;
964 arrayofEphiC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
965 arrayofdEzC[k] = new TMatrixD(rows,columns) ;
966
967 arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
968 arrayofEphioverEzA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
969 arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
970 arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
971 arrayofEphioverEzC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
972 arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
973
974 // Set the values to zero the lookup tables
975 // not necessary, it is done in the constructor of TMatrix - code deleted
976
977 }
7d855b04 978
cc3e558a 979 // list of points as used in the interpolation (during sum up)
9f3b99e2 980 Double_t rlist[kNRows], zedlist[kNColumns] , philist[kNPhiSlices];
cc3e558a 981 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
982 philist[k] = gridSizePhi * k;
983 for ( Int_t i = 0 ; i < rows ; i++ ) {
984 rlist[i] = fgkIFCRadius + i*gridSizeR ;
7d855b04 985 for ( Int_t j = 0 ; j < columns ; j++ ) {
cc3e558a 986 zedlist[j] = j * gridSizeZ ;
987 }
988 }
989 } // only done once
7d855b04 990
991
cc3e558a 992 TTree *treePOC = (TTree*)f->Get("POCall");
993
994 TVector *bEr = 0; TVector *bEphi= 0; TVector *bEz = 0;
7d855b04 995
cc3e558a 996 treePOC->SetBranchAddress("Er",&bEr);
997 if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
998 treePOC->SetBranchAddress("Ez",&bEz);
999
cc3e558a 1000
1001 // Read the complete tree and do a weighted sum-up over the POC configurations
1002 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7d855b04 1003
cc3e558a 1004 Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
1005 Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
1006
15687d71 1007 for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
7d855b04 1008
cc3e558a 1009 treePOC->GetEntry(itreepC);
1010
cc3e558a 1011 // center of the POC voxel in [meter]
1012 Double_t r0 = coordPOC(ipC,0);
1013 Double_t phi0 = coordPOC(ipC,1);
1014 Double_t z0 = coordPOC(ipC,2);
1015
15687d71 1016 ipC++; // POC configuration counter
cc3e558a 1017
1018 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1019 // note: coordinates are in [cm]
92a85338 1020 Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100,0);
1021 Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100,0);
7d855b04 1022
cc3e558a 1023 // Summing up the vector components according to their weight
1024
1025 Int_t ip = 0;
7d855b04 1026 for ( Int_t j = 0 ; j < columns ; j++ ) {
cc3e558a 1027 for ( Int_t i = 0 ; i < rows ; i++ ) {
1028 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 1029
cc3e558a 1030 // check wether the coordinates were screwed
7d855b04 1031 if (TMath::Abs((coord(0,ip)*100-rlist[i]))>1 ||
1032 TMath::Abs((coord(1,ip)-philist[k]))>1 ||
1033 TMath::Abs((coord(2,ip)*100-zedlist[j]))>1) {
cc3e558a 1034 AliError("internal error: coordinate system was screwed during the sum-up");
9f3b99e2 1035 printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
1036 printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist[i],philist[k],zedlist[j]);
15687d71 1037 AliError("Don't trust the results of the space charge calculation!");
cc3e558a 1038 }
7d855b04 1039
cc3e558a 1040 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1041 // This will be the most frequent usage (hopefully)
1042 // That's why we have to do this here ...
1043
1044 TMatrixD &erA = *arrayofErA[k] ;
1045 TMatrixD &ephiA = *arrayofEphiA[k];
15687d71 1046 TMatrixD &dEzA = *arrayofdEzA[k] ;
7d855b04 1047
cc3e558a 1048 TMatrixD &erC = *arrayofErC[k] ;
1049 TMatrixD &ephiC = *arrayofEphiC[k];
15687d71 1050 TMatrixD &dEzC = *arrayofdEzC[k] ;
7d855b04 1051
cc3e558a 1052 // Sum up - Efield values in [V/m] -> transition to [V/cm]
1053 erA(i,j) += ((*bEr)(ip)) * weightA /100;
1054 erC(i,j) += ((*bEr)(ip)) * weightC /100;
1055 if (!flagRadSym) {
15687d71 1056 ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
1057 ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
cc3e558a 1058 }
1059 dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
1060 dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
1061
1062 // increase the counter
1063 ip++;
1064 }
1065 }
15687d71 1066 } // end coordinate loop
7d855b04 1067
1068
cc3e558a 1069 // Rotation and summation in the rest of the dPhiSteps
15687d71 1070 // which were not stored in the this tree due to storage & symmetry reasons
cc3e558a 1071
1072 Int_t phiPoints = (Int_t) grid(1);
1073 Int_t phiPOC = (Int_t) grid(4);
7d855b04 1074
15687d71 1075 // printf("%d %d\n",phiPOC,flagRadSym);
7d855b04 1076
1077 for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
1078
cc3e558a 1079 r0 = coordPOC(ipC,0);
1080 phi0 = coordPOC(ipC,1);
1081 z0 = coordPOC(ipC,2);
7d855b04 1082
cc3e558a 1083 ipC++; // POC conf. counter
7d855b04 1084
cc3e558a 1085 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1086 // note: coordinates are in [cm]
7d855b04 1087 weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100,0);
92a85338 1088 weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100,0);
7d855b04 1089
9f3b99e2 1090 // printf("%f %f %f: %e %e\n",r0,phi0,z0,weightA,weightC);
7d855b04 1091
cc3e558a 1092 // Summing up the vector components according to their weight
1093 ip = 0;
7d855b04 1094 for ( Int_t j = 0 ; j < columns ; j++ ) {
cc3e558a 1095 for ( Int_t i = 0 ; i < rows ; i++ ) {
1096 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
7d855b04 1097
cc3e558a 1098 // Note: rotating the coordinated during the sum up
7d855b04 1099
cc3e558a 1100 Int_t rotVal = (phiPoints/phiPOC)*phiiC;
1101 Int_t ipR = -1;
7d855b04 1102
cc3e558a 1103 if ((ip%phiPoints)>=rotVal) {
1104 ipR = ip-rotVal;
1105 } else {
1106 ipR = ip+(phiPoints-rotVal);
1107 }
7d855b04 1108
cc3e558a 1109 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
7d855b04 1110 // This will be the most frequent usage
cc3e558a 1111 // That's why we have to do this here and not outside the loop ...
7d855b04 1112
cc3e558a 1113 TMatrixD &erA = *arrayofErA[k] ;
1114 TMatrixD &ephiA = *arrayofEphiA[k];
1115 TMatrixD &dEzA = *arrayofdEzA[k] ;
7d855b04 1116
cc3e558a 1117 TMatrixD &erC = *arrayofErC[k] ;
1118 TMatrixD &ephiC = *arrayofEphiC[k];
1119 TMatrixD &dEzC = *arrayofdEzC[k] ;
7d855b04 1120
cc3e558a 1121 // Sum up - Efield values in [V/m] -> transition to [V/cm]
1122 erA(i,j) += ((*bEr)(ipR)) * weightA /100;
1123 erC(i,j) += ((*bEr)(ipR)) * weightC /100;
1124 if (!flagRadSym) {
15687d71 1125 ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
1126 ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
cc3e558a 1127 }
1128 dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
1129 dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
1130
1131 // increase the counter
1132 ip++;
1133 }
1134 }
1135 } // end coordinate loop
1136
1137 } // end phi-POC summation (phiiC)
7d855b04 1138
cc3e558a 1139
9f3b99e2 1140 // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
7d855b04 1141
15687d71 1142 }
1143
cc3e558a 1144
cc3e558a 1145
1146 // -------------------------------------------------------------------------------
1147 // Division by the Ez (drift) field and integration along z
1148
15687d71 1149 AliInfo("Division and integration");
1150
cc3e558a 1151 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
1152
1153 for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
1154
15687d71 1155 // matrices holding the solution - summation of POC charges // see above
cc3e558a 1156 TMatrixD &erA = *arrayofErA[k] ;
1157 TMatrixD &ephiA = *arrayofEphiA[k];
1158 TMatrixD &dezA = *arrayofdEzA[k] ;
1159 TMatrixD &erC = *arrayofErC[k] ;
1160 TMatrixD &ephiC = *arrayofEphiC[k];
1161 TMatrixD &dezC = *arrayofdEzC[k] ;
1162
15687d71 1163 // matrices which will contain the integrated fields (divided by the drift field)
cc3e558a 1164 TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
1165 TMatrixD &ephiOverEzA = *arrayofEphioverEzA[k];
1166 TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
1167 TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
1168 TMatrixD &ephiOverEzC = *arrayofEphioverEzC[k];
7d855b04 1169 TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
1170
cc3e558a 1171 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
7d855b04 1172 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
cc3e558a 1173 // Count backwards to facilitate integration over Z
1174
7d855b04 1175 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
cc3e558a 1176
1177 erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0; deltaEzA(i,j) = 0;
1178 erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0; deltaEzC(i,j) = 0;
1179
1180 for ( Int_t m = j ; m < columns ; m++ ) { // integration
1181
1182 erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
1183 erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
1184 if (!flagRadSym) {
1185 ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
1186 ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
1187 }
15687d71 1188 deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
1189 deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
cc3e558a 1190
1191 if ( index != 4 ) index = 4; else index = 2 ;
1192
1193 }
1194
1195 if ( index == 4 ) {
1196 erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
1197 erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
1198 if (!flagRadSym) {
1199 ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
1200 ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
1201 }
15687d71 1202 deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
1203 deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
cc3e558a 1204 }
1205 if ( index == 2 ) {
1206 erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
1207 erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
1208 if (!flagRadSym) {
1209 ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
1210 ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
1211 }
15687d71 1212 deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
1213 deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
cc3e558a 1214 }
1215 if ( j == columns-2 ) {
1216 erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
1217 erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
1218 if (!flagRadSym) {
1219 ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
1220 ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
1221 }
15687d71 1222 deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
1223 deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
cc3e558a 1224 }
1225 if ( j == columns-1 ) {
7d855b04 1226 erOverEzA(i,j) = 0;
15687d71 1227 erOverEzC(i,j) = 0;
cc3e558a 1228 if (!flagRadSym) {
7d855b04 1229 ephiOverEzA(i,j) = 0;
15687d71 1230 ephiOverEzC(i,j) = 0;
cc3e558a 1231 }
7d855b04 1232 deltaEzA(i,j) = 0;
15687d71 1233 deltaEzC(i,j) = 0;
cc3e558a 1234 }
1235 }
1236 }
1237
1238 }
15687d71 1239
7d855b04 1240
1241
cc3e558a 1242 AliInfo("Interpolation to Standard grid");
1243
1244 // -------------------------------------------------------------------------------
15687d71 1245 // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
cc3e558a 1246
7d855b04 1247 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
cc3e558a 1248
1249 Double_t r, phi, z ;
1250 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1251 phi = fgkPhiList[k] ;
7d855b04 1252
2bf29b72 1253 TMatrixF &erOverEz = *fLookUpErOverEz[k] ;
1254 TMatrixF &ephiOverEz = *fLookUpEphiOverEz[k];
1255 TMatrixF &deltaEz = *fLookUpDeltaEz[k] ;
7d855b04 1256
cc3e558a 1257 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1258
1259 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
7d855b04 1260
1261 for ( Int_t i = 0 ; i < kNR ; i++ ) {
cc3e558a 1262 r = fgkRList[i] ;
1263
1264 // Interpolate Lookup tables onto standard grid
1265 if (fgkZList[j]>0) {
7d855b04 1266 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
cc3e558a 1267 rlist, zedlist, philist, arrayofEroverEzA );
1268 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1269 rlist, zedlist, philist, arrayofEphioverEzA);
1270 deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1271 rlist, zedlist, philist, arrayofDeltaEzA );
1272 } else {
7d855b04 1273 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
cc3e558a 1274 rlist, zedlist, philist, arrayofEroverEzC );
1275 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1276 rlist, zedlist, philist, arrayofEphioverEzC);
1277 deltaEz(i,j) = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1278 rlist, zedlist, philist, arrayofDeltaEzC );
1279 // negative coordinate system on C side
1280 }
1281
1282 } // end r loop
1283 } // end z loop
1284 } // end phi loop
1285
7d855b04 1286
cc3e558a 1287 // clear the temporary arrays lists
1288 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1289
7d855b04 1290 delete arrayofErA[k];
cc3e558a 1291 delete arrayofEphiA[k];
1292 delete arrayofdEzA[k];
7d855b04 1293 delete arrayofErC[k];
cc3e558a 1294 delete arrayofEphiC[k];
1295 delete arrayofdEzC[k];
1296
7d855b04 1297 delete arrayofEroverEzA[k];
cc3e558a 1298 delete arrayofEphioverEzA[k];
1299 delete arrayofDeltaEzA[k];
7d855b04 1300 delete arrayofEroverEzC[k];
cc3e558a 1301 delete arrayofEphioverEzC[k];
1302 delete arrayofDeltaEzC[k];
1303
1304 }
1305
cc3e558a 1306 fInitLookUp = kTRUE;
1307
1308}
1309
1310
15687d71 1311void AliTPCSpaceCharge3D::SetSCDataFileName(TString fname) {
7d855b04 1312 /// Set & load the Space charge density distribution from a file
1313 /// (linear interpolation onto a standard grid)
1314
cc3e558a 1315
1316 fSCDataFileName = fname;
1317
15687d71 1318 TFile *f = new TFile(fSCDataFileName.Data(),"READ");
7d855b04 1319 if (!f) {
cc3e558a 1320 AliError(Form("File %s, which should contain the space charge distribution, could not be found",
15687d71 1321 fSCDataFileName.Data()));
cc3e558a 1322 return;
1323 }
7d855b04 1324
15687d71 1325 TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
7d855b04 1326 if (!densityRZ) {
cc3e558a 1327 AliError(Form("The indicated file (%s) does not contain a histogram called %s",
15687d71 1328 fSCDataFileName.Data(),"SpaceChargeInRZ"));
1329 return;
1330 }
1331
1332 TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
7d855b04 1333 if (!densityRPhi) {
15687d71 1334 AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1335 fSCDataFileName.Data(),"SpaceChargeInRPhi"));
cc3e558a 1336 return;
1337 }
7d855b04 1338
cc3e558a 1339
1340 Double_t r, phi, z ;
15687d71 1341
1342 TMatrixD &scDensityInRZ = *fSCdensityInRZ;
1343 TMatrixD &scDensityInRPhiA = *fSCdensityInRPhiA;
1344 TMatrixD &scDensityInRPhiC = *fSCdensityInRPhiC;
cc3e558a 1345 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1346 phi = fgkPhiList[k] ;
2bf29b72 1347 TMatrixF &scDensity = *fSCdensityDistribution[k] ;
cc3e558a 1348 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
7d855b04 1349 z = fgkZList[j] ;
1350 for ( Int_t i = 0 ; i < kNR ; i++ ) {
cc3e558a 1351 r = fgkRList[i] ;
15687d71 1352
1353 // partial load in (r,z)
1354 if (k==0) // do just once
1355 scDensityInRZ(i,j) = densityRZ->Interpolate(r,z);
1356
1357 // partial load in (r,phi)
1358 if ( j==0 || j == kNZ/2 ) {
7d855b04 1359 if (z>0)
15687d71 1360 scDensityInRPhiA(i,k) = densityRPhi->Interpolate(r,phi,0.499); // A side
7d855b04 1361 else
15687d71 1362 scDensityInRPhiC(i,k) = densityRPhi->Interpolate(r,phi,-0.499); // C side
1363 }
1364
1365 // Full 3D configuration ...
7d855b04 1366 if (z>0)
1367 scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
15687d71 1368 else
7d855b04 1369 scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
cc3e558a 1370 }
1371 }
1372 }
1373
cc3e558a 1374 f->Close();
1375
1376 fInitLookUp = kFALSE;
1377
7d855b04 1378
cc3e558a 1379}
1380
92a85338 1381void AliTPCSpaceCharge3D::SetInputSpaceCharge(TH3 * hisSpaceCharge3D, TH2 * hisRPhi, TH2* hisRZ, Double_t norm){
7d855b04 1382 /// Use 3D space charge map as an optional input
1383 /// The layout of the input histogram is assumed to be: (phi,r,z)
1384 /// Density histogram is expreseed is expected to bin in C/m^3
1385 ///
1386 /// Standard histogram interpolation is used in order to use the density at center of voxel
1387
92a85338 1388 fSpaceChargeHistogram3D = hisSpaceCharge3D;
1389 fSpaceChargeHistogramRPhi = hisRPhi;
1390 fSpaceChargeHistogramRZ = hisRZ;
1391
1392 Double_t r, phi, z ;
1393 TMatrixD &scDensityInRZ = *fSCdensityInRZ;
1394 TMatrixD &scDensityInRPhiA = *fSCdensityInRPhiA;
1395 TMatrixD &scDensityInRPhiC = *fSCdensityInRPhiC;
1396 //
1397 Double_t rmin=hisSpaceCharge3D->GetYaxis()->GetBinCenter(0);
1398 Double_t rmax=hisSpaceCharge3D->GetYaxis()->GetBinUpEdge(hisSpaceCharge3D->GetYaxis()->GetNbins());
1399 Double_t zmin=hisSpaceCharge3D->GetZaxis()->GetBinCenter(0);
1400 Double_t zmax=hisSpaceCharge3D->GetZaxis()->GetBinCenter(hisSpaceCharge3D->GetZaxis()->GetNbins());
1401
1402 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1403 phi = fgkPhiList[k] ;
1404 TMatrixF &scDensity = *fSCdensityDistribution[k] ;
1405 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
7d855b04 1406 z = fgkZList[j] ;
1407 for ( Int_t i = 0 ; i < kNR ; i++ ) {
92a85338 1408 // Full 3D configuration ...
1409 r = fgkRList[i] ;
7d855b04 1410 if (r>rmin && r<rmax && z>zmin && z< zmax){
92a85338 1411 // partial load in (r,z)
1412 if (k==0) {
1413 if (fSpaceChargeHistogramRZ) scDensityInRZ(i,j) = norm*fSpaceChargeHistogramRZ->Interpolate(r,z) ;
1414 }
1415 // partial load in (r,phi)
1416 if ( (j==0 || j == kNZ/2) && fSpaceChargeHistogramRPhi) {
7d855b04 1417 if (z>0)
92a85338 1418 scDensityInRPhiA(i,k) = norm*fSpaceChargeHistogramRPhi->Interpolate(phi,r); // A side
7d855b04 1419 else
92a85338 1420 scDensityInRPhiC(i,k) = norm*fSpaceChargeHistogramRPhi->Interpolate(phi+TMath::TwoPi(),r); // C side
1421 }
7d855b04 1422
1423 if (z>0)
1424 scDensity(i,j) = norm*fSpaceChargeHistogram3D->Interpolate(phi,r,z);
92a85338 1425 else
1426 scDensity(i,j) = norm*fSpaceChargeHistogram3D->Interpolate(phi,r,z);
1427 }
1428 }
1429 }
1430 }
7d855b04 1431
92a85338 1432 fInitLookUp = kFALSE;
1433
1434}
1435
cc3e558a 1436
15687d71 1437Float_t AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
7d855b04 1438 /// returns the (input) space charge density at a given point according
1439 /// Note: input in [cm], output in [C/m^3/e0] !!
cc3e558a 1440
cc3e558a 1441 while (phi<0) phi += TMath::TwoPi();
1442 while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
1443
1444
1445 // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
1446 Int_t order = 1; //
15687d71 1447 Float_t sc = 0;
cc3e558a 1448
15687d71 1449 if (mode == 0) { // return full load
7d855b04 1450 sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
15687d71 1451 fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
7d855b04 1452
15687d71 1453 } else if (mode == 1) { // return partial load in (r,z)
1454 TMatrixD &scDensityInRZ = *fSCdensityInRZ;
1455 sc = Interpolate2DTable(order, r, z, kNR, kNZ, fgkRList, fgkZList, scDensityInRZ );
7d855b04 1456
15687d71 1457 } else if (mode == 2) { // return partial load in (r,phi)
1458
1459 if (z>0) {
1460 TMatrixD &scDensityInRPhi = *fSCdensityInRPhiA;
1461 sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1462 } else {
1463 TMatrixD &scDensityInRPhi = *fSCdensityInRPhiC;
1464 sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1465 }
1466
1467 } else {
1468 // should i give a warning?
1469 sc = 0;
1470 }
7d855b04 1471
9f3b99e2 1472 // printf("%f %f %f: %f\n",r,phi,z,sc);
7d855b04 1473
cc3e558a 1474 return sc;
1475}
1476
1477
15687d71 1478TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
7d855b04 1479 /// return a simple histogramm containing the space charge distribution (input for the calculation)
cc3e558a 1480
1481 TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1482 nx,-250.,250.,ny,-250.,250.);
1483
1484 for (Int_t iy=1;iy<=ny;++iy) {
1485 Double_t yp = h->GetYaxis()->GetBinCenter(iy);
1486 for (Int_t ix=1;ix<=nx;++ix) {
1487 Double_t xp = h->GetXaxis()->GetBinCenter(ix);
7d855b04 1488
cc3e558a 1489 Float_t r = TMath::Sqrt(xp*xp+yp*yp);
7d855b04 1490 Float_t phi = TMath::ATan2(yp,xp);
1491
cc3e558a 1492 if (85.<=r && r<=250.) {
15687d71 1493 Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
7d855b04 1494 h->SetBinContent(ix,iy,sc);
cc3e558a 1495 } else {
1496 h->SetBinContent(ix,iy,0.);
1497 }
1498 }
1499 }
7d855b04 1500
cc3e558a 1501 return h;
7d855b04 1502}
cc3e558a 1503
15687d71 1504TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
7d855b04 1505 /// return a simple histogramm containing the space charge distribution (input for the calculation)
cc3e558a 1506
1507 TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1508 nz,-250.,250.,nr,85.,250.);
1509
1510 for (Int_t ir=1;ir<=nr;++ir) {
1511 Float_t r = h->GetYaxis()->GetBinCenter(ir);
1512 for (Int_t iz=1;iz<=nz;++iz) {
1513 Float_t z = h->GetXaxis()->GetBinCenter(iz);
15687d71 1514 Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
cc3e558a 1515 h->SetBinContent(iz,ir,sc);
1516 }
1517 }
1518
1519 return h;
7d855b04 1520}
cc3e558a 1521
1522void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
7d855b04 1523 /// Example on how to write a Space charge distribution into a File
1524 /// (see below: estimate from scaling STAR measurements to Alice)
1525 /// Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1526 /// the needed calculation time
cc3e558a 1527
1528 TFile *f = new TFile(fname,"RECREATE");
7d855b04 1529
cc3e558a 1530 // some grid, not too course
15687d71 1531 Int_t nr = 350;
cc3e558a 1532 Int_t nphi = 180;
15687d71 1533 Int_t nz = 500;
cc3e558a 1534
1535 Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1536 Double_t dphi = TMath::TwoPi()/(nphi+1);
1537 Double_t dz = 500./(nz+1);
1538 Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
1539
15687d71 1540
1541 // Charge distribution in ZR (rotational symmetric) ------------------
1542
1543 TH2F *histoZR = new TH2F("chargeZR","chargeZR",
1544 nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1545 nz,-250-dz-safty,250+dz+safty);
7d855b04 1546
cc3e558a 1547 for (Int_t ir=1;ir<=nr;++ir) {
15687d71 1548 Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1549 for (Int_t iz=1;iz<=nz;++iz) {
1550 Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
7d855b04 1551
15687d71 1552 // recalculation to meter
1553 Double_t lZ = 2.5; // approx. TPC drift length
1554 Double_t rpM = rp/100.; // in [m]
1555 Double_t zpM = TMath::Abs(zp/100.); // in [m]
7d855b04 1556
15687d71 1557 // setting of mb multiplicity and Interaction rate
1558 Double_t multiplicity = 950;
1559 Double_t intRate = 7800;
cc3e558a 1560
15687d71 1561 // calculation of "scaled" parameters
1562 Double_t a = multiplicity*intRate/79175;
1563 Double_t b = a/lZ;
7d855b04 1564
15687d71 1565 Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
7d855b04 1566
15687d71 1567 charge = charge*fgke0; // [C/m^3]
7d855b04 1568
15687d71 1569 if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
cc3e558a 1570
15687d71 1571 // charge = 0; // for tests
7d855b04 1572 histoZR->SetBinContent(ir,iz,charge);
15687d71 1573 }
1574 }
7d855b04 1575
15687d71 1576 histoZR->Write("SpaceChargeInRZ");
7d855b04 1577
15687d71 1578
1579 // Charge distribution in RPhi (e.g. Floating GG wire) ------------
7d855b04 1580
15687d71 1581 TH3F *histoRPhi = new TH3F("chargeRPhi","chargeRPhi",
1582 nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1583 nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
1584 2,-1,1); // z part - to allow A and C side differences
7d855b04 1585
15687d71 1586 // some 'arbitrary' GG leaks
1587 Int_t nGGleaks = 5;
1588 Double_t secPosA[5] = {3,6,6,11,13}; // sector
752b0cc7 1589 Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm
1590 Double_t secPosC[5] = {1,8,12,15,15}; // sector
1591 Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm
15687d71 1592
1593 for (Int_t ir=1;ir<=nr;++ir) {
1594 Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir);
1595 for (Int_t iphi=1;iphi<=nphi;++iphi) {
1596 Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi);
1597 for (Int_t iz=1;iz<=2;++iz) {
1598 Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz);
7d855b04 1599
15687d71 1600 Double_t charge = 0;
7d855b04 1601
15687d71 1602 for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks
7d855b04 1603
15687d71 1604 // A side
7d855b04 1605 Double_t secPos = secPosA[igg];
15687d71 1606 Double_t radialPos = radialPosA[igg];
1607
1608 if (zp<0) { // C side
7d855b04 1609 secPos = secPosC[igg];
15687d71 1610 radialPos = radialPosC[igg];
7d855b04 1611 }
15687d71 1612
1613 // some 'arbitrary' GG leaks
1614 if ( (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice
1615 if ( rp>(radialPos-2.5) && rp<(radialPos+2.5)) // 5 cm slice
7d855b04 1616 charge = 300;
15687d71 1617 }
7d855b04 1618
1619 }
1620
cc3e558a 1621 charge = charge*fgke0; // [C/m^3]
1622
7d855b04 1623 histoRPhi->SetBinContent(ir,iphi,iz,charge);
cc3e558a 1624 }
1625 }
1626 }
1627
15687d71 1628 histoRPhi->Write("SpaceChargeInRPhi");
cc3e558a 1629
cc3e558a 1630 f->Close();
7d855b04 1631
cc3e558a 1632}
1633
1634
1635void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
7d855b04 1636 /// Print function to check the settings of the boundary vectors
1637 /// option=="a" prints the C0 and C1 coefficents for calibration purposes
cc3e558a 1638
1639 TString opt = option; opt.ToLower();
1640 printf("%s\n",GetTitle());
1641 printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
1642 printf(" SC correction factor: %f \n",fCorrectionFactor);
1643
1644 if (opt.Contains("a")) { // Print all details
1645 printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
1646 printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
7d855b04 1647 }
1648
cc3e558a 1649 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
1650
7d855b04 1651}
92a85338 1652
1653
1654
1655void AliTPCSpaceCharge3D::InitSpaceCharge3DPoisson(Int_t kRows, Int_t kColumns, Int_t kPhiSlices, Int_t kIterations){
7d855b04 1656 /// MI extension - calculate E field
1657 /// - inspired by AliTPCROCVoltError3D::InitROCVoltError3D()
1658 /// Initialization of the Lookup table which contains the solutions of the
1659 /// Dirichlet boundary problem
1660 /// Calculation of the single 3D-Poisson solver is done just if needed
1661 /// (see basic lookup tables in header file)
1662
1663 Int_t kPhiSlicesPerSector = kPhiSlices/18;
92a85338 1664 //
7d855b04 1665 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
92a85338 1666 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
1667 const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
1668 const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
1669
1670 // temporary arrays to create the boundary conditions
7d855b04 1671 TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
1672 TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
92a85338 1673
1674 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1675 arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
1676 arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
1677 arrayofEroverEz[k] = new TMatrixD(kRows,kColumns) ;
1678 arrayofEphioverEz[k] = new TMatrixD(kRows,kColumns) ;
1679 arrayofDeltaEz[k] = new TMatrixD(kRows,kColumns) ;
1680 }
7d855b04 1681
92a85338 1682 // list of point as used in the poisson relation and the interpolation (during sum up)
1683 Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
1684 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1685 philist[k] = gridSizePhi * k;
1686 for ( Int_t i = 0 ; i < kRows ; i++ ) {
1687 rlist[i] = fgkIFCRadius + i*gridSizeR ;
1688 for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
1689 zedlist[j] = j * gridSizeZ ;
1690 }
1691 }
1692 }
1693
1694 // ==========================================================================
1695 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
1696 // Allow for different size grid spacing in R and Z directions
7d855b04 1697
92a85338 1698 const Int_t symmetry = 0;
7d855b04 1699
92a85338 1700 // Set bondaries and solve Poisson's equation --------------------------
7d855b04 1701
92a85338 1702 if ( !fInitLookUp ) {
7d855b04 1703
92a85338 1704 AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
7d855b04 1705
92a85338 1706 for ( Int_t side = 0 ; side < 2 ; side++ ) { // Solve Poisson3D twice; once for +Z and once for -Z
1707 AliSysInfo::AddStamp("RunSide", 1,side,0);
1708 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1709 TMatrixD &arrayV = *arrayofArrayV[k] ;
1710 TMatrixD &charge = *arrayofCharge[k] ;
7d855b04 1711
92a85338 1712 //Fill arrays with initial conditions. V on the boundary and Charge in the volume.
1713// for ( Int_t i = 0 ; i < kRows ; i++ ) {
1714// for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
7d855b04 1715// arrayV(i,j) = 0.0 ;
92a85338 1716// charge(i,j) = 0.0 ;
1717
1718// // Float_t radius0 = rlist[i] ;
1719// // Float_t phi0 = gridSizePhi * k ;
7d855b04 1720
92a85338 1721// // To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
1722// // if ( j == (kColumns-1) ) {
1723// // arrayV(i,j) = 0.5* ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
1724
1725// // if (side==1) // C side
1726// // arrayV(i,j) = -arrayV(i,j); // minus sign on the C side to allow a consistent usage of global z when setting the boundaries
1727// // }
1728// }
7d855b04 1729// }
1730
1731 for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
1732 for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
92a85338 1733 Float_t radius0 = rlist[i] ;
1734 Float_t phi0 = gridSizePhi * k ;
1735 Double_t z0 = zedlist[j];
1736 if (side==1) z0= -TMath::Abs(zedlist[j]);
7d855b04 1737 arrayV(i,j) = 0.0 ;
92a85338 1738 charge(i,j) = fSpaceChargeHistogram3D->Interpolate(phi0,radius0,z0);
1739 }
1740 }
7d855b04 1741 }
92a85338 1742 AliSysInfo::AddStamp("RunPoisson", 2,side,0);
7d855b04 1743
92a85338 1744 // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
1745 // Allow for different size grid spacing in R and Z directions
7d855b04 1746
1747 // PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
92a85338 1748 // arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
7d855b04 1749 // kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
92a85338 1750 // symmetry , fROCdisplacement) ;
7d855b04 1751 PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
92a85338 1752 arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
7d855b04 1753 kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
92a85338 1754 symmetry ) ;
7d855b04 1755
92a85338 1756 //Interpolate results onto a custom grid which is used just for these calculations.
1757 Double_t r, phi, z ;
1758 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1759 phi = fgkPhiList[k] ;
7d855b04 1760
92a85338 1761 TMatrixF &erOverEz = *fLookUpErOverEz[k] ;
1762 TMatrixF &ephiOverEz = *fLookUpEphiOverEz[k];
1763 TMatrixF &deltaEz = *fLookUpDeltaEz[k] ;
7d855b04 1764
92a85338 1765 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1766
1767 z = TMath::Abs(fgkZList[j]) ; // Symmetric solution in Z that depends only on ABS(Z)
7d855b04 1768
92a85338 1769 if ( side == 0 && fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
1770 if ( side == 1 && fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
7d855b04 1771
1772 for ( Int_t i = 0 ; i < kNR ; i++ ) {
92a85338 1773 r = fgkRList[i] ;
1774
1775 // Interpolate basicLookup tables; once for each rod, then sum the results
7d855b04 1776 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
92a85338 1777 rlist, zedlist, philist, arrayofEroverEz );
1778 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
1779 rlist, zedlist, philist, arrayofEphioverEz);
1780 deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
1781 rlist, zedlist, philist, arrayofDeltaEz );
1782
1783 if (side == 1) deltaEz(i,j) = - deltaEz(i,j); // negative coordinate system on C side
1784
1785 } // end r loop
1786 }// end z loop
1787 }// end phi loop
7d855b04 1788 AliSysInfo::AddStamp("Interpolate Poisson", 3,side,0);
92a85338 1789 if ( side == 0 ) AliInfo(" A side done");
1790 if ( side == 1 ) AliInfo(" C side done");
1791 } // end side loop
1792 }
7d855b04 1793
92a85338 1794 // clear the temporary arrays lists
1795 for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
1796 delete arrayofArrayV[k];
1797 delete arrayofCharge[k];
7d855b04 1798 delete arrayofEroverEz[k];
92a85338 1799 delete arrayofEphioverEz[k];
1800 delete arrayofDeltaEz[k];
1801 }
7d855b04 1802
92a85338 1803
1804 fInitLookUp = kTRUE;
1805
cc3e558a 1806}
92a85338 1807
8a94851d 1808
1809
1810AliTPCSpaceCharge3D * AliTPCSpaceCharge3D::MakeCorrection22D(const char* fileName, Double_t multiplicity, Double_t intRate, Double_t epsIROC, Double_t epsOROC,
7d855b04 1811 Double_t gasfactor,
8a94851d 1812 Double_t radialScaling){
7d855b04 1813 /// Origin: Christian Lippmann, CERN, Christian.Lippmann@cern.ch based on the internal note (xxx ...)
1814 /// adopted by Marian Ivanov (different epsilon in IROC and OROC)
1815 ///
1816 /// Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1817 /// the needed calculation time.
1818 ///
1819 /// Explanation of variables:
1820 /// 1) multiplicity: charghed particle dn/deta for top 80% centrality (660 for 2011,
1821 /// expect 950 for full energy)
1822 /// 2) intRate: Total interaction rate (e.g. 50kHz for the upgrade)
1823 /// 3) eps: Number of backdrifting ions per primary electron (0 for MWPC, e.g.10 for GEM)
1824 /// 4) gasfactor: Use different gas. E.g. Ar/CO2 has twice the primary ionization, ion drift
1825 /// velocity factor 2.5 slower, so gasfactor = 5.
1826
1827
1828 TFile *f = new TFile(fileName, "RECREATE");
8a94851d 1829 // some grid, not too coarse
1830 const Int_t nr = 350;
1831 const Int_t nphi = 180;
1832 const Int_t nz = 500;
7d855b04 1833 const Double_t kROROC=134; // hardwired OROC radius
8a94851d 1834
1835 Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1836 Double_t dphi = TMath::TwoPi()/(nphi+1);
1837 Double_t dz = 500./(nz+1);
1838 Double_t safty = 0.; // due to a root bug which does not interpolate the boundary ..
1839 // .. (first and last bin) correctly
1840
1841 // Charge distribution in ZR (rotational symmetric) ------------------
1842
1843 TH2F *histoZR = new TH2F("chargeZR", "chargeZR",
1844 nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty,
1845 nz, -250-dz-safty, 250+dz+safty);
1846
1847 // For the normalization to same integral as radial exponent = 2
1848 Double_t radialExponent = -2.; // reference = 2
1849 Double_t radiusInner = histoZR->GetXaxis()->GetBinCenter(1) / 100.;//in [m]
1850 Double_t radiusOuter = histoZR->GetXaxis()->GetBinCenter(nr) / 100.;//in [m]
7d855b04 1851 Double_t integralRadialExponent2 = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1)
8a94851d 1852 - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1);
7d855b04 1853
1854 radialExponent = -radialScaling; // user set
8a94851d 1855 Double_t integralRadialExponentUser = 0.;
1856 if(radialScaling > 1 + 0.000001 || radialScaling < 1 - 0.000001 ) // to avoid n = -1
7d855b04 1857 integralRadialExponentUser = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1)
8a94851d 1858 - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1);
1859 else
1860 integralRadialExponentUser = TMath::Log(radiusOuter) - TMath::Log(radiusInner);
7d855b04 1861
8a94851d 1862 Double_t normRadialExponent = integralRadialExponent2 / integralRadialExponentUser;
7d855b04 1863
8a94851d 1864 for (Int_t ir=1;ir<=nr;++ir) {
1865 Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1866 for (Int_t iz=1;iz<=nz;++iz) {
1867 Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
7d855b04 1868 Double_t eps = (rp <kROROC) ? epsIROC:epsOROC;
8a94851d 1869 // recalculation to meter
1870 Double_t lZ = 2.5; // approx. TPC drift length
1871 Double_t rpM = rp/100.; // in [m]
1872 Double_t zpM = TMath::Abs(zp/100.); // in [m]
7d855b04 1873
8a94851d 1874 // calculation of "scaled" parameters
1875 Double_t a = multiplicity*intRate/76628;
1876 //Double_t charge = gasfactor * ( a / (rpM*rpM) * (1 - zpM/lZ) ); // charge in [C/m^3/e0], no IBF
1877 Double_t charge = normRadialExponent * gasfactor * ( a / (TMath::Power(rpM,radialScaling)) * (1 - zpM/lZ + eps) ); // charge in [C/m^3/e0], with IBF
7d855b04 1878
8a94851d 1879 charge = charge*fgke0; // [C/m^3]
1880 if (zp<0) charge *= 0.9; // Slightly less on C side due to front absorber
1881
7d855b04 1882 histoZR->SetBinContent(ir, iz, charge);
8a94851d 1883 }
1884 }
7d855b04 1885
8a94851d 1886 histoZR->Write("SpaceChargeInRZ");
1887 //
1888 // Charge distribution in RPhi (e.g. Floating GG wire) ------------
1889 //
1890 TH3F *histoRPhi = new TH3F("chargeRPhi", "chargeRPhi",
1891 nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty,
1892 nphi, 0-dphi-safty, TMath::TwoPi()+dphi+safty,
7d855b04 1893 2, -1, 1); // z part - to allow A and C side differences
8a94851d 1894 histoRPhi->Write("SpaceChargeInRPhi");
1895 f->Close();
1896 //
1897 //
1898 //
1899 AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
1900 spaceCharge->SetSCDataFileName(fileName);
1901 spaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1902 spaceCharge->InitSpaceCharge3DDistortion();
1903 spaceCharge->AddVisualCorrection(spaceCharge,1);
1904 //
1905 f = new TFile(fileName, "update");
1906 spaceCharge->Write("spaceCharge");
1907 f->Close();
1908 return spaceCharge;
1909}
1910