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