]>
Commit | Line | Data |
---|---|---|
c9cbd2f2 | 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 | ////////////////////////////////////////////////////////////////////////////// | |
17 | // // | |
18 | // AliTPCSpaceCharge class // | |
19 | // The class calculates the space point distortions due to a space charge // | |
20 | // effect .... // | |
21 | // The class allows "effective Omega Tau" corrections. // | |
22 | // // | |
23 | // NOTE: This class is capable of calculating z distortions due to // | |
24 | // drift velocity change in dependence of the electric field!!! // | |
25 | // // | |
26 | // date: 23/08/2010 // | |
27 | // Authors: Jim Thomas, Stefan Rossegger // | |
28 | // // | |
29 | // Example usage: // | |
30 | ////////////////////////////////////////////////////////////////////////////// | |
31 | ||
32 | #include "AliMagF.h" | |
33 | #include "TGeoGlobalMagField.h" | |
34 | #include "AliTPCcalibDB.h" | |
35 | #include "AliTPCParam.h" | |
36 | #include "AliLog.h" | |
37 | #include "TMatrixD.h" | |
38 | ||
39 | #include "TMath.h" | |
40 | #include "AliTPCROC.h" | |
41 | #include "AliTPCSpaceCharge.h" | |
42 | ||
43 | ClassImp(AliTPCSpaceCharge) | |
44 | ||
45 | AliTPCSpaceCharge::AliTPCSpaceCharge() | |
46 | : AliTPCCorrection("SpaceCharge2D","Space Charge 2D"), | |
47 | fC0(0.),fC1(0.),fCorrectionFactor(0.001), | |
48 | fInitLookUp(kFALSE) | |
49 | { | |
50 | // | |
51 | // default constructor | |
52 | // | |
53 | ||
54 | } | |
55 | ||
56 | AliTPCSpaceCharge::~AliTPCSpaceCharge() { | |
57 | // | |
58 | // default destructor | |
59 | // | |
60 | } | |
61 | ||
62 | ||
63 | ||
64 | void AliTPCSpaceCharge::Init() { | |
65 | // | |
66 | // Initialization funtion | |
67 | // | |
68 | ||
69 | AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
70 | if (!magF) AliError("Magneticd field - not initialized"); | |
71 | Double_t bzField = magF->SolenoidField()/10.; //field in T | |
72 | AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters(); | |
73 | if (!param) AliError("Parameters - not initialized"); | |
74 | Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally) | |
75 | Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully) | |
76 | Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; | |
77 | // Correction Terms for effective omegaTau; obtained by a laser calibration run | |
78 | SetOmegaTauT1T2(wt,fT1,fT2); | |
79 | ||
80 | InitSpaceChargeDistortion(); // fill the look up table | |
81 | } | |
82 | ||
83 | void AliTPCSpaceCharge::Update(const TTimeStamp &/*timeStamp*/) { | |
84 | // | |
85 | // Update function | |
86 | // | |
87 | AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
88 | if (!magF) AliError("Magneticd field - not initialized"); | |
89 | Double_t bzField = magF->SolenoidField()/10.; //field in T | |
90 | AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters(); | |
91 | if (!param) AliError("Parameters - not initialized"); | |
92 | Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally) | |
93 | Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully) | |
94 | Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; | |
95 | // Correction Terms for effective omegaTau; obtained by a laser calibration run | |
96 | SetOmegaTauT1T2(wt,fT1,fT2); | |
97 | ||
98 | // SetCorrectionFactor(1.); // should come from some database | |
99 | ||
100 | } | |
101 | ||
102 | ||
103 | ||
104 | void AliTPCSpaceCharge::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) { | |
105 | // | |
106 | // Calculates the correction due the Space Charge effect within the TPC drift volume | |
107 | // | |
108 | ||
109 | if (!fInitLookUp) { | |
110 | AliInfo("Lookup table was not initialized! Perform the inizialisation now ..."); | |
111 | InitSpaceChargeDistortion(); | |
112 | } | |
113 | Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2 | |
114 | ||
115 | Double_t intEr, intEphi, intdEz; | |
116 | Double_t r, phi, z ; | |
117 | Int_t sign; | |
118 | ||
119 | r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ; | |
120 | phi = TMath::ATan2(x[1],x[0]) ; | |
121 | if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi | |
122 | z = x[2] ; // Create temporary copy of x[2] | |
123 | ||
124 | if ( (roc%36) < 18 ) { | |
125 | sign = 1; // (TPC A side) | |
126 | } else { | |
127 | sign = -1; // (TPC C side) | |
128 | } | |
129 | ||
130 | if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE | |
131 | if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE | |
132 | ||
133 | ||
134 | if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check | |
135 | AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); | |
136 | ||
137 | // Efield is symmetric in phi - 2D calculation | |
138 | intEphi = 0.0; | |
139 | // Get the E field integrals | |
140 | Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr ); | |
141 | // Get DeltaEz field integral | |
142 | Interpolate2DEdistortion( order, r, z, fLookUpDeltaEz, intdEz ); | |
143 | ||
144 | ||
145 | // Calculate distorted position | |
146 | if ( r > 0.0 ) { | |
147 | phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r; | |
148 | r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi ); | |
149 | } | |
150 | Double_t dz = intdEz*fCorrectionFactor; | |
151 | ||
152 | // Calculate correction in cartesian coordinates | |
153 | dx[0] = r * TMath::Cos(phi) - x[0]; | |
154 | dx[1] = r * TMath::Sin(phi) - x[1]; | |
155 | dx[2] = dz; // z distortion - (internally scaled with driftvelocity dependency | |
156 | // on the Ez field | |
157 | ||
158 | } | |
159 | ||
160 | void AliTPCSpaceCharge::InitSpaceChargeDistortion() { | |
161 | // | |
162 | // Initialization of the Lookup table which contains the solutions of the | |
163 | // poisson problem | |
164 | // | |
165 | ||
166 | const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ; | |
167 | const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ; | |
168 | ||
169 | TMatrixD voltArray(kRows,kColumns); // dummy boundary vectors | |
170 | TMatrixD chargeDensity(kRows,kColumns); // charge | |
171 | TMatrixD arrayErOverEz(kRows,kColumns); // solution in Er | |
172 | TMatrixD arrayDeltaEz(kRows,kColumns); // solution in Ez | |
173 | ||
174 | Double_t rList[kRows], zedList[kColumns] ; | |
175 | ||
176 | // Fill arrays with initial conditions. V on the boundary and ChargeDensity in the volume. | |
177 | for ( Int_t j = 0 ; j < kColumns ; j++ ) { | |
178 | Double_t zed = j*gridSizeZ ; | |
179 | zedList[j] = zed ; | |
180 | for ( Int_t i = 0 ; i < kRows ; i++ ) { | |
181 | Double_t radius = fgkIFCRadius + i*gridSizeR ; | |
182 | rList[i] = radius ; | |
183 | voltArray(i,j) = 0; // Initialize voltArray to zero - not used in this class | |
184 | chargeDensity(i,j) = 0; // Initialize ChargeDensity to zero | |
185 | } | |
186 | } | |
187 | ||
188 | // Fill the initial conditions | |
189 | for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) { | |
190 | Double_t zed = j*gridSizeZ ; | |
191 | for ( Int_t i = 1 ; i < kRows-1 ; i++ ) { | |
192 | Double_t radius = fgkIFCRadius + i*gridSizeR ; | |
193 | ||
194 | Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadius*fgkOFCRadius - fgkIFCRadius*fgkIFCRadius) / fgkTPCZ0 ; | |
195 | // for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions | |
196 | chargeDensity(i,j) = zterm / ( TMath::Log(fgkOFCRadius/fgkIFCRadius) * ( radius*radius ) ) ; | |
197 | } | |
198 | } | |
199 | ||
200 | ||
201 | // Solve the electrosatic problem in 2D | |
202 | ||
203 | PoissonRelaxation2D( voltArray, chargeDensity, arrayErOverEz, arrayDeltaEz, kRows, kColumns, kIterations ) ; | |
204 | ||
205 | //Interpolate results onto standard grid for Electric Fields | |
206 | Int_t ilow=0, jlow=0 ; | |
207 | Double_t z,r; | |
208 | Float_t saveEr[2], saveEz[2] ; | |
209 | for ( Int_t i = 0 ; i < kNZ ; ++i ) { | |
210 | z = TMath::Abs( fgkZList[i] ) ; // assume symmetric behaviour on A and C side | |
211 | for ( Int_t j = 0 ; j < kNR ; ++j ) { | |
212 | ||
213 | // Linear interpolation !! | |
214 | r = fgkRList[j] ; | |
215 | Search( kRows, rList, r, ilow ) ; // Note switch - R in rows and Z in columns | |
216 | Search( kColumns, zedList, z, jlow ) ; | |
217 | if ( ilow < 0 ) ilow = 0 ; // check if out of range | |
218 | if ( jlow < 0 ) jlow = 0 ; | |
219 | if ( ilow + 1 >= kRows - 1 ) ilow = kRows - 2 ; | |
220 | if ( jlow + 1 >= kColumns - 1 ) jlow = kColumns - 2 ; | |
221 | ||
222 | saveEr[0] = arrayErOverEz(ilow,jlow) + | |
223 | (arrayErOverEz(ilow,jlow+1)-arrayErOverEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ; | |
224 | saveEr[1] = arrayErOverEz(ilow+1,jlow) + | |
225 | (arrayErOverEz(ilow+1,jlow+1)-arrayErOverEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ; | |
226 | saveEz[0] = arrayDeltaEz(ilow,jlow) + | |
227 | (arrayDeltaEz(ilow,jlow+1)-arrayDeltaEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ; | |
228 | saveEz[1] = arrayDeltaEz(ilow+1,jlow) + | |
229 | (arrayDeltaEz(ilow+1,jlow+1)-arrayDeltaEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ; | |
230 | ||
231 | ||
232 | fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ; | |
233 | fLookUpDeltaEz[i][j] = saveEz[0] + (saveEz[1]-saveEz[0])*(r-rList[ilow])/gridSizeR ; | |
234 | } | |
235 | } | |
236 | ||
237 | fInitLookUp = kTRUE; | |
238 | ||
239 | } | |
240 | ||
241 | void AliTPCSpaceCharge::Print(const Option_t* option) const { | |
242 | // | |
243 | // Print function to check the settings of the boundary vectors | |
244 | // option=="a" prints the C0 and C1 coefficents for calibration purposes | |
245 | // | |
246 | ||
247 | TString opt = option; opt.ToLower(); | |
248 | printf("%s\n",GetTitle()); | |
249 | printf(" - Space Charge effects assuming a radial symmetric z over r^2 SC-distribution.\n"); | |
250 | printf(" SC correction factor: %f \n",fCorrectionFactor); | |
251 | ||
252 | if (opt.Contains("a")) { // Print all details | |
253 | printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2); | |
254 | printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0); | |
255 | } | |
256 | ||
257 | if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceChargeDistortion() ..."); | |
258 | ||
259 | } |