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