1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////////
18 // AliTPCSpaceCharge3D class //
19 // The class calculates the space point distortions due to a space charge //
21 // Method of calculation: //
22 // The analytical solution for the poisson problem in 3D (cylindrical coord)//
23 // is used in form of look up tables. PieceOfCake (POC) voxel were pre- //
24 // calculated and can be sumed up (weighted) according to the what is needed//
26 // The class allows "effective Omega Tau" corrections. //
28 // NOTE: This class is not capable of calculation z distortions due to //
29 // drift velocity change in dependence of the electric field!!! //
31 // date: 19/06/2010 //
32 // Authors: Stefan Rossegger //
35 //////////////////////////////////////////////////////////////////////////////
38 #include "TGeoGlobalMagField.h"
39 #include "AliTPCcalibDB.h"
40 #include "AliTPCParam.h"
50 #include "AliTPCROC.h"
51 #include "AliTPCSpaceCharge3D.h"
53 ClassImp(AliTPCSpaceCharge3D)
55 AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
56 : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
58 fCorrectionFactor(1.),
60 fSCDataFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root"),
61 fSCLookUpPOCsFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root")
64 // default constructor
67 // Array which will contain the solution according to the setted charge density distribution
68 // see InitSpaceCharge3DDistortion() function
69 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
70 fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
71 fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
72 fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
73 fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
76 printf("%s\n",fSCDataFileName);
77 printf("%s\n",fSCLookUpPOCsFileName);
78 SetSCDataFileName(fSCDataFileName);
82 AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
87 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
88 delete fLookUpErOverEz[k];
89 delete fLookUpEphiOverEz[k];
90 delete fLookUpDeltaEz[k];
91 delete fSCdensityDistribution[k];
94 // delete fSCdensityDistribution;
99 void AliTPCSpaceCharge3D::Init() {
101 // Initialization funtion
104 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
105 if (!magF) AliError("Magneticd field - not initialized");
106 Double_t bzField = magF->SolenoidField()/10.; //field in T
107 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
108 if (!param) AliError("Parameters - not initialized");
109 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
110 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
111 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
112 // Correction Terms for effective omegaTau; obtained by a laser calibration run
113 SetOmegaTauT1T2(wt,fT1,fT2);
115 InitSpaceCharge3DDistortion(); // fill the look up table
118 void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
122 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
123 if (!magF) AliError("Magneticd field - not initialized");
124 Double_t bzField = magF->SolenoidField()/10.; //field in T
125 AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
126 if (!param) AliError("Parameters - not initialized");
127 Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
128 Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
129 Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
130 // Correction Terms for effective omegaTau; obtained by a laser calibration run
131 SetOmegaTauT1T2(wt,fT1,fT2);
133 // SetCorrectionFactor(1.); // should come from some database
139 void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
141 // Calculates the correction due the Space Charge effect within the TPC drift volume
145 AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
146 InitSpaceCharge3DDistortion();
149 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
151 Double_t intEr, intEphi, intdEz ;
155 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
156 phi = TMath::ATan2(x[1],x[0]) ;
157 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
158 z = x[2] ; // Create temporary copy of x[2]
160 if ( (roc%36) < 18 ) {
161 sign = 1; // (TPC A side)
163 sign = -1; // (TPC C side)
166 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
167 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
170 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
171 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
173 // Get the Er and Ephi field integrals plus the integral over DeltaEz
174 intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
175 fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
176 intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
177 fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
178 intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
179 fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
181 // Calculate distorted position
183 phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
184 r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
186 Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
188 // Calculate correction in cartesian coordinates
189 dx[0] = r * TMath::Cos(phi) - x[0];
190 dx[1] = r * TMath::Sin(phi) - x[1];
191 dx[2] = dz; // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
196 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
198 // Initialization of the Lookup table which contains the solutions of the
199 // "space charge" (poisson) problem
201 // The sum-up uses a look-up table which contains different discretized Space charge fields
202 // in order to calculate the corresponding field deviations due to a given (discretized)
203 // space charge distribution ....
205 // Method of calculation: Weighted sum up of the different fields within the look up table
209 AliInfo("Lookup table was already initialized!");
213 AliInfo("Preparation of the weighted look-up table");
215 TFile *f = new TFile(fSCLookUpPOCsFileName);
217 AliError("Precalculated POC-looup-table could not be found");
222 TVector *gridf = (TVector*) f->Get("constants");
223 TVector &grid = *gridf;
224 TMatrix *coordf = (TMatrix*) f->Get("coordinates");
225 TMatrix &coord = *coordf;
226 TMatrix *coordPOCf = (TMatrix*) f->Get("POCcoord");
227 TMatrix &coordPOC = *coordPOCf;
229 Bool_t flagRadSym = 0;
230 if (grid(1)==1 && grid(4)==1) {
231 AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
235 Int_t rows = (Int_t)grid(0); // number of points in r direction
236 Int_t phiSlices = (Int_t)grid(1); // number of points in phi
237 Int_t columns = (Int_t)grid(2); // number of points in z direction
239 const Float_t gridSizeR = grid(6)*100; // unit in [cm]
240 const Float_t gridSizePhi = grid(7); // unit in [rad]
241 const Float_t gridSizeZ = grid(8)*100; // unit in [cm]
243 // temporary matrices needed for the calculation
244 TMatrixD *arrayofErA[phiSlices], *arrayofEphiA[phiSlices], *arrayofdEzA[phiSlices];
245 TMatrixD *arrayofErC[phiSlices], *arrayofEphiC[phiSlices], *arrayofdEzC[phiSlices];
247 TMatrixD *arrayofEroverEzA[phiSlices], *arrayofEphioverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
248 TMatrixD *arrayofEroverEzC[phiSlices], *arrayofEphioverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
251 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
253 arrayofErA[k] = new TMatrixD(rows,columns) ;
254 arrayofEphiA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
255 arrayofdEzA[k] = new TMatrixD(rows,columns) ;
256 arrayofErC[k] = new TMatrixD(rows,columns) ;
257 arrayofEphiC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
258 arrayofdEzC[k] = new TMatrixD(rows,columns) ;
260 arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
261 arrayofEphioverEzA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
262 arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
263 arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
264 arrayofEphioverEzC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
265 arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
267 // Set the values to zero the lookup tables
268 // not necessary, it is done in the constructor of TMatrix - code deleted
272 // list of points as used in the interpolation (during sum up)
273 Double_t rlist[rows], zedlist[columns] , philist[phiSlices];
274 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
275 philist[k] = gridSizePhi * k;
276 for ( Int_t i = 0 ; i < rows ; i++ ) {
277 rlist[i] = fgkIFCRadius + i*gridSizeR ;
278 for ( Int_t j = 0 ; j < columns ; j++ ) {
279 zedlist[j] = j * gridSizeZ ;
285 // Prepare summation Vectors
286 // Int_t numPosInd = (Int_t)(grid(0)*grid(1)*grid(2)); // Number of fulcrums in the look-up table
287 // Int_t numPOC = (Int_t)(grid(3)*grid(4)*grid(5)); // Number of POC conf. in the look-up table
289 TTree *treePOC = (TTree*)f->Get("POCall");
291 TVector *bEr = 0; TVector *bEphi= 0; TVector *bEz = 0;
293 treePOC->SetBranchAddress("Er",&bEr);
294 if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
295 treePOC->SetBranchAddress("Ez",&bEz);
297 /* TBranch *b2Er = treePOC->GetBranch("Er");
299 if (!flagRadSym) b2Ephi = treePOC->GetBranch("Ephi");
300 TBranch *b2Ez = treePOC->GetBranch("Ez");
303 // Read the complete tree and do a weighted sum-up over the POC configurations
304 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
306 Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
307 Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
309 for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // loop over POC configurations in tree
311 // if (!flagRadSym) cout<<ipC<<endl;
313 treePOC->GetEntry(itreepC);
315 /* b2Er->GetEntry(itreepC);
316 if (!flagRadSym) b2Ephi->GetEntry(itreepC);
317 b2Ez->GetEntry(itreepC);
320 // center of the POC voxel in [meter]
321 Double_t r0 = coordPOC(ipC,0);
322 Double_t phi0 = coordPOC(ipC,1);
323 Double_t z0 = coordPOC(ipC,2);
325 ipC++; // POC conf. counter
327 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
328 // note: coordinates are in [cm]
329 Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
330 Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
332 // Summing up the vector components according to their weight
335 for ( Int_t j = 0 ; j < columns ; j++ ) {
336 for ( Int_t i = 0 ; i < rows ; i++ ) {
337 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
339 // check wether the coordinates were screwed
340 if ((coord(0,ip)*100-rlist[i])>0.01 ||
341 (coord(1,ip)-philist[k])>0.01 ||
342 (coord(2,ip)*100-zedlist[j])>0.01) {
343 AliError("internal error: coordinate system was screwed during the sum-up");
344 printf("lookup: (r,phi,z)=(%lf,%lf,%lf)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
345 printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist[i],philist[k],zedlist[j]);
346 AliError("Don't trust the results of the space charge calibration!");
349 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
350 // This will be the most frequent usage (hopefully)
351 // That's why we have to do this here ...
353 TMatrixD &erA = *arrayofErA[k] ;
354 TMatrixD &ephiA = *arrayofEphiA[k];
355 TMatrixD &dEzA = *arrayofdEzA[k] ;
357 TMatrixD &erC = *arrayofErC[k] ;
358 TMatrixD &ephiC = *arrayofEphiC[k];
359 TMatrixD &dEzC = *arrayofdEzC[k] ;
361 // Sum up - Efield values in [V/m] -> transition to [V/cm]
362 erA(i,j) += ((*bEr)(ip)) * weightA /100;
363 erC(i,j) += ((*bEr)(ip)) * weightC /100;
365 ephiA(i,j) += ((*bEphi)(ip)) * weightA; // [V/rad]
366 ephiC(i,j) += ((*bEphi)(ip)) * weightC; // [V/rad]
368 dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
369 dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
371 // increase the counter
377 // Rotation and summation in the rest of the dPhiSteps
378 // which were not stored in the tree due to storage & symmetry reasons
380 Int_t phiPoints = (Int_t) grid(1);
381 Int_t phiPOC = (Int_t) grid(4);
383 // printf("%d %d\n",phiPOC,flagRadSym);
385 for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
387 r0 = coordPOC(ipC,0);
388 phi0 = coordPOC(ipC,1);
389 z0 = coordPOC(ipC,2);
391 ipC++; // POC conf. counter
393 // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
394 // note: coordinates are in [cm]
395 weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
396 weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
399 // Summing up the vector components according to their weight
401 for ( Int_t j = 0 ; j < columns ; j++ ) {
402 for ( Int_t i = 0 ; i < rows ; i++ ) {
403 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
405 // Note: rotating the coordinated during the sum up
407 Int_t rotVal = (phiPoints/phiPOC)*phiiC;
410 if ((ip%phiPoints)>=rotVal) {
413 ipR = ip+(phiPoints-rotVal);
416 // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
417 // This will be the most frequent usage (hopefully)
418 // That's why we have to do this here and not outside the loop ...
420 TMatrixD &erA = *arrayofErA[k] ;
421 TMatrixD &ephiA = *arrayofEphiA[k];
422 TMatrixD &dEzA = *arrayofdEzA[k] ;
424 TMatrixD &erC = *arrayofErC[k] ;
425 TMatrixD &ephiC = *arrayofEphiC[k];
426 TMatrixD &dEzC = *arrayofdEzC[k] ;
428 // Sum up - Efield values in [V/m] -> transition to [V/cm]
429 erA(i,j) += ((*bEr)(ipR)) * weightA /100;
430 erC(i,j) += ((*bEr)(ipR)) * weightC /100;
432 ephiA(i,j) += ((*bEphi)(ipR)) * weightA; // [V/rad]
433 ephiC(i,j) += ((*bEphi)(ipR)) * weightC; // [V/rad]
435 dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
436 dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
438 // increase the counter
442 } // end coordinate loop
444 } // end phi-POC summation (phiiC)
447 // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
451 AliInfo("Division and integration");
453 // -------------------------------------------------------------------------------
454 // Division by the Ez (drift) field and integration along z
456 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
458 for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
460 // matrices holding the solution - summation of POC charges
461 TMatrixD &erA = *arrayofErA[k] ;
462 TMatrixD &ephiA = *arrayofEphiA[k];
463 TMatrixD &dezA = *arrayofdEzA[k] ;
464 TMatrixD &erC = *arrayofErC[k] ;
465 TMatrixD &ephiC = *arrayofEphiC[k];
466 TMatrixD &dezC = *arrayofdEzC[k] ;
468 // matrices which contain the integrated fields (divided by the drift field)
469 TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
470 TMatrixD &ephiOverEzA = *arrayofEphioverEzA[k];
471 TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
472 TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
473 TMatrixD &ephiOverEzC = *arrayofEphioverEzC[k];
474 TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
476 for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
477 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
478 // Count backwards to facilitate integration over Z
480 Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
482 erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0; deltaEzA(i,j) = 0;
483 erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0; deltaEzC(i,j) = 0;
485 for ( Int_t m = j ; m < columns ; m++ ) { // integration
487 erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
488 erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
490 ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
491 ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
493 deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m) ;
494 deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m) ;
496 if ( index != 4 ) index = 4; else index = 2 ;
501 erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
502 erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
504 ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
505 ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
507 deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1) ;
508 deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1) ;
511 erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
512 erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
514 ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
515 ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
517 deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1)) ;
518 deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1)) ;
520 if ( j == columns-2 ) {
521 erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
522 erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
524 ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
525 ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
527 deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1)) ;
528 deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1)) ;
530 if ( j == columns-1 ) {
531 erOverEzA(i,j) = 0; erOverEzC(i,j) = 0;
533 ephiOverEzA(i,j) = 0; ephiOverEzC(i,j) = 0;
535 deltaEzA(i,j) = 0; deltaEzC(i,j) = 0;
542 AliInfo("Interpolation to Standard grid");
544 // -------------------------------------------------------------------------------
545 // Interpolate results onto the standard grid which is used for all AliTPCCorrections
547 const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
550 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
551 phi = fgkPhiList[k] ;
553 TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
554 TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
555 TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
557 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
559 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
561 for ( Int_t i = 0 ; i < kNR ; i++ ) {
564 // Interpolate Lookup tables onto standard grid
566 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
567 rlist, zedlist, philist, arrayofEroverEzA );
568 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
569 rlist, zedlist, philist, arrayofEphioverEzA);
570 deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
571 rlist, zedlist, philist, arrayofDeltaEzA );
573 erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
574 rlist, zedlist, philist, arrayofEroverEzC );
575 ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
576 rlist, zedlist, philist, arrayofEphioverEzC);
577 deltaEz(i,j) = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
578 rlist, zedlist, philist, arrayofDeltaEzC );
579 // negative coordinate system on C side
587 // clear the temporary arrays lists
588 for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
590 delete arrayofErA[k];
591 delete arrayofEphiA[k];
592 delete arrayofdEzA[k];
593 delete arrayofErC[k];
594 delete arrayofEphiC[k];
595 delete arrayofdEzC[k];
597 delete arrayofEroverEzA[k];
598 delete arrayofEphioverEzA[k];
599 delete arrayofDeltaEzA[k];
600 delete arrayofEroverEzC[k];
601 delete arrayofEphioverEzC[k];
602 delete arrayofDeltaEzC[k];
612 void AliTPCSpaceCharge3D::SetSCDataFileName(char *const fname) {
614 // Set & load the Space charge density distribution from a file
615 // (linear interpolation onto a standard grid)
618 fSCDataFileName = fname;
620 TFile *f = new TFile(fSCDataFileName,"READ");
622 AliError(Form("File %s, which should contain the space charge distribution, could not be found",
627 TH3F *density = (TH3F*) f->Get("SpaceChargeDistribution");
629 AliError(Form("The indicated file (%s) does not contain a histogram called %s",
630 fSCDataFileName,"SpaceChargeDistribution"));
636 for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
637 phi = fgkPhiList[k] ;
638 TMatrixD &scdensity = *fSCdensityDistribution[k] ;
639 for ( Int_t j = 0 ; j < kNZ ; j++ ) {
640 z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
641 for ( Int_t i = 0 ; i < kNR ; i++ ) {
643 scdensity(i,j) = density->Interpolate(r,phi,z); // quite slow
644 // printf("%lf %lf %lf: %lf\n",r,phi,z,scdensity(i,j));
652 fInitLookUp = kFALSE;
658 Float_t AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z) {
660 // returns the (input) space charge density at a given point according
661 // Note: input in [cm], output in [C/m^3/e0] !!
664 if (!fSCdensityDistribution) {
665 printf("Scheiss is nuul - argg\n");
669 while (phi<0) phi += TMath::TwoPi();
670 while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
673 // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
675 Float_t sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
676 fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
678 // printf("%lf %lf %lf: %lf\n",r,phi,z,sc);
684 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny) {
686 // return a simple histogramm containing the space charge distribution (input for the calculation)
689 TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
690 nx,-250.,250.,ny,-250.,250.);
692 for (Int_t iy=1;iy<=ny;++iy) {
693 Double_t yp = h->GetYaxis()->GetBinCenter(iy);
694 for (Int_t ix=1;ix<=nx;++ix) {
695 Double_t xp = h->GetXaxis()->GetBinCenter(ix);
697 Float_t r = TMath::Sqrt(xp*xp+yp*yp);
698 Float_t phi = TMath::ATan2(yp,xp);
700 if (85.<=r && r<=250.) {
701 Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
702 h->SetBinContent(ix,iy,sc);
704 h->SetBinContent(ix,iy,0.);
712 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr) {
714 // return a simple histogramm containing the space charge distribution (input for the calculation)
717 TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
718 nz,-250.,250.,nr,85.,250.);
720 for (Int_t ir=1;ir<=nr;++ir) {
721 Float_t r = h->GetYaxis()->GetBinCenter(ir);
722 for (Int_t iz=1;iz<=nz;++iz) {
723 Float_t z = h->GetXaxis()->GetBinCenter(iz);
724 Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
725 h->SetBinContent(iz,ir,sc);
732 void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
734 // Example on how to write a Space charge distribution into a File
735 // (see below: estimate from scaling STAR measurements to Alice)
738 TFile *f = new TFile(fname,"RECREATE");
740 // some grid, not too course
745 Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
746 Double_t dphi = TMath::TwoPi()/(nphi+1);
747 Double_t dz = 500./(nz+1);
748 Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
750 TH3F *histo = new TH3F("charge","charge",
751 nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
752 nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
753 nz,-250-dz-safty,250+dz+safty);
755 for (Int_t ir=1;ir<=nr;++ir) {
756 Double_t rp = histo->GetXaxis()->GetBinCenter(ir);
757 for (Int_t iphi=1;iphi<=nphi;++iphi) {
758 Double_t phip = histo->GetYaxis()->GetBinCenter(iphi);
759 for (Int_t iz=1;iz<=nz;++iz) {
760 Double_t zp = histo->GetZaxis()->GetBinCenter(iz);
763 zp = TMath::Abs(zp); // symmetric on both sides of the TPC
765 // recalculation to meter
766 Double_t lZ = 2.5; // approx. TPC drift length
767 Double_t rpM = rp/100.; // in [m]
768 Double_t zpM = zp/100.; // in [m]
770 // setting of mb multiplicity and Interaction rate
771 Double_t multiplicity = 950;
772 Double_t intRate = 8000;
774 // calculation of "scaled" parameters
775 Double_t a = multiplicity*intRate/79175;
778 Double_t charge = (a - b * zpM)/(rpM*rpM); // charge in [C/m^3/e0]
780 // some 'arbitrary' GG leaks
781 if ( (phip<(TMath::Pi()/9*3) && phip>(TMath::Pi()/9*2) ) ) { //||
782 // (phip<(TMath::Pi()/9*7) && phip>(TMath::Pi()/9*6) ) ||
783 // (phip<(TMath::Pi()/9*16) && phip>(TMath::Pi()/9*13) ) ){
788 // printf("%lf %lf %lf: %lf\n",rp,phip,zp,charge);
790 charge = charge*fgke0; // [C/m^3]
792 histo->SetBinContent(ir,iphi,iz,charge);
798 histo->Write("SpaceChargeDistribution");
804 void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
806 // Print function to check the settings of the boundary vectors
807 // option=="a" prints the C0 and C1 coefficents for calibration purposes
810 TString opt = option; opt.ToLower();
811 printf("%s\n",GetTitle());
812 printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
813 printf(" SC correction factor: %f \n",fCorrectionFactor);
815 if (opt.Contains("a")) { // Print all details
816 printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
817 printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
820 if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");