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 }