TPC/AliTPCSpaceCharge3D.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>  AliTPCSpaceCharge3D class   </h2>
  20 //   The class calculates the space point distortions due to an arbitrary space
  21 //   charge distribution in 3D.
  22 //   <p>
  23 //   The method of calculation is based on the analytical solution for the Poisson
  24 //   problem in 3D (cylindrical coordinates). The solution is used in form of
  25 //   look up tables, where the pre calculated solutions for different voxel
  26 //   positions are stored. These voxel solutions can be summed up according
  27 //   to the weight of the position of the applied space charge distribution.
  28 //   Further details can be found in \cite[chap.5]{PhD-thesis_S.Rossegger}.
  29 //   <p>
  30 //   The class also allows a simple scaling of the resulting distortions
  31 //   via the function SetCorrectionFactor. This speeds up the calculation
  32 //   time under the assumption, that the distortions scales linearly with the
  33 //   magnitude of the space charge distribution $\rho(r,z)$ and the shape stays
  34 //   the same at higher luminosities.
  35 //   <p>
  36 //   In contrast to the implementation in 2D (see the class AliTPCSpaceChargeabove),
  37 //   the input charge distribution can be of arbitrary character. An example on how
  38 //   to produce a corresponding charge distribution can be found in the function
  39 //   WriteChargeDistributionToFile. In there, a $\rho(r,z) = (A-B\,z)/r^2$,
  40 //   with slightly different magnitude on the A and C side (due to the muon absorber),
  41 //   is superpositioned with a few leaking wires at arbitrary positions.
  42 // End_Html
  43 //
  44 // Begin_Macro(source)
  45 //   {
  46 //   gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
  47 //   TCanvas *c2 = new TCanvas("c2","c2",500,400);
  48 //   AliTPCSpaceCharge3D sc;
  49 //   sc.WriteChargeDistributionToFile("SC_zr2_GGleaks.root");
  50 //   sc.SetSCDataFileName("SC_zr2_GGleaks.root");
  51 //   sc.SetOmegaTauT1T2(0,1,1); // B=0
  52 //   sc.InitSpaceCharge3DDistortion();
  53 //   sc.CreateHistoDRinXY(15,300,300)->Draw("colz");
  54 //   return c2;
  55 //   }
  56 // End_Macro
  57 //
  58 // Begin_Html
  59 //   <p>
  60 //   Date: 19/06/2010  <br>
  61 //   Authors: Stefan Rossegger
  62 // End_Html
  63 // _________________________________________________________________
  64
  65
  66 #include "AliMagF.h"
  67 #include "TGeoGlobalMagField.h"
  68 #include "AliTPCcalibDB.h"
  69 #include "AliTPCParam.h"
  70 #include "AliLog.h"
  71 #include "TH2F.h"
  72 #include "TH3F.h"
  73 #include "TFile.h"
  74 #include "TVector.h"
  75 #include "TMatrix.h"
  76 #include "TMatrixD.h"
  77
  78 #include "TMath.h"
  79 #include "AliTPCROC.h"
  80 #include "AliTPCSpaceCharge3D.h"
  81
  82 ClassImp(AliTPCSpaceCharge3D)
  83
  84 AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
  85   : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
  86     fC0(0.),fC1(0.),
  87     fCorrectionFactor(1.),
  88     fInitLookUp(kFALSE),
  89     fSCDataFileName(""),
  90     fSCLookUpPOCsFileName3D(""),
  91     fSCLookUpPOCsFileNameRZ(""),
  92     fSCLookUpPOCsFileNameRPhi(""),
  93     fSCdensityInRZ(0),
  94     fSCdensityInRPhiA(0),
  95     fSCdensityInRPhiC(0)
  96 {
  97   //
  98   // default constructor
  99   //
 100
 101   // Array which will contain the solution according to the setted charge density distribution
 102   // see InitSpaceCharge3DDistortion() function
 103   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 104     fLookUpErOverEz[k]   =  new TMatrixD(kNR,kNZ);
 105     fLookUpEphiOverEz[k] =  new TMatrixD(kNR,kNZ);
 106     fLookUpDeltaEz[k]    =  new TMatrixD(kNR,kNZ);
 107     fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
 108   }
 109   fSCdensityInRZ   = new TMatrixD(kNR,kNZ);
 110   fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
 111   fSCdensityInRPhiC = new TMatrixD(kNR,kNPhi);
 112
 113   // location of the precalculated look up tables
 114
 115   fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
 116   fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/Calib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
 117   fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
 118   //  fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
 119
 120
 121
 122   // standard location of the space charge distibution ... can be changes
 123   fSCDataFileName="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root";
 124
 125   //  SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
 126
 127
 128 }
 129
 130 AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
 131   //
 132   // default destructor
 133   //
 134
 135   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 136     delete fLookUpErOverEz[k];
 137     delete fLookUpEphiOverEz[k];
 138     delete fLookUpDeltaEz[k];
 139     delete fSCdensityDistribution[k];
 140   }
 141   delete fSCdensityInRZ;
 142   delete fSCdensityInRPhiA;
 143   delete fSCdensityInRPhiC;
 144
 145 }
 146
 147
 148 void AliTPCSpaceCharge3D::Init() {
 149   //
 150   // Initialization funtion
 151   //
 152
 153   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 154   if (!magF) AliError("Magneticd field - not initialized");
 155   Double_t bzField = magF->SolenoidField()/10.; //field in T
 156   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 157   if (!param) AliError("Parameters - not initialized");
 158   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 159   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 160   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 161   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 162   SetOmegaTauT1T2(wt,fT1,fT2);
 163
 164   InitSpaceCharge3DDistortion(); // fill the look up table
 165 }
 166
 167 void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
 168   //
 169   // Update function
 170   //
 171   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 172   if (!magF) AliError("Magneticd field - not initialized");
 173   Double_t bzField = magF->SolenoidField()/10.; //field in T
 174   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 175   if (!param) AliError("Parameters - not initialized");
 176   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
 177   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 178   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 179   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 180   SetOmegaTauT1T2(wt,fT1,fT2);
 181
 182   //  SetCorrectionFactor(1.); // should come from some database
 183
 184 }
 185
 186
 187 void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 188   //
 189   // Calculates the correction due the Space Charge effect within the TPC drift volume
 190   //
 191
 192   if (!fInitLookUp) {
 193     AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
 194     InitSpaceCharge3DDistortion();
 195   }
 196
 197   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 198
 199   Double_t intEr, intEphi, intdEz ;
 200   Double_t r, phi, z ;
 201   Int_t    sign;
 202
 203   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 204   phi    =  TMath::ATan2(x[1],x[0]) ;
 205   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 206   z      =  x[2] ;                                         // Create temporary copy of x[2]
 207
 208   if ( (roc%36) < 18 ) {
 209     sign =  1;       // (TPC A side)
 210   } else {
 211     sign = -1;       // (TPC C side)
 212   }
 213
 214   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 215   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 216
 217
 218   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 219     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 220
 221   // Get the Er and Ephi field integrals plus the integral over DeltaEz
 222   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 223                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
 224   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 225                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
 226   intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 227                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );
 228
 229   // Calculate distorted position
 230   if ( r > 0.0 ) {
 231     phi =  phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
 232     r   =  r   + fCorrectionFactor *( fC0*intEr   + fC1*intEphi );
 233   }
 234   Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
 235
 236   // Calculate correction in cartesian coordinates
 237   dx[0] = - (r * TMath::Cos(phi) - x[0]);
 238   dx[1] = - (r * TMath::Sin(phi) - x[1]);
 239   dx[2] = - dz;  // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
 240
 241 }
 242
 243 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
 244  //
 245   // Initialization of the Lookup table which contains the solutions of the
 246   // "space charge" (poisson) problem - Faster and more accureate
 247   //
 248   // Method: Weighted sum-up of the different fields within the look up table
 249   // but using two lookup tables with higher granularity in the (r,z) and the (rphi)- plane to emulate
 250   // more realistic space charges. (r,z) from primary ionisation. (rphi) for possible Gating leaks
 251
 252   if (fInitLookUp) {
 253     AliInfo("Lookup table was already initialized!  Doing it again anyway ...");
 254     //    return;
 255   }
 256
 257   // ------------------------------------------------------------------------------------------------------
 258   // step 1: lookup table in rz, fine grid, radial symetric, to emulate primary ionization
 259
 260   AliInfo("Step 1: Preparation of the weighted look-up tables.");
 261
 262   // lookup table in rz, fine grid
 263
 264   TFile *fZR = new TFile(fSCLookUpPOCsFileNameRZ.Data(),"READ");
 265   if ( !fZR ) {
 266     AliError("Precalculated POC-looup-table in ZR could not be found");
 267     return;
 268   }
 269
 270   // units are in [m]
 271   TVector *gridf1  = (TVector*) fZR->Get("constants");
 272   TVector &grid1 = *gridf1;
 273   TMatrix *coordf1  = (TMatrix*) fZR->Get("coordinates");
 274   TMatrix &coord1 = *coordf1;
 275   TMatrix *coordPOCf1  = (TMatrix*) fZR->Get("POCcoord");
 276   TMatrix &coordPOC1 = *coordPOCf1;
 277
 278   Int_t rows      = (Int_t)grid1(0);   // number of points in r direction  - from RZ or RPhi table
 279   Int_t phiSlices = (Int_t)grid1(1);   // number of points in phi          - from RPhi table
 280   Int_t columns   = (Int_t)grid1(2);   // number of points in z direction  - from RZ table
 281
 282   Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 283   Float_t gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 284
 285   // temporary matrices needed for the calculation // for rotational symmetric RZ table, phislices is 1
 286
 287   TMatrixD *arrayofErA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
 288   TMatrixD *arrayofErC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
 289
 290   TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
 291   TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
 292
 293
 294   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 295
 296     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 297     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 298     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 299     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 300
 301     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 302     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 303     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 304     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 305
 306     // zero initialization not necessary, it is done in the constructor of TMatrix
 307
 308   }
 309
 310   // list of points as used during sum up
 311   Double_t  rlist1[kNRows], zedlist1[kNColumns];// , philist1[phiSlices];
 312   for ( Int_t i = 0 ; i < rows ; i++ )    {
 313     rlist1[i] = fgkIFCRadius + i*gridSizeR ;
 314     for ( Int_t j = 0 ; j < columns ; j++ ) {
 315       zedlist1[j]  = j * gridSizeZ ;
 316     }
 317   }
 318
 319   TTree *treePOC = (TTree*)fZR->Get("POCall");
 320
 321   TVector *bEr  = 0;   //TVector *bEphi= 0;
 322   TVector *bEz  = 0;
 323
 324   treePOC->SetBranchAddress("Er",&bEr);
 325   treePOC->SetBranchAddress("Ez",&bEz);
 326
 327
 328   // Read the complete tree and do a weighted sum-up over the POC configurations
 329   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 330
 331   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 332   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 333
 334   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 335
 336     treePOC->GetEntry(itreepC);
 337
 338     // center of the POC voxel in [meter]
 339     Double_t r0 = coordPOC1(ipC,0);
 340     Double_t phi0 = coordPOC1(ipC,1);
 341     Double_t z0 = coordPOC1(ipC,2);
 342
 343     ipC++; // POC configuration counter
 344
 345     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 346     // note: coordinates are in [cm]
 347     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100, 1);  // partial load in r,z
 348     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100, 1);  // partial load in r,z
 349
 350     // Summing up the vector components according to their weight
 351
 352     Int_t ip = 0;
 353     for ( Int_t j = 0 ; j < columns ; j++ ) {
 354       for ( Int_t i = 0 ; i < rows ; i++ )    {
 355         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 356
 357           // check wether the coordinates were screwed
 358           if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
 359               TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
 360             AliError("internal error: coordinate system was screwed during the sum-up");
 361             printf("sum-up: (r,z)=(%f,%f)\n",rlist1[i],zedlist1[j]);
 362             printf("lookup: (r,z)=(%f,%f)\n",coord1(0,ip)*100,coord1(2,ip)*100);
 363             AliError("Don't trust the results of the space charge calculation!");
 364           }
 365
 366           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 367           // This will be the most frequent usage (hopefully)
 368           // That's why we have to do this here ...
 369
 370           TMatrixD &erA   =  *arrayofErA[k]  ;
 371           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
 372
 373           TMatrixD &erC   =  *arrayofErC[k]  ;
 374           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
 375
 376           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 377           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 378           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 379           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 380           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 381
 382           // increase the counter
 383           ip++;
 384         }
 385       }
 386     } // end coordinate loop
 387   } // end POC loop
 388
 389
 390   // -------------------------------------------------------------------------------
 391   // Division by the Ez (drift) field and integration along z
 392
 393   //  AliInfo("Step 1: Division and integration");
 394
 395   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
 396
 397   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 398
 399     // matrices holding the solution - summation of POC charges // see above
 400     TMatrixD &erA   =  *arrayofErA[k]  ;
 401     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
 402     TMatrixD &erC   =  *arrayofErC[k]  ;
 403     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
 404
 405     // matrices which will contain the integrated fields (divided by the drift field)
 406     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
 407     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
 408     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
 409     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
 410
 411     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 412       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 413         // Count backwards to facilitate integration over Z
 414
 415         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
 416                           // by trapezoidal rule.
 417
 418         erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
 419         deltaEzA(i,j) = 0;
 420         erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
 421         deltaEzC(i,j) = 0;
 422
 423         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 424
 425           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 426           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 427           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
 428           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
 429
 430           if ( index != 4 )  index = 4; else index = 2 ;
 431
 432         }
 433
 434         if ( index == 4 ) {
 435           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 436           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 437           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
 438           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
 439         }
 440         if ( index == 2 ) {
 441           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 442           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 443           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
 444           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
 445         }
 446         if ( j == columns-2 ) {
 447           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 448           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 449           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
 450           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
 451         }
 452         if ( j == columns-1 ) {
 453           erOverEzA(i,j)   = 0;
 454           erOverEzC(i,j)   = 0;
 455           deltaEzA(i,j)    = 0;
 456           deltaEzC(i,j)    = 0;
 457         }
 458       }
 459     }
 460
 461   }
 462
 463   //  AliInfo("Step 1: Interpolation to Standard grid");
 464
 465   // -------------------------------------------------------------------------------
 466   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
 467
 468   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
 469
 470   Double_t  r, z;//phi, z ;
 471   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 472     //    phi = fgkPhiList[k] ;
 473
 474     // final lookup table
 475     TMatrixD &erOverEzFinal   =  *fLookUpErOverEz[k]  ;
 476     TMatrixD &deltaEzFinal    =  *fLookUpDeltaEz[k]   ;
 477
 478     // calculated and integrated tables - just one phi slice
 479     TMatrixD &erOverEzA   =  *arrayofEroverEzA[0]  ;
 480     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[0];
 481     TMatrixD &erOverEzC   =  *arrayofEroverEzC[0]  ;
 482     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[0];
 483
 484
 485     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 486
 487       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 488
 489       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 490         r = fgkRList[i] ;
 491
 492         // Interpolate Lookup tables onto standard grid
 493         if (fgkZList[j]>0) {
 494           erOverEzFinal(i,j)   = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzA  );
 495           deltaEzFinal(i,j)    = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzA   );
 496         } else {
 497           erOverEzFinal(i,j)   = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzC  );
 498           deltaEzFinal(i,j)    = - Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzC   );
 499           // negative coordinate system on C side
 500         }
 501
 502       } // end r loop
 503     } // end z loop
 504   } // end phi loop
 505
 506
 507   // clear the temporary arrays lists
 508   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 509
 510     delete arrayofErA[k];
 511     delete arrayofdEzA[k];
 512     delete arrayofErC[k];
 513     delete arrayofdEzC[k];
 514
 515     delete arrayofEroverEzA[k];
 516     delete arrayofDeltaEzA[k];
 517     delete arrayofEroverEzC[k];
 518     delete arrayofDeltaEzC[k];
 519
 520   }
 521
 522   fZR->Close();
 523
 524   // ------------------------------------------------------------------------------------------------------
 525   // Step 2: Load and sum up lookup table in rphi, fine grid, to emulate for example a GG leak
 526
 527   //  AliInfo("Step 2: Preparation of the weighted look-up table");
 528
 529   TFile *fRPhi = new TFile(fSCLookUpPOCsFileNameRPhi.Data(),"READ");
 530   if ( !fRPhi ) {
 531     AliError("Precalculated POC-looup-table in RPhi could not be found");
 532     return;
 533   }
 534
 535   // units are in [m]
 536   TVector *gridf2  = (TVector*) fRPhi->Get("constants");
 537   TVector &grid2 = *gridf2;
 538   TMatrix *coordf2  = (TMatrix*) fRPhi->Get("coordinates");
 539   TMatrix &coord2 = *coordf2;
 540   TMatrix *coordPOCf2  = (TMatrix*) fRPhi->Get("POCcoord");
 541   TMatrix &coordPOC2 = *coordPOCf2;
 542
 543   rows      = (Int_t)grid2(0);   // number of points in r direction
 544   phiSlices = (Int_t)grid2(1);   // number of points in phi
 545   columns   = (Int_t)grid2(2);   // number of points in z direction
 546
 547   gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 548   Float_t gridSizePhi =  TMath::TwoPi()/phiSlices;         // unit in [rad]
 549   gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 550
 551   // list of points as used during sum up
 552   Double_t  rlist2[kNRows], philist2[kNPhiSlices], zedlist2[kNColumns];
 553   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 554     philist2[k] =  gridSizePhi * k;
 555     for ( Int_t i = 0 ; i < rows ; i++ )    {
 556       rlist2[i] = fgkIFCRadius + i*gridSizeR ;
 557       for ( Int_t j = 0 ; j < columns ; j++ ) {
 558         zedlist2[j]  = j * gridSizeZ ;
 559       }
 560     }
 561   } // only done once
 562
 563   // temporary matrices needed for the calculation
 564
 565   TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
 566   TMatrixD *arrayofErC2[kNPhiSlices], *arrayofEphiC2[kNPhiSlices], *arrayofdEzC2[kNPhiSlices];
 567
 568   TMatrixD *arrayofEroverEzA2[kNPhiSlices], *arrayofEphioverEzA2[kNPhiSlices], *arrayofDeltaEzA2[kNPhiSlices];
 569   TMatrixD *arrayofEroverEzC2[kNPhiSlices], *arrayofEphioverEzC2[kNPhiSlices], *arrayofDeltaEzC2[kNPhiSlices];
 570
 571
 572   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 573
 574     arrayofErA2[k]   =   new TMatrixD(rows,columns) ;
 575     arrayofEphiA2[k] =   new TMatrixD(rows,columns) ;
 576     arrayofdEzA2[k]  =   new TMatrixD(rows,columns) ;
 577     arrayofErC2[k]   =   new TMatrixD(rows,columns) ;
 578     arrayofEphiC2[k] =   new TMatrixD(rows,columns) ;
 579     arrayofdEzC2[k]  =   new TMatrixD(rows,columns) ;
 580
 581     arrayofEroverEzA2[k]   =   new TMatrixD(rows,columns) ;
 582     arrayofEphioverEzA2[k] =   new TMatrixD(rows,columns) ;
 583     arrayofDeltaEzA2[k]    =   new TMatrixD(rows,columns) ;
 584     arrayofEroverEzC2[k]   =   new TMatrixD(rows,columns) ;
 585     arrayofEphioverEzC2[k] =   new TMatrixD(rows,columns) ;
 586     arrayofDeltaEzC2[k]    =   new TMatrixD(rows,columns) ;
 587
 588     // zero initialization not necessary, it is done in the constructor of TMatrix
 589
 590   }
 591
 592
 593   treePOC = (TTree*)fRPhi->Get("POCall");
 594
 595   //  TVector *bEr  = 0;   // done above
 596   TVector *bEphi= 0;
 597   //  TVector *bEz  = 0;   // done above
 598
 599   treePOC->SetBranchAddress("Er",&bEr);
 600   treePOC->SetBranchAddress("Ephi",&bEphi);
 601   treePOC->SetBranchAddress("Ez",&bEz);
 602
 603   // Read the complete tree and do a weighted sum-up over the POC configurations
 604   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 605
 606   treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 607   ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 608
 609   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
 610
 611     treePOC->GetEntry(itreepC);
 612
 613     // center of the POC voxel in [meter]
 614     Double_t r0 = coordPOC2(ipC,0);
 615     Double_t phi0 = coordPOC2(ipC,1);
 616     //    Double_t z0 = coordPOC2(ipC,2);
 617
 618      // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 619     // note: coordinates are in [cm]
 620     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, 0.499, 2);  // partial load in r,phi
 621     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-0.499, 2);  // partial load in r,phi
 622
 623     //    printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
 624
 625     // Summing up the vector components according to their weight
 626
 627     Int_t ip = 0;
 628     for ( Int_t j = 0 ; j < columns ; j++ ) {
 629       for ( Int_t i = 0 ; i < rows ; i++ )    {
 630         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 631
 632           // check wether the coordinates were screwed
 633           if (TMath::Abs((coord2(0,ip)*100-rlist2[i]))>1 ||
 634               TMath::Abs((coord2(1,ip)-philist2[k]))>1 ||
 635               TMath::Abs((coord2(2,ip)*100-zedlist2[j]))>1) {
 636             AliError("internal error: coordinate system was screwed during the sum-up");
 637             printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
 638             printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist2[i],philist2[k],zedlist2[j]);
 639             AliError("Don't trust the results of the space charge calculation!");
 640           }
 641
 642           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 643           // This will be the most frequent usage (hopefully)
 644           // That's why we have to do this here ...
 645
 646           TMatrixD &erA   =  *arrayofErA2[k]  ;
 647           TMatrixD &ephiA =  *arrayofEphiA2[k];
 648           TMatrixD &dEzA  =  *arrayofdEzA2[k] ;
 649
 650           TMatrixD &erC   =  *arrayofErC2[k]  ;
 651           TMatrixD &ephiC =  *arrayofEphiC2[k];
 652           TMatrixD &dEzC  =  *arrayofdEzC2[k]   ;
 653
 654           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 655           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 656           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 657           ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
 658           ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
 659           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 660           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 661
 662           // increase the counter
 663           ip++;
 664         }
 665       }
 666     } // end coordinate loop
 667
 668
 669     // Rotation and summation in the rest of the dPhiSteps
 670     // which were not stored in the this tree due to storage & symmetry reasons
 671
 672
 673     Int_t phiPoints = (Int_t) grid2(1);
 674     Int_t phiPOC    = (Int_t) grid2(4);
 675
 676     //   printf("%d %d\n",phiPOC,flagRadSym);
 677
 678     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
 679
 680       Double_t phi0R = phiiC*phi0*2 + phi0; // rotate further
 681
 682       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 683       // note: coordinates are in [cm] // ecxept z
 684       weightA = GetSpaceChargeDensity(r0*100,phi0R, 0.499, 2);  // partial load in r,phi
 685       weightC = GetSpaceChargeDensity(r0*100,phi0R,-0.499, 2);  // partial load in r,phi
 686
 687       // printf("%f %f : %e %e\n",r0,phi0R,weightA,weightC);
 688
 689       // Summing up the vector components according to their weight
 690       ip = 0;
 691       for ( Int_t j = 0 ; j < columns ; j++ ) {
 692         for ( Int_t i = 0 ; i < rows ; i++ )    {
 693           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 694
 695             // Note: rotating the coordinated during the sum up
 696
 697             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
 698             Int_t ipR = -1;
 699
 700             if ((ip%phiPoints)>=rotVal) {
 701               ipR = ip-rotVal;
 702             } else {
 703               ipR = ip+(phiPoints-rotVal);
 704             }
 705
 706             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 707             // This will be the most frequent usage
 708             // That's why we have to do this here and not outside the loop ...
 709
 710             TMatrixD &erA   =  *arrayofErA2[k]  ;
 711             TMatrixD &ephiA =  *arrayofEphiA2[k];
 712             TMatrixD &dEzA  =  *arrayofdEzA2[k]   ;
 713
 714             TMatrixD &erC   =  *arrayofErC2[k]  ;
 715             TMatrixD &ephiC =  *arrayofEphiC2[k];
 716             TMatrixD &dEzC  =  *arrayofdEzC2[k]   ;
 717
 718             // Sum up - Efield values in [V/m] -> transition to [V/cm]
 719             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
 720             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
 721             ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
 722             ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
 723             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
 724             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
 725
 726             // increase the counter
 727             ip++;
 728           }
 729         }
 730       } // end coordinate loop
 731
 732     } // end phi-POC summation (phiiC)
 733
 734     ipC++; // POC configuration counter
 735
 736     //   printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
 737
 738   }
 739
 740
 741
 742
 743   // -------------------------------------------------------------------------------
 744   // Division by the Ez (drift) field and integration along z
 745
 746   //  AliInfo("Step 2: Division and integration");
 747
 748
 749   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 750
 751     // matrices holding the solution - summation of POC charges // see above
 752     TMatrixD &erA   =  *arrayofErA2[k]  ;
 753     TMatrixD &ephiA =  *arrayofEphiA2[k];
 754     TMatrixD &dezA  =  *arrayofdEzA2[k]   ;
 755     TMatrixD &erC   =  *arrayofErC2[k]  ;
 756     TMatrixD &ephiC =  *arrayofEphiC2[k];
 757     TMatrixD &dezC  =  *arrayofdEzC2[k]   ;
 758
 759     // matrices which will contain the integrated fields (divided by the drift field)
 760     TMatrixD &erOverEzA   =  *arrayofEroverEzA2[k]  ;
 761     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA2[k];
 762     TMatrixD &deltaEzA    =  *arrayofDeltaEzA2[k];
 763     TMatrixD &erOverEzC   =  *arrayofEroverEzC2[k]  ;
 764     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC2[k];
 765     TMatrixD &deltaEzC    =  *arrayofDeltaEzC2[k];
 766
 767     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 768       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 769         // Count backwards to facilitate integration over Z
 770
 771         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
 772
 773         erOverEzA(i,j) = 0;
 774         ephiOverEzA(i,j) = 0;
 775         deltaEzA(i,j) = 0;
 776         erOverEzC(i,j) = 0;
 777         ephiOverEzC(i,j) = 0;
 778         deltaEzC(i,j) = 0;
 779
 780         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 781
 782           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 783           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 784           ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
 785           ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
 786           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
 787           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
 788
 789           if ( index != 4 )  index = 4; else index = 2 ;
 790
 791         }
 792
 793         if ( index == 4 ) {
 794           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 795           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 796           ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
 797           ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
 798           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
 799           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
 800         }
 801         if ( index == 2 ) {
 802           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 803           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 804           ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
 805           ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
 806           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
 807           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
 808         }
 809         if ( j == columns-2 ) {
 810           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 811           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 812           ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
 813           ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
 814           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
 815           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
 816         }
 817         if ( j == columns-1 ) {
 818           erOverEzA(i,j)   = 0;
 819           erOverEzC(i,j)   = 0;
 820           ephiOverEzA(i,j) = 0;
 821           ephiOverEzC(i,j) = 0;
 822           deltaEzA(i,j)    = 0;
 823           deltaEzC(i,j)    = 0;
 824         }
 825       }
 826     }
 827
 828   }
 829
 830   AliInfo("Step 2: Interpolation to Standard grid");
 831
 832   // -------------------------------------------------------------------------------
 833   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
 834
 835
 836   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 837     Double_t phi = fgkPhiList[k] ;
 838
 839     // final lookup table
 840     TMatrixD &erOverEzFinal   =  *fLookUpErOverEz[k]  ;
 841     TMatrixD &ephiOverEzFinal =  *fLookUpEphiOverEz[k];
 842     TMatrixD &deltaEzFinal    =  *fLookUpDeltaEz[k]   ;
 843
 844     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 845
 846       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 847
 848       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 849         r = fgkRList[i] ;
 850
 851         // Interpolate Lookup tables onto standard grid
 852         if (fgkZList[j]>0) {
 853           erOverEzFinal(i,j)   += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 854                                                rlist2, zedlist2, philist2, arrayofEroverEzA2  );
 855           ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 856                                                rlist2, zedlist2, philist2, arrayofEphioverEzA2);
 857           deltaEzFinal(i,j)    += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 858                                                rlist2, zedlist2, philist2, arrayofDeltaEzA2   );
 859         } else {
 860           erOverEzFinal(i,j)   += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 861                                                rlist2, zedlist2, philist2, arrayofEroverEzC2  );
 862           ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 863                                                rlist2, zedlist2, philist2, arrayofEphioverEzC2);
 864           deltaEzFinal(i,j)  -=  Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 865                                                rlist2, zedlist2, philist2, arrayofDeltaEzC2   );
 866         }
 867
 868       } // end r loop
 869     } // end z loop
 870   } // end phi loop
 871
 872
 873   // clear the temporary arrays lists
 874   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 875
 876     delete arrayofErA2[k];
 877     delete arrayofEphiA2[k];
 878     delete arrayofdEzA2[k];
 879     delete arrayofErC2[k];
 880     delete arrayofEphiC2[k];
 881     delete arrayofdEzC2[k];
 882
 883     delete arrayofEroverEzA2[k];
 884     delete arrayofEphioverEzA2[k];
 885     delete arrayofDeltaEzA2[k];
 886     delete arrayofEroverEzC2[k];
 887     delete arrayofEphioverEzC2[k];
 888     delete arrayofDeltaEzC2[k];
 889
 890   }
 891
 892   fRPhi->Close();
 893
 894   // FINISHED
 895
 896   fInitLookUp = kTRUE;
 897
 898 }
 899
 900 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortionCourse() {
 901   //
 902   // Initialization of the Lookup table which contains the solutions of the
 903   // "space charge" (poisson) problem
 904   //
 905   // The sum-up uses a look-up table which contains different discretized Space charge fields
 906   // in order to calculate the corresponding field deviations due to a given (discretized)
 907   // space charge distribution ....
 908   //
 909   // Method of calculation: Weighted sum-up of the different fields within the look up table
 910   // Note: Full 3d version: Course and slow ...
 911
 912   if (fInitLookUp) {
 913     AliInfo("Lookup table was already initialized!");
 914     //    return;
 915   }
 916
 917   AliInfo("Preparation of the weighted look-up table");
 918
 919   TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
 920   if ( !f ) {
 921     AliError("Precalculated POC-looup-table could not be found");
 922     return;
 923   }
 924
 925   // units are in [m]
 926   TVector *gridf  = (TVector*) f->Get("constants");
 927   TVector &grid = *gridf;
 928   TMatrix *coordf  = (TMatrix*) f->Get("coordinates");
 929   TMatrix &coord = *coordf;
 930   TMatrix *coordPOCf  = (TMatrix*) f->Get("POCcoord");
 931   TMatrix &coordPOC = *coordPOCf;
 932
 933   Bool_t flagRadSym = 0;
 934   if (grid(1)==1 && grid(4)==1) {
 935     //   AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
 936     flagRadSym=1;
 937   }
 938
 939   Int_t rows      = (Int_t)grid(0);   // number of points in r direction
 940   Int_t phiSlices = (Int_t)grid(1);   // number of points in phi
 941   Int_t columns   = (Int_t)grid(2);   // number of points in z direction
 942
 943   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius)/(rows-1);   // unit in [cm]
 944   const Float_t gridSizePhi =  TMath::TwoPi()/phiSlices;         // unit in [rad]
 945   const Float_t gridSizeZ   =  fgkTPCZ0/(columns-1);                  // unit in [cm]
 946
 947   // temporary matrices needed for the calculation
 948   TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
 949   TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
 950
 951   TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
 952   TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
 953
 954
 955   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 956
 957     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 958     arrayofEphiA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 959     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 960     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 961     arrayofEphiC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 962     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 963
 964     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 965     arrayofEphioverEzA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 966     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 967     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 968     arrayofEphioverEzC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 969     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 970
 971     // Set the values to zero the lookup tables
 972     // not necessary, it is done in the constructor of TMatrix - code deleted
 973
 974   }
 975
 976   // list of points as used in the interpolation (during sum up)
 977   Double_t  rlist[kNRows], zedlist[kNColumns] , philist[kNPhiSlices];
 978   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 979     philist[k] =  gridSizePhi * k;
 980     for ( Int_t i = 0 ; i < rows ; i++ )    {
 981       rlist[i] = fgkIFCRadius + i*gridSizeR ;
 982       for ( Int_t j = 0 ; j < columns ; j++ ) {
 983         zedlist[j]  = j * gridSizeZ ;
 984       }
 985     }
 986   } // only done once
 987
 988
 989   TTree *treePOC = (TTree*)f->Get("POCall");
 990
 991   TVector *bEr  = 0;   TVector *bEphi= 0;   TVector *bEz  = 0;
 992
 993   treePOC->SetBranchAddress("Er",&bEr);
 994   if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
 995   treePOC->SetBranchAddress("Ez",&bEz);
 996
 997
 998   // Read the complete tree and do a weighted sum-up over the POC configurations
 999   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1000
