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       if (fgkZList[i]<0)  fLookUpDeltaEz[i][j] *= -1; // C side is negative z
 236     }
 237   }
 238
 239   fInitLookUp = kTRUE;
 240
 241 }
 242
 243 void AliTPCSpaceCharge::Print(const Option_t* option) const {
 244   //
 245   // Print function to check the settings of the boundary vectors
 246   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 247   //
 248
 249   TString opt = option; opt.ToLower();
 250   printf("%s\n",GetTitle());
 251   printf(" - Space Charge effects assuming a radial symmetric z over r^2 SC-distribution.\n");
 252   printf("   SC correction factor: %f \n",fCorrectionFactor);
 253
 254   if (opt.Contains("a")) { // Print all details
 255     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 256     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
 257   }
 258
 259   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceChargeDistortion() ...");
 260
 261 }