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