Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /// \class AliTPCSpaceCharge
17 : /// \brief The class calculates the space point distortions due to a rotational
18 : ///
19 : /// symmetric space charge distribution with the TPC drift volume.
20 : ///
21 : /// The class uses the PoissonRelaxation2D to calculate the resulting
22 : /// electrical field inhomogeneities in the (r,z)-plane. Then, the
23 : /// Langevin-integral formalism is used to calculate the space point distortions.
24 : ///
25 : /// The class assumes, that the distortions scales linearly with the magnitude
26 : /// of the space charge distribution $\rho(r,z)$. The in here assumed distribution is
27 : /// $$\rho(r,z) = \frac{(A-B\,z)}{r^2} $$ wherein the factors A and B scale with the
28 : /// event multiplicity and the interaction rate.
29 : ///
30 : /// The scaling factor can be set via the function SetCorrectionFactor. An example of
31 : /// the shape of the distortions is given below.
32 : ///
33 : /// MI modification - 22.05.2013
34 : /// As an optional input the Space charge histogram RZ is used in case it is provided
35 : /// - using the SetInputSpaceCharge function
36 : /// 
37 : ///
38 : /// \author Jim Thomas, Stefan Rossegger
39 : /// \date 23/08/2010
40 :
41 :
42 :
43 : #include "AliMagF.h"
44 : #include "TGeoGlobalMagField.h"
45 : #include "AliTPCcalibDB.h"
46 : #include "AliTPCParam.h"
47 : #include "AliLog.h"
48 : #include "TMatrixD.h"
49 : #include "TH2.h"
50 :
51 : #include "TMath.h"
52 : #include "AliTPCROC.h"
53 : #include "AliTPCSpaceCharge.h"
54 :
55 : /// \cond CLASSIMP
56 24 : ClassImp(AliTPCSpaceCharge)
57 : /// \endcond
58 :
59 : AliTPCSpaceCharge::AliTPCSpaceCharge()
60 0 : : AliTPCCorrection("SpaceCharge2D","Space Charge 2D"),
61 0 : fC0(0.),fC1(0.),fCorrectionFactor(0.001),fSpaceChargeHistogram(0),
62 0 : fInitLookUp(kFALSE)
63 0 : {
64 : //
65 : // default constructor
66 : //
67 :
68 0 : }
69 :
70 0 : AliTPCSpaceCharge::~AliTPCSpaceCharge() {
71 : /// default destructor
72 :
73 0 : }
74 :
75 :
76 :
77 : void AliTPCSpaceCharge::Init() {
78 : /// Initialization funtion
79 :
80 0 : AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
81 0 : if (!magF) AliError("Magneticd field - not initialized");
82 0 : Double_t bzField = magF->SolenoidField()/10.; //field in T
83 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
84 0 : if (!param) AliError("Parameters - not initialized");
85 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
86 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
87 0 : Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
88 : // Correction Terms for effective omegaTau; obtained by a laser calibration run
89 0 : SetOmegaTauT1T2(wt,fT1,fT2);
90 :
91 0 : InitSpaceChargeDistortion(); // fill the look up table
92 0 : }
93 :
94 : void AliTPCSpaceCharge::Update(const TTimeStamp &/*timeStamp*/) {
95 : /// Update function
96 :
97 0 : AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
98 0 : if (!magF) AliError("Magneticd field - not initialized");
99 0 : Double_t bzField = magF->SolenoidField()/10.; //field in T
100 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
101 0 : if (!param) AliError("Parameters - not initialized");
102 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
103 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
104 0 : Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
105 : // Correction Terms for effective omegaTau; obtained by a laser calibration run
106 0 : SetOmegaTauT1T2(wt,fT1,fT2);
107 :
108 : // SetCorrectionFactor(1.); // should come from some database
109 :
110 0 : }
111 :
112 :
113 :
114 : void AliTPCSpaceCharge::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
115 : /// Calculates the correction due the Space Charge effect within the TPC drift volume
116 :
117 0 : if (!fInitLookUp) {
118 0 : AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
119 0 : InitSpaceChargeDistortion();
120 0 : }
121 : Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
122 :
123 0 : Double_t intEr, intEphi, intdEz;
124 : Double_t r, phi, z ;
125 : Int_t sign;
126 :
127 0 : r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
128 0 : phi = TMath::ATan2(x[1],x[0]) ;
129 0 : if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
130 0 : z = x[2] ; // Create temporary copy of x[2]
131 :
132 0 : if ( (roc%36) < 18 ) {
133 : sign = 1; // (TPC A side)
134 0 : } else {
135 : sign = -1; // (TPC C side)
136 : }
137 :
138 0 : if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
139 0 : if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
140 :
141 :
142 0 : if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
143 0 : AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
144 :
145 : // Efield is symmetric in phi - 2D calculation
146 : intEphi = 0.0;
147 : // Get the E field integrals
148 0 : Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr );
149 : // Get DeltaEz field integral
150 0 : Interpolate2DEdistortion( order, r, z, fLookUpDeltaEz, intdEz );
151 :
152 :
153 : // Calculate distorted position
154 0 : if ( r > 0.0 ) {
155 0 : phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
156 0 : r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
157 0 : }
158 0 : Double_t dz = intdEz*fCorrectionFactor;
159 :
160 : // Calculate correction in cartesian coordinates
161 0 : dx[0] = - (r * TMath::Cos(phi) - x[0]);
162 0 : dx[1] = - (r * TMath::Sin(phi) - x[1]);
163 0 : dx[2] = - dz; // z distortion - (internally scaled with driftvelocity dependency
164 : // on the Ez field
165 :
166 0 : }
167 :
168 : void AliTPCSpaceCharge::InitSpaceChargeDistortion() {
169 : /// Initialization of the Lookup table which contains the solutions of the
170 : /// poisson problem
171 :
172 0 : const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
173 0 : const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
174 :
175 0 : TMatrixD voltArray(kRows,kColumns); // dummy boundary vectors
176 0 : TMatrixD chargeDensity(kRows,kColumns); // charge
177 0 : TMatrixD arrayErOverEz(kRows,kColumns); // solution in Er
178 0 : TMatrixD arrayDeltaEz(kRows,kColumns); // solution in Ez
179 :
180 0 : Double_t rList[kRows], zedList[kColumns] ;
181 :
182 : // Fill arrays with initial conditions. V on the boundary and ChargeDensity in the volume.
183 0 : for ( Int_t j = 0 ; j < kColumns ; j++ ) {
184 0 : Double_t zed = j*gridSizeZ ;
185 0 : zedList[j] = zed ;
186 0 : for ( Int_t i = 0 ; i < kRows ; i++ ) {
187 0 : Double_t radius = fgkIFCRadius + i*gridSizeR ;
188 0 : rList[i] = radius ;
189 0 : voltArray(i,j) = 0; // Initialize voltArray to zero - not used in this class
190 0 : chargeDensity(i,j) = 0; // Initialize ChargeDensity to zero
191 : }
192 : }
193 :
194 : // Fill the initial conditions
195 0 : for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
196 0 : Double_t zed = j*gridSizeZ ;
197 0 : for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
198 0 : Double_t radius = fgkIFCRadius + i*gridSizeR ;
199 :
200 0 : Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadius*fgkOFCRadius - fgkIFCRadius*fgkIFCRadius) / fgkTPCZ0 ;
201 : // for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
202 0 : chargeDensity(i,j) = zterm / ( TMath::Log(fgkOFCRadius/fgkIFCRadius) * ( radius*radius ) ) ;
203 : }
204 : }
205 : // Fill the initial space charge in case histogram exist
206 0 : if (fSpaceChargeHistogram){
207 0 : for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
208 0 : Double_t zed = j*gridSizeZ ;
209 0 : for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
210 0 : Double_t radius = fgkIFCRadius + i*gridSizeR ;
211 :
212 0 : Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadius*fgkOFCRadius - fgkIFCRadius*fgkIFCRadius) / fgkTPCZ0 ;
213 : // for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
214 0 : chargeDensity(i,j) = fSpaceChargeHistogram->Interpolate(radius,zed);
215 : }
216 : }
217 0 : }
218 :
219 :
220 : // Solve the electrosatic problem in 2D
221 :
222 0 : PoissonRelaxation2D( voltArray, chargeDensity, arrayErOverEz, arrayDeltaEz, kRows, kColumns, kIterations ) ;
223 :
224 : //Interpolate results onto standard grid for Electric Fields
225 0 : Int_t ilow=0, jlow=0 ;
226 : Double_t z,r;
227 : Float_t saveEr[2], saveEz[2] ;
228 0 : for ( Int_t i = 0 ; i < kNZ ; ++i ) {
229 0 : z = TMath::Abs( fgkZList[i] ) ; // assume symmetric behaviour on A and C side
230 0 : for ( Int_t j = 0 ; j < kNR ; ++j ) {
231 :
232 : // Linear interpolation !!
233 0 : r = fgkRList[j] ;
234 0 : Search( kRows, rList, r, ilow ) ; // Note switch - R in rows and Z in columns
235 0 : Search( kColumns, zedList, z, jlow ) ;
236 0 : if ( ilow < 0 ) ilow = 0 ; // check if out of range
237 0 : if ( jlow < 0 ) jlow = 0 ;
238 0 : if ( ilow + 1 >= kRows - 1 ) ilow = kRows - 2 ;
239 0 : if ( jlow + 1 >= kColumns - 1 ) jlow = kColumns - 2 ;
240 :
241 0 : saveEr[0] = arrayErOverEz(ilow,jlow) +
242 0 : (arrayErOverEz(ilow,jlow+1)-arrayErOverEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
243 0 : saveEr[1] = arrayErOverEz(ilow+1,jlow) +
244 0 : (arrayErOverEz(ilow+1,jlow+1)-arrayErOverEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
245 0 : saveEz[0] = arrayDeltaEz(ilow,jlow) +
246 0 : (arrayDeltaEz(ilow,jlow+1)-arrayDeltaEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
247 0 : saveEz[1] = arrayDeltaEz(ilow+1,jlow) +
248 0 : (arrayDeltaEz(ilow+1,jlow+1)-arrayDeltaEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
249 :
250 :
251 0 : fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ;
252 0 : fLookUpDeltaEz[i][j] = saveEz[0] + (saveEz[1]-saveEz[0])*(r-rList[ilow])/gridSizeR ;
253 :
254 0 : if (fgkZList[i]<0) fLookUpDeltaEz[i][j] *= -1; // C side is negative z
255 : }
256 : }
257 :
258 0 : fInitLookUp = kTRUE;
259 :
260 0 : }
261 :
262 : void AliTPCSpaceCharge::Print(const Option_t* option) const {
263 : /// Print function to check the settings of the boundary vectors
264 : /// option=="a" prints the C0 and C1 coefficents for calibration purposes
265 :
266 0 : TString opt = option; opt.ToLower();
267 0 : printf("%s\n",GetTitle());
268 0 : printf(" - Space Charge effects assuming a radial symmetric z over r^2 SC-distribution.\n");
269 0 : printf(" SC correction factor: %f \n",fCorrectionFactor);
270 :
271 0 : if (opt.Contains("a")) { // Print all details
272 0 : printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
273 0 : printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
274 : }
275 :
276 0 : if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceChargeDistortion() ...");
277 :
278 0 : }
|
