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 // AliTPCSpaceCharge3D class                                                //
  19 // The class calculates the space point distortions due to a space charge   //
  20 // effect ....                                                              //
  21 // Method of calculation:                                                   //
  22 // The analytical solution for the poisson problem in 3D (cylindrical coord)//
  23 // is used in form of look up tables. PieceOfCake (POC) voxel were pre-     //
  24 // calculated and can be sumed up (weighted) according to the what is needed//
  25 //                                                                          //
  26 // The class allows "effective Omega Tau" corrections.                      //
  27 //                                                                          //
  28 // NOTE: This class is not  capable of calculation z distortions due to     //
  29 //       drift velocity change in dependence of the electric field!!!       //
  30 //                                                                          //
  31 // date: 19/06/2010                                                         //
  32 // Authors: Stefan Rossegger                                                //
  33 //                                                                          //
  34 // Example usage:                                                           //
  35 //////////////////////////////////////////////////////////////////////////////
  36
  37 #include "AliMagF.h"
  38 #include "TGeoGlobalMagField.h"
  39 #include "AliTPCcalibDB.h"
  40 #include "AliTPCParam.h"
  41 #include "AliLog.h"
  42 #include "TH2F.h"
  43 #include "TH3F.h"
  44 #include "TFile.h"
  45 #include "TVector.h"
  46 #include "TMatrix.h"
  47 #include "TMatrixD.h"
  48
  49 #include "TMath.h"
  50 #include "AliTPCROC.h"
  51 #include "AliTPCSpaceCharge3D.h"
  52
  53 ClassImp(AliTPCSpaceCharge3D)
  54
  55 AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
  56   : AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
  57     fC0(0.),fC1(0.),
  58     fCorrectionFactor(1.),
  59     fInitLookUp(kFALSE),
  60     fSCDataFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root"),
  61     fSCLookUpPOCsFileName((char*)"$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root")
  62 {
  63   //
  64   // default constructor
  65   //
  66
  67   // Array which will contain the solution according to the setted charge density distribution
  68   // see InitSpaceCharge3DDistortion() function
  69   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
  70     fLookUpErOverEz[k]   =  new TMatrixD(kNR,kNZ);
  71     fLookUpEphiOverEz[k] =  new TMatrixD(kNR,kNZ);
  72     fLookUpDeltaEz[k]    =  new TMatrixD(kNR,kNZ);
  73     fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
  74   }
  75
  76   printf("%s\n",fSCDataFileName);
  77   printf("%s\n",fSCLookUpPOCsFileName);
  78   SetSCDataFileName(fSCDataFileName);
  79
  80 }
  81
  82 AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
  83   //
  84   // default destructor
  85   //
  86
  87   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
  88     delete fLookUpErOverEz[k];
  89     delete fLookUpEphiOverEz[k];
  90     delete fLookUpDeltaEz[k];
  91     delete fSCdensityDistribution[k];
  92   }
  93
  94   //  delete fSCdensityDistribution;
  95 }
  96
  97
  98
  99 void AliTPCSpaceCharge3D::Init() {
 100   //
 101   // Initialization funtion
 102   //
 103
 104   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 105   if (!magF) AliError("Magneticd field - not initialized");
 106   Double_t bzField = magF->SolenoidField()/10.; //field in T
 107   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 108   if (!param) AliError("Parameters - not initialized");
 109   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 110   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 111   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 112   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 113   SetOmegaTauT1T2(wt,fT1,fT2);
 114
 115   InitSpaceCharge3DDistortion(); // fill the look up table
 116 }
 117
 118 void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
 119   //
 120   // Update function
 121   //
 122   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 123   if (!magF) AliError("Magneticd field - not initialized");
 124   Double_t bzField = magF->SolenoidField()/10.; //field in T
 125   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 126   if (!param) AliError("Parameters - not initialized");
 127   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
 128   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 129   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 130   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 131   SetOmegaTauT1T2(wt,fT1,fT2);
 132
 133   //  SetCorrectionFactor(1.); // should come from some database
 134
 135 }
 136
 137
 138
 139 void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 140   //
 141   // Calculates the correction due the Space Charge effect within the TPC drift volume
 142   //
 143
 144   if (!fInitLookUp) {
 145     AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
 146     InitSpaceCharge3DDistortion();
 147   }
 148
 149   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 150
 151   Double_t intEr, intEphi, intdEz ;
 152   Double_t r, phi, z ;
 153   Int_t    sign;
 154
 155   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 156   phi    =  TMath::ATan2(x[1],x[0]) ;
 157   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 158   z      =  x[2] ;                                         // Create temporary copy of x[2]
 159
 160   if ( (roc%36) < 18 ) {
 161     sign =  1;       // (TPC A side)
 162   } else {
 163     sign = -1;       // (TPC C side)
 164   }
 165
 166   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 167   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 168
 169
 170   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 171     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 172
 173   // Get the Er and Ephi field integrals plus the integral over DeltaEz
 174   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 175                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
 176   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 177                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
 178   intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 179                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );
 180
 181   // Calculate distorted position
 182   if ( r > 0.0 ) {
 183     phi =  phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
 184     r   =  r   + fCorrectionFactor *( fC0*intEr   + fC1*intEphi );
 185   }
 186   Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
 187
 188   // Calculate correction in cartesian coordinates
 189   dx[0] = r * TMath::Cos(phi) - x[0];
 190   dx[1] = r * TMath::Sin(phi) - x[1];
 191   dx[2] = dz;  // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
 192
 193 }
 194
 195
 196 void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
 197   //
 198   // Initialization of the Lookup table which contains the solutions of the
 199   // "space charge" (poisson) problem
 200   //
 201   // The sum-up uses a look-up table which contains different discretized Space charge fields
 202   // in order to calculate the corresponding field deviations due to a given (discretized)
 203   // space charge distribution ....
 204   //
 205   // Method of calculation: Weighted sum up of the different fields within the look up table
 206   //
 207
 208   if (fInitLookUp) {
 209     AliInfo("Lookup table was already initialized!");
 210     //    return;
 211   }
 212
 213   AliInfo("Preparation of the weighted look-up table");
 214
 215   TFile *f = new TFile(fSCLookUpPOCsFileName);
 216   if (!f) {
 217     AliError("Precalculated POC-looup-table could not be found");
 218     return;
 219   }
 220
 221   // units are in [m]
 222   TVector *gridf  = (TVector*) f->Get("constants");
 223   TVector &grid = *gridf;
 224   TMatrix *coordf  = (TMatrix*) f->Get("coordinates");
 225   TMatrix &coord = *coordf;
 226   TMatrix *coordPOCf  = (TMatrix*) f->Get("POCcoord");
 227   TMatrix &coordPOC = *coordPOCf;
 228
 229   Bool_t flagRadSym = 0;
 230   if (grid(1)==1 && grid(4)==1) {
 231     AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
 232     flagRadSym=1;
 233   }
 234
 235   Int_t rows      = (Int_t)grid(0);   // number of points in r direction
 236   Int_t phiSlices = (Int_t)grid(1);   // number of points in phi
 237   Int_t columns   = (Int_t)grid(2);   // number of points in z direction
 238
 239   const Float_t gridSizeR   =  grid(6)*100;    // unit in [cm]
 240   const Float_t gridSizePhi   =  grid(7);  // unit in [rad]
 241   const Float_t gridSizeZ =  grid(8)*100;      // unit in [cm]
 242
 243   // temporary matrices needed for the calculation
 244   TMatrixD *arrayofErA[phiSlices], *arrayofEphiA[phiSlices], *arrayofdEzA[phiSlices];
 245   TMatrixD *arrayofErC[phiSlices], *arrayofEphiC[phiSlices], *arrayofdEzC[phiSlices];
 246
 247   TMatrixD *arrayofEroverEzA[phiSlices], *arrayofEphioverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
 248   TMatrixD *arrayofEroverEzC[phiSlices], *arrayofEphioverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
 249
 250
 251   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 252
 253     arrayofErA[k]   =   new TMatrixD(rows,columns) ;
 254     arrayofEphiA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 255     arrayofdEzA[k]  =   new TMatrixD(rows,columns) ;
 256     arrayofErC[k]   =   new TMatrixD(rows,columns) ;
 257     arrayofEphiC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 258     arrayofdEzC[k]  =   new TMatrixD(rows,columns) ;
 259
 260     arrayofEroverEzA[k]   =   new TMatrixD(rows,columns) ;
 261     arrayofEphioverEzA[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 262     arrayofDeltaEzA[k]    =   new TMatrixD(rows,columns) ;
 263     arrayofEroverEzC[k]   =   new TMatrixD(rows,columns) ;
 264     arrayofEphioverEzC[k] =   new TMatrixD(rows,columns) ; // zero if radial symmetric
 265     arrayofDeltaEzC[k]    =   new TMatrixD(rows,columns) ;
 266
 267     // Set the values to zero the lookup tables
 268     // not necessary, it is done in the constructor of TMatrix - code deleted
 269
 270   }
 271
 272   // list of points as used in the interpolation (during sum up)
 273   Double_t  rlist[rows], zedlist[columns] , philist[phiSlices];
 274   for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 275     philist[k] =  gridSizePhi * k;
 276     for ( Int_t i = 0 ; i < rows ; i++ )    {
 277       rlist[i] = fgkIFCRadius + i*gridSizeR ;
 278       for ( Int_t j = 0 ; j < columns ; j++ ) {
 279         zedlist[j]  = j * gridSizeZ ;
 280       }
 281     }
 282   } // only done once
 283
 284
 285   // Prepare summation Vectors
 286   //  Int_t numPosInd = (Int_t)(grid(0)*grid(1)*grid(2)); // Number of fulcrums in the look-up table
 287   //  Int_t numPOC = (Int_t)(grid(3)*grid(4)*grid(5));    // Number of POC conf. in the look-up table
 288
 289   TTree *treePOC = (TTree*)f->Get("POCall");
 290
 291   TVector *bEr  = 0;   TVector *bEphi= 0;   TVector *bEz  = 0;
 292
 293   treePOC->SetBranchAddress("Er",&bEr);
 294   if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
 295   treePOC->SetBranchAddress("Ez",&bEz);
 296
 297   /*  TBranch *b2Er = treePOC->GetBranch("Er");
 298       TBranch *b2Ephi =0;
 299       if (!flagRadSym) b2Ephi = treePOC->GetBranch("Ephi");
 300       TBranch *b2Ez = treePOC->GetBranch("Ez");
 301   */
 302
 303   // Read the complete tree and do a weighted sum-up over the POC configurations
 304   // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 305
 306   Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
 307   Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
 308
 309   for (Int_t itreepC=0; itreepC<treeNumPOC; itreepC++) { // loop over POC configurations in tree
 310
 311     // if (!flagRadSym) cout<<ipC<<endl;
 312
 313     treePOC->GetEntry(itreepC);
 314
 315     /*    b2Er->GetEntry(itreepC);
 316           if (!flagRadSym) b2Ephi->GetEntry(itreepC);
 317           b2Ez->GetEntry(itreepC);
 318     */
 319
 320     // center of the POC voxel in [meter]
 321     Double_t r0 = coordPOC(ipC,0);
 322     Double_t phi0 = coordPOC(ipC,1);
 323     Double_t z0 = coordPOC(ipC,2);
 324
 325     ipC++; // POC conf. counter
 326
 327     // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 328     // note: coordinates are in [cm]
 329     Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
 330     Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
 331
 332     // Summing up the vector components according to their weight
 333
 334     Int_t ip = 0;
 335     for ( Int_t j = 0 ; j < columns ; j++ ) {
 336       for ( Int_t i = 0 ; i < rows ; i++ )    {
 337         for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 338
 339           // check wether the coordinates were screwed
 340           if ((coord(0,ip)*100-rlist[i])>0.01 ||
 341               (coord(1,ip)-philist[k])>0.01 ||
 342               (coord(2,ip)*100-zedlist[j])>0.01) {
 343             AliError("internal error: coordinate system was screwed during the sum-up");
 344             printf("lookup: (r,phi,z)=(%lf,%lf,%lf)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
 345             printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist[i],philist[k],zedlist[j]);
 346             AliError("Don't trust the results of the space charge calibration!");
 347           }
 348
 349           // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 350           // This will be the most frequent usage (hopefully)
 351           // That's why we have to do this here ...
 352
 353           TMatrixD &erA   =  *arrayofErA[k]  ;
 354           TMatrixD &ephiA =  *arrayofEphiA[k];
 355           TMatrixD &dEzA    =  *arrayofdEzA[k]   ;
 356
 357           TMatrixD &erC   =  *arrayofErC[k]  ;
 358           TMatrixD &ephiC =  *arrayofEphiC[k];
 359           TMatrixD &dEzC    =  *arrayofdEzC[k]   ;
 360
 361           // Sum up - Efield values in [V/m] -> transition to [V/cm]
 362           erA(i,j) += ((*bEr)(ip)) * weightA /100;
 363           erC(i,j) += ((*bEr)(ip)) * weightC /100;
 364           if (!flagRadSym) {
 365             ephiA(i,j) += ((*bEphi)(ip)) * weightA; // [V/rad]
 366             ephiC(i,j) += ((*bEphi)(ip)) * weightC; // [V/rad]
 367           }
 368           dEzA(i,j)  += ((*bEz)(ip)) * weightA /100;
 369           dEzC(i,j)  += ((*bEz)(ip)) * weightC /100;
 370
 371           // increase the counter
 372           ip++;
 373         }
 374       }
 375     }
 376
 377     // Rotation and summation in the rest of the dPhiSteps
 378     // which were not stored in the tree due to storage & symmetry reasons
 379
 380     Int_t phiPoints = (Int_t) grid(1);
 381     Int_t phiPOC    = (Int_t) grid(4);
 382
 383     // printf("%d %d\n",phiPOC,flagRadSym);
 384
 385     for (Int_t phiiC = 1; phiiC<phiPOC; phiiC++) { // just used for non-radial symetric table
 386
 387       r0 = coordPOC(ipC,0);
 388       phi0 = coordPOC(ipC,1);
 389       z0 = coordPOC(ipC,2);
 390
 391       ipC++; // POC conf. counter
 392
 393       // weights (charge density) at POC position on the A and C side (in C/m^3/e0)
 394       // note: coordinates are in [cm]
 395       weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
 396       weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
 397
 398
 399       // Summing up the vector components according to their weight
 400       ip = 0;
 401       for ( Int_t j = 0 ; j < columns ; j++ ) {
 402         for ( Int_t i = 0 ; i < rows ; i++ )    {
 403           for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
 404
 405             // Note: rotating the coordinated during the sum up
 406
 407             Int_t rotVal = (phiPoints/phiPOC)*phiiC;
 408             Int_t ipR = -1;
 409
 410             if ((ip%phiPoints)>=rotVal) {
 411               ipR = ip-rotVal;
 412             } else {
 413               ipR = ip+(phiPoints-rotVal);
 414             }
 415
 416             // unfortunately, the lookup tables were produced to be faster for phi symmetric charges
 417             // This will be the most frequent usage (hopefully)
 418             // That's why we have to do this here and not outside the loop ...
 419
 420             TMatrixD &erA   =  *arrayofErA[k]  ;
 421             TMatrixD &ephiA =  *arrayofEphiA[k];
 422             TMatrixD &dEzA  =  *arrayofdEzA[k]   ;
 423
 424             TMatrixD &erC   =  *arrayofErC[k]  ;
 425             TMatrixD &ephiC =  *arrayofEphiC[k];
 426             TMatrixD &dEzC  =  *arrayofdEzC[k]   ;
 427
 428             // Sum up - Efield values in [V/m] -> transition to [V/cm]
 429             erA(i,j) += ((*bEr)(ipR)) * weightA /100;
 430             erC(i,j) += ((*bEr)(ipR)) * weightC /100;
 431             if (!flagRadSym) {
 432               ephiA(i,j) += ((*bEphi)(ipR)) * weightA; // [V/rad]
 433               ephiC(i,j) += ((*bEphi)(ipR)) * weightC; // [V/rad]
 434             }
 435             dEzA(i,j)  += ((*bEz)(ipR)) * weightA /100;
 436             dEzC(i,j)  += ((*bEz)(ipR)) * weightC /100;
 437
 438             // increase the counter
 439             ip++;
 440           }
 441         }
 442       } // end coordinate loop
 443
 444     } // end phi-POC summation (phiiC)
 445
 446
 447     // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
 448
 449   } // end POC loop
 450
 451   AliInfo("Division and integration");
 452
 453   // -------------------------------------------------------------------------------
 454   // Division by the Ez (drift) field and integration along z
 455
 456   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
 457
 458   for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
 459
 460     // matrices holding the solution - summation of POC charges
 461     TMatrixD &erA   =  *arrayofErA[k]  ;
 462     TMatrixD &ephiA =  *arrayofEphiA[k];
 463     TMatrixD &dezA  =  *arrayofdEzA[k]   ;
 464     TMatrixD &erC   =  *arrayofErC[k]  ;
 465     TMatrixD &ephiC =  *arrayofEphiC[k];
 466     TMatrixD &dezC  =  *arrayofdEzC[k]   ;
 467
 468     // matrices which contain the integrated fields (divided by the drift field)
 469     TMatrixD &erOverEzA   =  *arrayofEroverEzA[k]  ;
 470     TMatrixD &ephiOverEzA =  *arrayofEphioverEzA[k];
 471     TMatrixD &deltaEzA    =  *arrayofDeltaEzA[k];
 472     TMatrixD &erOverEzC   =  *arrayofEroverEzC[k]  ;
 473     TMatrixD &ephiOverEzC =  *arrayofEphioverEzC[k];
 474     TMatrixD &deltaEzC    =  *arrayofDeltaEzC[k];
 475
 476     for ( Int_t i = 0 ; i < rows ; i++ )    { // r loop
 477       for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
 478         // Count backwards to facilitate integration over Z
 479
 480         Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
 481
 482         erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0;  deltaEzA(i,j) = 0;
 483         erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0;  deltaEzC(i,j) = 0;
 484
 485         for ( Int_t m = j ; m < columns ; m++ ) { // integration
 486
 487           erOverEzA(i,j)   += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
 488           erOverEzC(i,j)   += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField)  ;
 489           if (!flagRadSym) {
 490             ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField)  ;
 491             ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField)  ;
 492           }
 493           deltaEzA(i,j)    += index*(gridSizeZ/3.0)*dezA(i,m) ;
 494           deltaEzC(i,j)    += index*(gridSizeZ/3.0)*dezC(i,m) ;
 495
 496           if ( index != 4 )  index = 4; else index = 2 ;
 497
 498         }
 499
 500         if ( index == 4 ) {
 501           erOverEzA(i,j)   -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
 502           erOverEzC(i,j)   -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
 503           if (!flagRadSym) {
 504             ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
 505             ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
 506           }
 507           deltaEzA(i,j)    -= (gridSizeZ/3.0)*dezA(i,columns-1) ;
 508           deltaEzC(i,j)    -= (gridSizeZ/3.0)*dezC(i,columns-1) ;
 509         }
 510         if ( index == 2 ) {
 511           erOverEzA(i,j)   += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
 512           erOverEzC(i,j)   += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
 513           if (!flagRadSym) {
 514             ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
 515             ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
 516           }
 517           deltaEzA(i,j)    += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1)) ;
 518           deltaEzC(i,j)    += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1)) ;
 519         }
 520         if ( j == columns-2 ) {
 521           erOverEzA(i,j)   = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
 522           erOverEzC(i,j)   = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
 523           if (!flagRadSym) {
 524             ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
 525             ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
 526           }
 527           deltaEzA(i,j)    = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1)) ;
 528           deltaEzC(i,j)    = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1)) ;
 529         }
 530         if ( j == columns-1 ) {
 531           erOverEzA(i,j)   = 0;    erOverEzC(i,j)   = 0;
 532           if (!flagRadSym) {
 533             ephiOverEzA(i,j) = 0;  ephiOverEzC(i,j) = 0;
 534           }
 535           deltaEzA(i,j)    = 0;    deltaEzC(i,j)    = 0;
 536         }
 537       }
 538     }
 539
 540   }
 541
 542   AliInfo("Interpolation to Standard grid");
 543
 544   // -------------------------------------------------------------------------------
 545   // Interpolate results onto the standard grid which is used for all AliTPCCorrections
 546
 547   const Int_t order  = 1  ;  // Linear interpolation = 1, Quadratic = 2
 548
 549   Double_t  r, phi, z ;
 550   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 551     phi = fgkPhiList[k] ;
 552
 553     TMatrixD &erOverEz   =  *fLookUpErOverEz[k]  ;
 554     TMatrixD &ephiOverEz =  *fLookUpEphiOverEz[k];
 555     TMatrixD &deltaEz    =  *fLookUpDeltaEz[k]   ;
 556
 557     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 558
 559       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 560
 561       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 562         r = fgkRList[i] ;
 563
 564         // Interpolate Lookup tables onto standard grid
 565         if (fgkZList[j]>0) {
 566           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 567                                                rlist, zedlist, philist, arrayofEroverEzA  );
 568           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 569                                                rlist, zedlist, philist, arrayofEphioverEzA);
 570           deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 571                                                rlist, zedlist, philist, arrayofDeltaEzA   );
 572         } else {
 573           erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 574                                                rlist, zedlist, philist, arrayofEroverEzC  );
 575           ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 576                                                rlist, zedlist, philist, arrayofEphioverEzC);
 577           deltaEz(i,j)  = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
 578                                                rlist, zedlist, philist, arrayofDeltaEzC   );
 579           // negative coordinate system on C side
 580         }
 581
 582       } // end r loop
 583     } // end z loop
 584   } // end phi loop
 585
 586
 587   // clear the temporary arrays lists
 588   for ( Int_t k = 0 ; k < phiSlices ; k++ )  {
 589
 590     delete arrayofErA[k];
 591     delete arrayofEphiA[k];
 592     delete arrayofdEzA[k];
 593     delete arrayofErC[k];
 594     delete arrayofEphiC[k];
 595     delete arrayofdEzC[k];
 596
 597     delete arrayofEroverEzA[k];
 598     delete arrayofEphioverEzA[k];
 599     delete arrayofDeltaEzA[k];
 600     delete arrayofEroverEzC[k];
 601     delete arrayofEphioverEzC[k];
 602     delete arrayofDeltaEzC[k];
 603
 604   }
 605
 606
 607   fInitLookUp = kTRUE;
 608
 609 }
 610
 611
 612 void AliTPCSpaceCharge3D::SetSCDataFileName(char *const fname) {
 613   //
 614   // Set & load the Space charge density distribution from a file
 615   // (linear interpolation onto a standard grid)
 616   //
 617
 618   fSCDataFileName = fname;
 619
 620   TFile *f = new TFile(fSCDataFileName,"READ");
 621   if (!f) {
 622     AliError(Form("File %s, which should contain the space charge distribution, could not be found",
 623                   fSCDataFileName));
 624     return;
 625   }
 626
 627   TH3F *density = (TH3F*) f->Get("SpaceChargeDistribution");
 628   if (!density) {
 629     AliError(Form("The indicated file (%s) does not contain a histogram called %s",
 630                   fSCDataFileName,"SpaceChargeDistribution"));
 631     return;
 632   }
 633
 634
 635   Double_t  r, phi, z ;
 636   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 637     phi = fgkPhiList[k] ;
 638     TMatrixD &scdensity   =  *fSCdensityDistribution[k]  ;
 639     for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 640       z = TMath::Abs(fgkZList[j]) ;  // z position is symmetric
 641       for ( Int_t i = 0 ; i < kNR ; i++ ) {
 642         r = fgkRList[i] ;
 643         scdensity(i,j) = density->Interpolate(r,phi,z); // quite slow
 644         //      printf("%lf %lf %lf: %lf\n",r,phi,z,scdensity(i,j));
 645       }
 646     }
 647   }
 648
 649
 650   f->Close();
 651
 652   fInitLookUp = kFALSE;
 653
 654
 655 }
 656
 657
 658 Float_t  AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z) {
 659   //
 660   // returns the (input) space charge density at a given point according
 661   // Note: input in [cm], output in [C/m^3/e0] !!
 662   //
 663
 664   if (!fSCdensityDistribution) {
 665     printf("Scheiss is nuul - argg\n");
 666     return 0.;
 667   }
 668
 669   while (phi<0) phi += TMath::TwoPi();
 670   while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
 671
 672
 673   // Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
 674   Int_t order = 1; //
 675   Float_t sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 676                                   fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
 677
 678   //  printf("%lf %lf %lf: %lf\n",r,phi,z,sc);
 679
 680   return sc;
 681 }
 682
 683
 684 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny) {
 685   //
 686   // return a simple histogramm containing the space charge distribution (input for the calculation)
 687   //
 688
 689   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
 690                      nx,-250.,250.,ny,-250.,250.);
 691
 692   for (Int_t iy=1;iy<=ny;++iy) {
 693     Double_t yp = h->GetYaxis()->GetBinCenter(iy);
 694     for (Int_t ix=1;ix<=nx;++ix) {
 695       Double_t xp = h->GetXaxis()->GetBinCenter(ix);
 696
 697       Float_t r = TMath::Sqrt(xp*xp+yp*yp);
 698       Float_t phi = TMath::ATan2(yp,xp);
 699
 700       if (85.<=r && r<=250.) {
 701         Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
 702         h->SetBinContent(ix,iy,sc);
 703       } else {
 704         h->SetBinContent(ix,iy,0.);
 705       }
 706     }
 707   }
 708
 709   return h;
 710 }
 711
 712 TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr) {
 713   //
 714   // return a simple histogramm containing the space charge distribution (input for the calculation)
 715   //
 716
 717   TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
 718                      nz,-250.,250.,nr,85.,250.);
 719
 720   for (Int_t ir=1;ir<=nr;++ir) {
 721     Float_t r = h->GetYaxis()->GetBinCenter(ir);
 722     for (Int_t iz=1;iz<=nz;++iz) {
 723       Float_t z = h->GetXaxis()->GetBinCenter(iz);
 724       Float_t sc = GetSpaceChargeDensity(r,phi,z)/fgke0; // in [C/m^3/e0]
 725       h->SetBinContent(iz,ir,sc);
 726     }
 727   }
 728
 729   return h;
 730 }
 731
 732 void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
 733   //
 734   // Example on how to write a Space charge distribution into a File
 735   //  (see below: estimate from scaling STAR measurements to Alice)
 736   //
 737
 738   TFile *f = new TFile(fname,"RECREATE");
 739
 740   // some grid, not too course
 741   Int_t nr = 72;
 742   Int_t nphi = 180;
 743   Int_t nz = 250;
 744
 745   Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
 746   Double_t dphi = TMath::TwoPi()/(nphi+1);
 747   Double_t dz = 500./(nz+1);
 748   Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
 749
 750   TH3F *histo = new TH3F("charge","charge",
 751                          nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
 752                          nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
 753                          nz,-250-dz-safty,250+dz+safty);
 754
 755   for (Int_t ir=1;ir<=nr;++ir) {
 756     Double_t rp = histo->GetXaxis()->GetBinCenter(ir);
 757     for (Int_t iphi=1;iphi<=nphi;++iphi) {
 758       Double_t phip = histo->GetYaxis()->GetBinCenter(iphi);
 759       for (Int_t iz=1;iz<=nz;++iz) {
 760         Double_t zp = histo->GetZaxis()->GetBinCenter(iz);
 761
 762
 763         zp = TMath::Abs(zp); // symmetric on both sides of the TPC
 764
 765         // recalculation to meter
 766         Double_t lZ = 2.5; // approx. TPC drift length
 767         Double_t rpM = rp/100.; // in [m]
 768         Double_t zpM = zp/100.; // in [m]
 769
 770         // setting of mb multiplicity and Interaction rate
 771         Double_t multiplicity = 950;
 772         Double_t intRate = 8000;
 773
 774         // calculation of "scaled" parameters
 775         Double_t a = multiplicity*intRate/79175;
 776         Double_t b = a/lZ;
 777
 778         Double_t charge = (a - b * zpM)/(rpM*rpM); // charge in [C/m^3/e0]
 779
 780         // some 'arbitrary' GG leaks
 781         if (  (phip<(TMath::Pi()/9*3) && phip>(TMath::Pi()/9*2) ) ) { //||
 782           //          (phip<(TMath::Pi()/9*7) && phip>(TMath::Pi()/9*6) ) ||
 783           //          (phip<(TMath::Pi()/9*16) && phip>(TMath::Pi()/9*13) ) ){
 784           if (rp>130&&rp<135)
 785             charge = 300;
 786         }
 787
 788         //      printf("%lf %lf %lf: %lf\n",rp,phip,zp,charge);
 789
 790         charge = charge*fgke0; // [C/m^3]
 791
 792         histo->SetBinContent(ir,iphi,iz,charge);
 793       }
 794     }
 795   }
 796
 797
 798   histo->Write("SpaceChargeDistribution");
 799   f->Close();
 800
 801 }
 802
 803
 804 void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
 805   //
 806   // Print function to check the settings of the boundary vectors
 807   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 808   //
 809
 810   TString opt = option; opt.ToLower();
 811   printf("%s\n",GetTitle());
 812   printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
 813   printf("   SC correction factor: %f \n",fCorrectionFactor);
 814
 815   if (opt.Contains("a")) { // Print all details
 816     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 817     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
 818   }
 819
 820   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
 821
 822 }