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 }