15 minutes interval for calculation of drift correction
[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 // 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 }