1001   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
1002   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
1003
1004   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // ------------- loop over POC configurations in tree
1005
1006     treePOC->GetEntry(itreepC);
1007
1008     // center of the POC voxel in [meter]
1009     Double_t r0 = coordPOC(ipC,0);
1010     Double_t phi0 = coordPOC(ipC,1);
1011     Double_t z0 = coordPOC(ipC,2);
1012
1013     ipC++; // POC configuration counter
1014
1015     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1016     // note: coordinates are in [cm]
1017     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
1018     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
1019
1020     // Summing up the vector components according to their weight
1021
1022     Int_t ip = 0;
1023     for ( Int_t j = 0 ; j < columns ; j++ ) {
1024       for ( Int_t i = 0 ; i < rows ; i++ )    {
1025         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1026
1027           // check wether the coordinates were screwed
1028           if (TMath::Abs((coord(0,ip)*100-rlist[i]))>1 ||
1029               TMath::Abs((coord(1,ip)-philist[k]))>1 ||
1030               TMath::Abs((coord(2,ip)*100-zedlist[j]))>1) {
1031             AliError("internal error: coordinate system was screwed during the sum-up");
1032             printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
1033             printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist[i],philist[k],zedlist[j]);
1034             AliError("Don't trust the results of the space charge calculation!");
1035           }
1036
1037           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1038           // This will be the most frequent usage (hopefully)
1039           // That's why we have to do this here ...
1040
1041           TMatrixD &erA   =  *arrayofErA[k]  ;
1042           TMatrixD &ephiA =  *arrayofEphiA[k];
1043           TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
1044
1045           TMatrixD &erC   =  *arrayofErC[k]  ;
1046           TMatrixD &ephiC =  *arrayofEphiC[k];
1047           TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
1048
1049           // Sum up - Efield values in [V/m] -> transition to [V/cm]
1050           erA(i,j) += ((*bEr)(ip)) * weightA /100;
1051           erC(i,j) += ((*bEr)(ip)) * weightC /100;
1052           if (!flagRadSym) {
1053             ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
1054             ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
1055           }
1056           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
1057           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
1058
1059           // increase the counter
1060           ip++;
1061         }
1062       }
1063     } // end coordinate loop
1064
1065
1066     // Rotation and summation in the rest of the dPhiSteps
1067     // which were not stored in the this tree due to storage & symmetry reasons
1068
1069     Int_t phiPoints = (Int_t) grid(1);
1070     Int_t phiPOC    = (Int_t) grid(4);
1071
1072     //   printf("%d %d\n",phiPOC,flagRadSym);
1073
1074     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
1075
1076       r0 = coordPOC(ipC,0);
1077       phi0 = coordPOC(ipC,1);
1078       z0 = coordPOC(ipC,2);
1079
1080       ipC++; // POC conf. counter
1081
1082       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
1083       // note: coordinates are in [cm]
1084       weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
1085       weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
1086
1087       //     printf("%f %f %f: %e %e\n",r0,phi0,z0,weightA,weightC);
1088
1089       // Summing up the vector components according to their weight
1090       ip = 0;
1091       for ( Int_t j = 0 ; j < columns ; j++ ) {
1092         for ( Int_t i = 0 ; i < rows ; i++ )    {
1093           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
1094
1095             // Note: rotating the coordinated during the sum up
1096
1097             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
1098             Int_t ipR = -1;
1099
1100             if ((ip%phiPoints)>=rotVal) {
1101               ipR = ip-rotVal;
1102             } else {
1103               ipR = ip+(phiPoints-rotVal);
1104             }
1105
1106             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
1107             // This will be the most frequent usage
1108             // That's why we have to do this here and not outside the loop ...
1109
1110             TMatrixD &erA   =  *arrayofErA[k]  ;
1111             TMatrixD &ephiA =  *arrayofEphiA[k];
1112             TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
1113
1114             TMatrixD &erC   =  *arrayofErC[k]  ;
1115             TMatrixD &ephiC =  *arrayofEphiC[k];
1116             TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
1117
1118             // Sum up - Efield values in [V/m] -> transition to [V/cm]
1119             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
1120             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
1121             if (!flagRadSym) {
1122               ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
1123               ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
1124             }
1125             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
1126             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
1127
1128             // increase the counter
1129             ip++;
1130           }
1131         }
1132       } // end coordinate loop
1133
1134     } // end phi-POC summation (phiiC)
1135
1136
1137     // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
1138
1139   }
1140
1141
1142
1143   // -------------------------------------------------------------------------------
1144   // Division by the Ez (drift) field and integration along z
1145
1146   AliInfo("Division and integration");
1147
1148   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
1149
1150   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
1151
1152     // matrices holding the solution - summation of POC charges // see above
1153     TMatrixD &erA   =  *arrayofErA[k]  ;
1154     TMatrixD &ephiA =  *arrayofEphiA[k];
1155     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
1156     TMatrixD &erC   =  *arrayofErC[k]  ;
1157     TMatrixD &ephiC =  *arrayofEphiC[k];
1158     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
1159
1160     // matrices which will contain the integrated fields (divided by the drift field)
1161     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
1162     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA[k];
1163     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
1164     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
1165     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC[k];
1166     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
1167
1168     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
1169       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
1170         // Count backwards to facilitate integration over Z
1171
1172         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
1173
1174         erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0;  deltaEzA(i,j) = 0;
1175         erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0;  deltaEzC(i,j) = 0;
1176
1177         for ( Int_t m = j ; m < columns ; m++ ) { // integration
1178
1179           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
1180           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
1181           if (!flagRadSym) {
1182             ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
1183             ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
1184           }
1185           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
1186           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
1187
1188           if ( index != 4 )  index = 4; else index = 2 ;
1189
1190         }
1191
1192         if ( index == 4 ) {
1193           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
1194           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
1195           if (!flagRadSym) {
1196             ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
1197             ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
1198           }
1199           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
1200           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
1201         }
1202         if ( index == 2 ) {
1203           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
1204           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
1205           if (!flagRadSym) {
1206             ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
1207             ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
1208           }
1209           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
1210           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
1211         }
1212         if ( j == columns-2 ) {
1213           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
1214           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
1215           if (!flagRadSym) {
1216             ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
1217             ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
1218           }
1219           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
1220           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
1221         }
1222         if ( j == columns-1 ) {
1223           erOverEzA(i,j)   = 0;
1224           erOverEzC(i,j)   = 0;
1225           if (!flagRadSym) {
1226             ephiOverEzA(i,j) = 0;
1227             ephiOverEzC(i,j) = 0;
1228           }
1229           deltaEzA(i,j)    = 0;
1230           deltaEzC(i,j)    = 0;
1231         }
1232       }
1233     }
1234
1235   }
1236
1237
1238
1239   AliInfo("Interpolation to Standard grid");
1240
1241   // -------------------------------------------------------------------------------
1242   // Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
1243
1244   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
1245
1246   Double_t  r, phi, z ;
1247   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1248     phi = fgkPhiList[k] ;
1249
1250     TMatrixD &erOverEz   =  *fLookUpErOverEz[k]  ;
1251     TMatrixD &ephiOverEz =  *fLookUpEphiOverEz[k];
1252     TMatrixD &deltaEz    =  *fLookUpDeltaEz[k]   ;
1253
1254     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1255
1256       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
1257
1258       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1259         r = fgkRList[i] ;
1260
1261         // Interpolate Lookup tables onto standard grid
1262         if (fgkZList[j]>0) {
1263           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1264                                                rlist, zedlist, philist, arrayofEroverEzA  );
1265           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1266                                                rlist, zedlist, philist, arrayofEphioverEzA);
1267           deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1268                                                rlist, zedlist, philist, arrayofDeltaEzA   );
1269         } else {
1270           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1271                                                rlist, zedlist, philist, arrayofEroverEzC  );
1272           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1273                                                rlist, zedlist, philist, arrayofEphioverEzC);
1274           deltaEz(i,j)  = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
1275                                                rlist, zedlist, philist, arrayofDeltaEzC   );
1276           // negative coordinate system on C side
1277         }
1278
1279       } // end r loop
1280     } // end z loop
1281   } // end phi loop
1282
1283
1284   // clear the temporary arrays lists
1285   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
1286
1287     delete arrayofErA[k];
1288     delete arrayofEphiA[k];
1289     delete arrayofdEzA[k];
1290     delete arrayofErC[k];
1291     delete arrayofEphiC[k];
1292     delete arrayofdEzC[k];
1293
1294     delete arrayofEroverEzA[k];
1295     delete arrayofEphioverEzA[k];
1296     delete arrayofDeltaEzA[k];
1297     delete arrayofEroverEzC[k];
1298     delete arrayofEphioverEzC[k];
1299     delete arrayofDeltaEzC[k];
1300
1301   }
1302
1303   fInitLookUp = kTRUE;
1304
1305 }
1306
1307
1308 void AliTPCSpaceCharge3D::SetSCDataFileName(TString fname) {
1309   //
1310   // Set & load the Space charge density distribution from a file
1311   // (linear interpolation onto a standard grid)
1312   //
1313
1314
1315   fSCDataFileName = fname;
1316
1317   TFile *f = new TFile(fSCDataFileName.Data(),"READ");
1318   if (!f) {
1319     AliError(Form("File %s, which should contain the space charge distribution, could not be found",
1320                   fSCDataFileName.Data()));
1321     return;
1322   }
1323
1324   TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
1325   if (!densityRZ) {
1326     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1327                   fSCDataFileName.Data(),"SpaceChargeInRZ"));
1328     return;
1329   }
1330
1331   TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
1332   if (!densityRPhi) {
1333     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
1334                   fSCDataFileName.Data(),"SpaceChargeInRPhi"));
1335     return;
1336   }
1337
1338
1339   Double_t  r, phi, z ;
1340
1341   TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1342   TMatrixD &scDensityInRPhiA   =  *fSCdensityInRPhiA;
1343   TMatrixD &scDensityInRPhiC   =  *fSCdensityInRPhiC;
1344   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
1345     phi = fgkPhiList[k] ;
1346     TMatrixD &scDensity   =  *fSCdensityDistribution[k]  ;
1347     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
1348       z = fgkZList[j] ;
1349       for ( Int_t i = 0 ; i < kNR ; i++ ) {
1350         r = fgkRList[i] ;
1351
1352         // partial load in (r,z)
1353         if (k==0) // do just once
1354           scDensityInRZ(i,j) =  densityRZ->Interpolate(r,z);
1355
1356         // partial load in (r,phi)
1357         if ( j==0 || j == kNZ/2 ) {
1358           if (z>0)
1359             scDensityInRPhiA(i,k) =  densityRPhi->Interpolate(r,phi,0.499);  // A side
1360           else
1361             scDensityInRPhiC(i,k) =  densityRPhi->Interpolate(r,phi,-0.499); // C side
1362         }
1363
1364         // Full 3D configuration ...
1365         if (z>0)
1366            scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
1367         else
1368            scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
1369       }
1370     }
1371   }
1372
1373   f->Close();
1374
1375   fInitLookUp = kFALSE;
1376
1377
1378 }
1379
1380
1381 Float_t  AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
1382   //
1383   // returns the (input) space charge density at a given point according
1384   // Note: input in [cm], output in [C/m^3/e0] !!
1385   //
1386
1387   while (phi<0) phi += TMath::TwoPi();
1388   while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
1389
1390
1391   // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
1392   Int_t order = 1; //
1393   Float_t sc = 0;
1394
1395   if (mode == 0) { // return full load
1396     sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
1397                             fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
1398
1399   } else if (mode == 1) { // return partial load in (r,z)
1400     TMatrixD &scDensityInRZ   =  *fSCdensityInRZ;
1401     sc = Interpolate2DTable(order, r, z, kNR, kNZ, fgkRList, fgkZList, scDensityInRZ );
1402
1403   } else if (mode == 2) { // return partial load in (r,phi)
1404
1405     if (z>0) {
1406       TMatrixD &scDensityInRPhi   =  *fSCdensityInRPhiA;
1407       sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1408     } else {
1409       TMatrixD &scDensityInRPhi   =  *fSCdensityInRPhiC;
1410       sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
1411     }
1412
1413   } else {
1414     // should i give a warning?
1415     sc = 0;
1416   }
1417
1418   //  printf("%f %f %f: %f\n",r,phi,z,sc);
1419
1420   return sc;
1421 }
1422
1423
1424 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
1425   //
1426   // return a simple histogramm containing the space charge distribution (input for the calculation)
1427   //
1428
1429   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1430                      nx,-250.,250.,ny,-250.,250.);
1431
1432   for (Int_t iy=1;iy<=ny;++iy) {
1433     Double_t yp = h->GetYaxis()->GetBinCenter(iy);
1434     for (Int_t ix=1;ix<=nx;++ix) {
1435       Double_t xp = h->GetXaxis()->GetBinCenter(ix);
1436
1437       Float_t r = TMath::Sqrt(xp*xp+yp*yp);
1438       Float_t phi = TMath::ATan2(yp,xp);
1439
1440       if (85.<=r && r<=250.) {
1441         Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1442         h->SetBinContent(ix,iy,sc);
1443       } else {
1444         h->SetBinContent(ix,iy,0.);
1445       }
1446     }
1447   }
1448
1449   return h;
1450 }
1451
1452 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
1453   //
1454   // return a simple histogramm containing the space charge distribution (input for the calculation)
1455   //
1456
1457   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
1458                      nz,-250.,250.,nr,85.,250.);
1459
1460   for (Int_t ir=1;ir<=nr;++ir) {
1461     Float_t r = h->GetYaxis()->GetBinCenter(ir);
1462     for (Int_t iz=1;iz<=nz;++iz) {
1463       Float_t z = h->GetXaxis()->GetBinCenter(iz);
1464       Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
1465       h->SetBinContent(iz,ir,sc);
1466     }
1467   }
1468
1469   return h;
1470 }
1471
1472 void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
1473   //
1474   // Example on how to write a Space charge distribution into a File
1475   //  (see below: estimate from scaling STAR measurements to Alice)
1476   // Charge distribution is splitted into two (RZ and RPHI) in order to speed up
1477   // the needed calculation time
1478   //
1479
1480   TFile *f = new TFile(fname,"RECREATE");
1481
1482   // some grid, not too course
1483   Int_t nr = 350;
1484   Int_t nphi = 180;
1485   Int_t nz = 500;
1486
1487   Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
1488   Double_t dphi = TMath::TwoPi()/(nphi+1);
1489   Double_t dz = 500./(nz+1);
1490   Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
1491
1492
1493   // Charge distribution in ZR (rotational symmetric) ------------------
1494
1495   TH2F *histoZR = new TH2F("chargeZR","chargeZR",
1496                            nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1497                            nz,-250-dz-safty,250+dz+safty);
1498
1499   for (Int_t ir=1;ir<=nr;++ir) {
1500     Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
1501     for (Int_t iz=1;iz<=nz;++iz) {
1502       Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
1503
1504       // recalculation to meter
1505       Double_t lZ = 2.5; // approx. TPC drift length
1506       Double_t rpM = rp/100.; // in [m]
1507       Double_t zpM = TMath::Abs(zp/100.); // in [m]
1508
1509       // setting of mb multiplicity and Interaction rate
1510       Double_t multiplicity = 950;
1511       Double_t intRate = 7800;
1512
1513       // calculation of "scaled" parameters
1514       Double_t a = multiplicity*intRate/79175;
1515       Double_t b = a/lZ;
1516
1517       Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
1518
1519       charge = charge*fgke0; // [C/m^3]
1520
1521       if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
1522
1523       //  charge = 0; // for tests
1524       histoZR->SetBinContent(ir,iz,charge);
1525     }
1526   }
1527
1528   histoZR->Write("SpaceChargeInRZ");
1529
1530
1531   // Charge distribution in RPhi (e.g. Floating GG wire) ------------
1532
1533   TH3F *histoRPhi = new TH3F("chargeRPhi","chargeRPhi",
1534                              nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
1535                              nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
1536                              2,-1,1); // z part - to allow A and C side differences
1537
1538   // some 'arbitrary' GG leaks
1539   Int_t   nGGleaks = 5;
1540   Double_t secPosA[5]    = {3,6,6,11,13};         // sector
1541   Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm
1542   Double_t secPosC[5]    = {1,8,12,15,15};        // sector
1543   Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm
1544
1545   for (Int_t ir=1;ir<=nr;++ir) {
1546     Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir);
1547     for (Int_t iphi=1;iphi<=nphi;++iphi) {
1548       Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi);
1549       for (Int_t iz=1;iz<=2;++iz) {
1550         Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz);
1551
1552         Double_t charge = 0;
1553
1554         for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks
1555
1556           // A side
1557           Double_t secPos = secPosA[igg];
1558           Double_t radialPos = radialPosA[igg];
1559
1560           if (zp<0) { // C side
1561             secPos = secPosC[igg];
1562             radialPos = radialPosC[igg];
1563           }
1564
1565           // some 'arbitrary' GG leaks
1566           if (  (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice
1567             if ( rp>(radialPos-2.5) && rp<(radialPos+2.5))  // 5 cm slice
1568               charge = 300;
1569           }
1570
1571         }
1572
1573         charge = charge*fgke0; // [C/m^3]
1574
1575         histoRPhi->SetBinContent(ir,iphi,iz,charge);
1576       }
1577     }
1578   }
1579
1580   histoRPhi->Write("SpaceChargeInRPhi");
1581
1582   f->Close();
1583
1584 }
1585
1586
1587 void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
1588   //
1589   // Print function to check the settings of the boundary vectors
1590   // option=="a" prints the C0 and C1 coefficents for calibration purposes
1591   //
1592
1593   TString opt = option; opt.ToLower();
1594   printf("%s\n",GetTitle());
1595   printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
1596   printf("   SC correction factor: %f \n",fCorrectionFactor);
1597
1598   if (opt.Contains("a")) { // Print all details
1599     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
1600     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
1601   }
1602
1603   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
1604
1605 }