/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
* Contributors are mentioned in the code where appropriate. *
* *
* Permission to use, copy, modify and distribute this software and its *
* documentation strictly for non-commercial purposes is hereby granted *
* without fee, provided that the above copyright notice appears in all *
* copies and that both the copyright notice and this permission notice *
* appear in the supporting documentation. The authors make no claims *
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
// _________________________________________________________________
//
// Begin_Html
//
AliTPCSpaceCharge3D class
// The class calculates the space point distortions due to an arbitrary space
// charge distribution in 3D.
//
// The method of calculation is based on the analytical solution for the Poisson
// problem in 3D (cylindrical coordinates). The solution is used in form of
// look up tables, where the pre calculated solutions for different voxel
// positions are stored. These voxel solutions can be summed up according
// to the weight of the position of the applied space charge distribution.
// Further details can be found in \cite[chap.5]{PhD-thesis_S.Rossegger}.
//
// The class also allows a simple scaling of the resulting distortions
// via the function SetCorrectionFactor. This speeds up the calculation
// time under the assumption, that the distortions scales linearly with the
// magnitude of the space charge distribution $\rho(r,z)$ and the shape stays
// the same at higher luminosities.
//
// In contrast to the implementation in 2D (see the class AliTPCSpaceChargeabove),
// the input charge distribution can be of arbitrary character. An example on how
// to produce a corresponding charge distribution can be found in the function
// WriteChargeDistributionToFile. In there, a $\rho(r,z) = (A-B\,z)/r^2$,
// with slightly different magnitude on the A and C side (due to the muon absorber),
// is superpositioned with a few leaking wires at arbitrary positions.
// End_Html
//
// Begin_Macro(source)
// {
// gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
// TCanvas *c2 = new TCanvas("c2","c2",500,400);
// AliTPCSpaceCharge3D sc;
// sc.WriteChargeDistributionToFile("SC_zr2_GGleaks.root");
// sc.SetSCDataFileName("SC_zr2_GGleaks.root");
// sc.SetOmegaTauT1T2(0,1,1); // B=0
// sc.InitSpaceCharge3DDistortion();
// sc.CreateHistoDRinXY(15,300,300)->Draw("colz");
// return c2;
// }
// End_Macro
//
// Begin_Html
//
// Date: 19/06/2010
// Authors: Stefan Rossegger
// End_Html
// _________________________________________________________________
#include "AliMagF.h"
#include "TGeoGlobalMagField.h"
#include "AliTPCcalibDB.h"
#include "AliTPCParam.h"
#include "AliLog.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TFile.h"
#include "TVector.h"
#include "TMatrix.h"
#include "TMatrixD.h"
#include "TMath.h"
#include "AliTPCROC.h"
#include "AliTPCSpaceCharge3D.h"
ClassImp(AliTPCSpaceCharge3D)
AliTPCSpaceCharge3D::AliTPCSpaceCharge3D()
: AliTPCCorrection("SpaceCharge3D","Space Charge - 3D"),
fC0(0.),fC1(0.),
fCorrectionFactor(1.),
fInitLookUp(kFALSE),
fSCDataFileName(""),
fSCLookUpPOCsFileName3D(""),
fSCLookUpPOCsFileNameRZ(""),
fSCLookUpPOCsFileNameRPhi(""),
fSCdensityInRZ(0),
fSCdensityInRPhiA(0),
fSCdensityInRPhiC(0)
{
//
// default constructor
//
// Array which will contain the solution according to the setted charge density distribution
// see InitSpaceCharge3DDistortion() function
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
fLookUpErOverEz[k] = new TMatrixD(kNR,kNZ);
fLookUpEphiOverEz[k] = new TMatrixD(kNR,kNZ);
fLookUpDeltaEz[k] = new TMatrixD(kNR,kNZ);
fSCdensityDistribution[k] = new TMatrixD(kNR,kNZ);
}
fSCdensityInRZ = new TMatrixD(kNR,kNZ);
fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
fSCdensityInRPhiC = new TMatrixD(kNR,kNPhi);
// location of the precalculated look up tables
fSCLookUpPOCsFileName3D="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_raw_18-18-26_17p-18p-25p-MN30.root"; // rough estimate
fSCLookUpPOCsFileNameRZ="$(ALICE_ROOT)/TPC/Calib/maps/sc_radSym_35-01-51_34p-01p-50p_MN60.root";
fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-144-26_34p-18p-01p-MN30.root";
// fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-36-26_34p-18p-01p-MN40.root";
// standard location of the space charge distibution ... can be changes
fSCDataFileName="$(ALICE_ROOT)/TPC/Calib/maps/sc_3D_distribution_Sim.root";
// SetSCDataFileName(fSCDataFileName.Data()); // should be done by the user
}
AliTPCSpaceCharge3D::~AliTPCSpaceCharge3D() {
//
// default destructor
//
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
delete fLookUpErOverEz[k];
delete fLookUpEphiOverEz[k];
delete fLookUpDeltaEz[k];
delete fSCdensityDistribution[k];
}
delete fSCdensityInRZ;
delete fSCdensityInRPhiA;
delete fSCdensityInRPhiC;
}
void AliTPCSpaceCharge3D::Init() {
//
// Initialization funtion
//
AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
if (!magF) AliError("Magneticd field - not initialized");
Double_t bzField = magF->SolenoidField()/10.; //field in T
AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
if (!param) AliError("Parameters - not initialized");
Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
// Correction Terms for effective omegaTau; obtained by a laser calibration run
SetOmegaTauT1T2(wt,fT1,fT2);
InitSpaceCharge3DDistortion(); // fill the look up table
}
void AliTPCSpaceCharge3D::Update(const TTimeStamp &/*timeStamp*/) {
//
// Update function
//
AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
if (!magF) AliError("Magneticd field - not initialized");
Double_t bzField = magF->SolenoidField()/10.; //field in T
AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
if (!param) AliError("Parameters - not initialized");
Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
// Correction Terms for effective omegaTau; obtained by a laser calibration run
SetOmegaTauT1T2(wt,fT1,fT2);
// SetCorrectionFactor(1.); // should come from some database
}
void AliTPCSpaceCharge3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
//
// Calculates the correction due the Space Charge effect within the TPC drift volume
//
if (!fInitLookUp) {
AliInfo("Lookup table was not initialized! Performing the inizialisation now ...");
InitSpaceCharge3DDistortion();
}
Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
Double_t intEr, intEphi, intdEz ;
Double_t r, phi, z ;
Int_t sign;
r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
phi = TMath::ATan2(x[1],x[0]) ;
if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
z = x[2] ; // Create temporary copy of x[2]
if ( (roc%36) < 18 ) {
sign = 1; // (TPC A side)
} else {
sign = -1; // (TPC C side)
}
if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
// Get the Er and Ephi field integrals plus the integral over DeltaEz
intEr = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz );
intEphi = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
intdEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
// Calculate distorted position
if ( r > 0.0 ) {
phi = phi + fCorrectionFactor *( fC0*intEphi - fC1*intEr ) / r;
r = r + fCorrectionFactor *( fC0*intEr + fC1*intEphi );
}
Double_t dz = intdEz * fCorrectionFactor * fgkdvdE;
// Calculate correction in cartesian coordinates
dx[0] = - (r * TMath::Cos(phi) - x[0]);
dx[1] = - (r * TMath::Sin(phi) - x[1]);
dx[2] = - dz; // z distortion - (scaled with driftvelocity dependency on the Ez field and the overall scaling factor)
}
void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortion() {
//
// Initialization of the Lookup table which contains the solutions of the
// "space charge" (poisson) problem - Faster and more accureate
//
// Method: Weighted sum-up of the different fields within the look up table
// but using two lookup tables with higher granularity in the (r,z) and the (rphi)- plane to emulate
// more realistic space charges. (r,z) from primary ionisation. (rphi) for possible Gating leaks
if (fInitLookUp) {
AliInfo("Lookup table was already initialized! Doing it again anyway ...");
// return;
}
// ------------------------------------------------------------------------------------------------------
// step 1: lookup table in rz, fine grid, radial symetric, to emulate primary ionization
AliInfo("Step 1: Preparation of the weighted look-up tables.");
// lookup table in rz, fine grid
TFile *fZR = new TFile(fSCLookUpPOCsFileNameRZ.Data(),"READ");
if ( !fZR ) {
AliError("Precalculated POC-looup-table in ZR could not be found");
return;
}
// units are in [m]
TVector *gridf1 = (TVector*) fZR->Get("constants");
TVector &grid1 = *gridf1;
TMatrix *coordf1 = (TMatrix*) fZR->Get("coordinates");
TMatrix &coord1 = *coordf1;
TMatrix *coordPOCf1 = (TMatrix*) fZR->Get("POCcoord");
TMatrix &coordPOC1 = *coordPOCf1;
Int_t rows = (Int_t)grid1(0); // number of points in r direction - from RZ or RPhi table
Int_t phiSlices = (Int_t)grid1(1); // number of points in phi - from RPhi table
Int_t columns = (Int_t)grid1(2); // number of points in z direction - from RZ table
Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
// temporary matrices needed for the calculation // for rotational symmetric RZ table, phislices is 1
TMatrixD *arrayofErA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
TMatrixD *arrayofErC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
arrayofErA[k] = new TMatrixD(rows,columns) ;
arrayofdEzA[k] = new TMatrixD(rows,columns) ;
arrayofErC[k] = new TMatrixD(rows,columns) ;
arrayofdEzC[k] = new TMatrixD(rows,columns) ;
arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
// zero initialization not necessary, it is done in the constructor of TMatrix
}
// list of points as used during sum up
Double_t rlist1[kNRows], zedlist1[kNColumns];// , philist1[phiSlices];
for ( Int_t i = 0 ; i < rows ; i++ ) {
rlist1[i] = fgkIFCRadius + i*gridSizeR ;
for ( Int_t j = 0 ; j < columns ; j++ ) {
zedlist1[j] = j * gridSizeZ ;
}
}
TTree *treePOC = (TTree*)fZR->Get("POCall");
TVector *bEr = 0; //TVector *bEphi= 0;
TVector *bEz = 0;
treePOC->SetBranchAddress("Er",&bEr);
treePOC->SetBranchAddress("Ez",&bEz);
// Read the complete tree and do a weighted sum-up over the POC configurations
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
for (Int_t itreepC=0; itreepCGetEntry(itreepC);
// center of the POC voxel in [meter]
Double_t r0 = coordPOC1(ipC,0);
Double_t phi0 = coordPOC1(ipC,1);
Double_t z0 = coordPOC1(ipC,2);
ipC++; // POC configuration counter
// weights (charge density) at POC position on the A and C side (in C/m^3/e0)
// note: coordinates are in [cm]
Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100, 1); // partial load in r,z
Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100, 1); // partial load in r,z
// Summing up the vector components according to their weight
Int_t ip = 0;
for ( Int_t j = 0 ; j < columns ; j++ ) {
for ( Int_t i = 0 ; i < rows ; i++ ) {
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
// check wether the coordinates were screwed
if (TMath::Abs((coord1(0,ip)*100-rlist1[i]))>1 ||
TMath::Abs((coord1(2,ip)*100-zedlist1[j])>1)) {
AliError("internal error: coordinate system was screwed during the sum-up");
printf("sum-up: (r,z)=(%f,%f)\n",rlist1[i],zedlist1[j]);
printf("lookup: (r,z)=(%f,%f)\n",coord1(0,ip)*100,coord1(2,ip)*100);
AliError("Don't trust the results of the space charge calculation!");
}
// unfortunately, the lookup tables were produced to be faster for phi symmetric charges
// This will be the most frequent usage (hopefully)
// That's why we have to do this here ...
TMatrixD &erA = *arrayofErA[k] ;
TMatrixD &dEzA = *arrayofdEzA[k] ;
TMatrixD &erC = *arrayofErC[k] ;
TMatrixD &dEzC = *arrayofdEzC[k] ;
// Sum up - Efield values in [V/m] -> transition to [V/cm]
erA(i,j) += ((*bEr)(ip)) * weightA /100;
erC(i,j) += ((*bEr)(ip)) * weightC /100;
dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
// increase the counter
ip++;
}
}
} // end coordinate loop
} // end POC loop
// -------------------------------------------------------------------------------
// Division by the Ez (drift) field and integration along z
// AliInfo("Step 1: Division and integration");
Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
// matrices holding the solution - summation of POC charges // see above
TMatrixD &erA = *arrayofErA[k] ;
TMatrixD &dezA = *arrayofdEzA[k] ;
TMatrixD &erC = *arrayofErC[k] ;
TMatrixD &dezC = *arrayofdEzC[k] ;
// matrices which will contain the integrated fields (divided by the drift field)
TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
// Count backwards to facilitate integration over Z
Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point
// by trapezoidal rule.
erOverEzA(i,j) = 0; //ephiOverEzA(i,j) = 0;
deltaEzA(i,j) = 0;
erOverEzC(i,j) = 0; //ephiOverEzC(i,j) = 0;
deltaEzC(i,j) = 0;
for ( Int_t m = j ; m < columns ; m++ ) { // integration
erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
if ( index != 4 ) index = 4; else index = 2 ;
}
if ( index == 4 ) {
erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
}
if ( index == 2 ) {
erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
}
if ( j == columns-2 ) {
erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
}
if ( j == columns-1 ) {
erOverEzA(i,j) = 0;
erOverEzC(i,j) = 0;
deltaEzA(i,j) = 0;
deltaEzC(i,j) = 0;
}
}
}
}
// AliInfo("Step 1: Interpolation to Standard grid");
// -------------------------------------------------------------------------------
// Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
Double_t r, z;//phi, z ;
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
// phi = fgkPhiList[k] ;
// final lookup table
TMatrixD &erOverEzFinal = *fLookUpErOverEz[k] ;
TMatrixD &deltaEzFinal = *fLookUpDeltaEz[k] ;
// calculated and integrated tables - just one phi slice
TMatrixD &erOverEzA = *arrayofEroverEzA[0] ;
TMatrixD &deltaEzA = *arrayofDeltaEzA[0];
TMatrixD &erOverEzC = *arrayofEroverEzC[0] ;
TMatrixD &deltaEzC = *arrayofDeltaEzC[0];
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
for ( Int_t i = 0 ; i < kNR ; i++ ) {
r = fgkRList[i] ;
// Interpolate Lookup tables onto standard grid
if (fgkZList[j]>0) {
erOverEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzA );
deltaEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzA );
} else {
erOverEzFinal(i,j) = Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, erOverEzC );
deltaEzFinal(i,j) = - Interpolate2DTable(order, r, z, rows, columns, rlist1, zedlist1, deltaEzC );
// negative coordinate system on C side
}
} // end r loop
} // end z loop
} // end phi loop
// clear the temporary arrays lists
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
delete arrayofErA[k];
delete arrayofdEzA[k];
delete arrayofErC[k];
delete arrayofdEzC[k];
delete arrayofEroverEzA[k];
delete arrayofDeltaEzA[k];
delete arrayofEroverEzC[k];
delete arrayofDeltaEzC[k];
}
fZR->Close();
// ------------------------------------------------------------------------------------------------------
// Step 2: Load and sum up lookup table in rphi, fine grid, to emulate for example a GG leak
// AliInfo("Step 2: Preparation of the weighted look-up table");
TFile *fRPhi = new TFile(fSCLookUpPOCsFileNameRPhi.Data(),"READ");
if ( !fRPhi ) {
AliError("Precalculated POC-looup-table in RPhi could not be found");
return;
}
// units are in [m]
TVector *gridf2 = (TVector*) fRPhi->Get("constants");
TVector &grid2 = *gridf2;
TMatrix *coordf2 = (TMatrix*) fRPhi->Get("coordinates");
TMatrix &coord2 = *coordf2;
TMatrix *coordPOCf2 = (TMatrix*) fRPhi->Get("POCcoord");
TMatrix &coordPOC2 = *coordPOCf2;
rows = (Int_t)grid2(0); // number of points in r direction
phiSlices = (Int_t)grid2(1); // number of points in phi
columns = (Int_t)grid2(2); // number of points in z direction
gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
Float_t gridSizePhi = TMath::TwoPi()/phiSlices; // unit in [rad]
gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
// list of points as used during sum up
Double_t rlist2[kNRows], philist2[kNPhiSlices], zedlist2[kNColumns];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
philist2[k] = gridSizePhi * k;
for ( Int_t i = 0 ; i < rows ; i++ ) {
rlist2[i] = fgkIFCRadius + i*gridSizeR ;
for ( Int_t j = 0 ; j < columns ; j++ ) {
zedlist2[j] = j * gridSizeZ ;
}
}
} // only done once
// temporary matrices needed for the calculation
TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
TMatrixD *arrayofErC2[kNPhiSlices], *arrayofEphiC2[kNPhiSlices], *arrayofdEzC2[kNPhiSlices];
TMatrixD *arrayofEroverEzA2[kNPhiSlices], *arrayofEphioverEzA2[kNPhiSlices], *arrayofDeltaEzA2[kNPhiSlices];
TMatrixD *arrayofEroverEzC2[kNPhiSlices], *arrayofEphioverEzC2[kNPhiSlices], *arrayofDeltaEzC2[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
arrayofErA2[k] = new TMatrixD(rows,columns) ;
arrayofEphiA2[k] = new TMatrixD(rows,columns) ;
arrayofdEzA2[k] = new TMatrixD(rows,columns) ;
arrayofErC2[k] = new TMatrixD(rows,columns) ;
arrayofEphiC2[k] = new TMatrixD(rows,columns) ;
arrayofdEzC2[k] = new TMatrixD(rows,columns) ;
arrayofEroverEzA2[k] = new TMatrixD(rows,columns) ;
arrayofEphioverEzA2[k] = new TMatrixD(rows,columns) ;
arrayofDeltaEzA2[k] = new TMatrixD(rows,columns) ;
arrayofEroverEzC2[k] = new TMatrixD(rows,columns) ;
arrayofEphioverEzC2[k] = new TMatrixD(rows,columns) ;
arrayofDeltaEzC2[k] = new TMatrixD(rows,columns) ;
// zero initialization not necessary, it is done in the constructor of TMatrix
}
treePOC = (TTree*)fRPhi->Get("POCall");
// TVector *bEr = 0; // done above
TVector *bEphi= 0;
// TVector *bEz = 0; // done above
treePOC->SetBranchAddress("Er",&bEr);
treePOC->SetBranchAddress("Ephi",&bEphi);
treePOC->SetBranchAddress("Ez",&bEz);
// Read the complete tree and do a weighted sum-up over the POC configurations
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
for (Int_t itreepC=0; itreepCGetEntry(itreepC);
// center of the POC voxel in [meter]
Double_t r0 = coordPOC2(ipC,0);
Double_t phi0 = coordPOC2(ipC,1);
// Double_t z0 = coordPOC2(ipC,2);
// weights (charge density) at POC position on the A and C side (in C/m^3/e0)
// note: coordinates are in [cm]
Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, 0.499, 2); // partial load in r,phi
Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-0.499, 2); // partial load in r,phi
// printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
// Summing up the vector components according to their weight
Int_t ip = 0;
for ( Int_t j = 0 ; j < columns ; j++ ) {
for ( Int_t i = 0 ; i < rows ; i++ ) {
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
// check wether the coordinates were screwed
if (TMath::Abs((coord2(0,ip)*100-rlist2[i]))>1 ||
TMath::Abs((coord2(1,ip)-philist2[k]))>1 ||
TMath::Abs((coord2(2,ip)*100-zedlist2[j]))>1) {
AliError("internal error: coordinate system was screwed during the sum-up");
printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist2[i],philist2[k],zedlist2[j]);
AliError("Don't trust the results of the space charge calculation!");
}
// unfortunately, the lookup tables were produced to be faster for phi symmetric charges
// This will be the most frequent usage (hopefully)
// That's why we have to do this here ...
TMatrixD &erA = *arrayofErA2[k] ;
TMatrixD &ephiA = *arrayofEphiA2[k];
TMatrixD &dEzA = *arrayofdEzA2[k] ;
TMatrixD &erC = *arrayofErC2[k] ;
TMatrixD &ephiC = *arrayofEphiC2[k];
TMatrixD &dEzC = *arrayofdEzC2[k] ;
// Sum up - Efield values in [V/m] -> transition to [V/cm]
erA(i,j) += ((*bEr)(ip)) * weightA /100;
erC(i,j) += ((*bEr)(ip)) * weightC /100;
ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
// increase the counter
ip++;
}
}
} // end coordinate loop
// Rotation and summation in the rest of the dPhiSteps
// which were not stored in the this tree due to storage & symmetry reasons
Int_t phiPoints = (Int_t) grid2(1);
Int_t phiPOC = (Int_t) grid2(4);
// printf("%d %d\n",phiPOC,flagRadSym);
for (Int_t phiiC = 1; phiiC=rotVal) {
ipR = ip-rotVal;
} else {
ipR = ip+(phiPoints-rotVal);
}
// unfortunately, the lookup tables were produced to be faster for phi symmetric charges
// This will be the most frequent usage
// That's why we have to do this here and not outside the loop ...
TMatrixD &erA = *arrayofErA2[k] ;
TMatrixD &ephiA = *arrayofEphiA2[k];
TMatrixD &dEzA = *arrayofdEzA2[k] ;
TMatrixD &erC = *arrayofErC2[k] ;
TMatrixD &ephiC = *arrayofEphiC2[k];
TMatrixD &dEzC = *arrayofdEzC2[k] ;
// Sum up - Efield values in [V/m] -> transition to [V/cm]
erA(i,j) += ((*bEr)(ipR)) * weightA /100;
erC(i,j) += ((*bEr)(ipR)) * weightC /100;
ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
// increase the counter
ip++;
}
}
} // end coordinate loop
} // end phi-POC summation (phiiC)
ipC++; // POC configuration counter
// printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
}
// -------------------------------------------------------------------------------
// Division by the Ez (drift) field and integration along z
// AliInfo("Step 2: Division and integration");
for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
// matrices holding the solution - summation of POC charges // see above
TMatrixD &erA = *arrayofErA2[k] ;
TMatrixD &ephiA = *arrayofEphiA2[k];
TMatrixD &dezA = *arrayofdEzA2[k] ;
TMatrixD &erC = *arrayofErC2[k] ;
TMatrixD &ephiC = *arrayofEphiC2[k];
TMatrixD &dezC = *arrayofdEzC2[k] ;
// matrices which will contain the integrated fields (divided by the drift field)
TMatrixD &erOverEzA = *arrayofEroverEzA2[k] ;
TMatrixD &ephiOverEzA = *arrayofEphioverEzA2[k];
TMatrixD &deltaEzA = *arrayofDeltaEzA2[k];
TMatrixD &erOverEzC = *arrayofEroverEzC2[k] ;
TMatrixD &ephiOverEzC = *arrayofEphioverEzC2[k];
TMatrixD &deltaEzC = *arrayofDeltaEzC2[k];
for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
// Count backwards to facilitate integration over Z
Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
erOverEzA(i,j) = 0;
ephiOverEzA(i,j) = 0;
deltaEzA(i,j) = 0;
erOverEzC(i,j) = 0;
ephiOverEzC(i,j) = 0;
deltaEzC(i,j) = 0;
for ( Int_t m = j ; m < columns ; m++ ) { // integration
erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
if ( index != 4 ) index = 4; else index = 2 ;
}
if ( index == 4 ) {
erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
}
if ( index == 2 ) {
erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
}
if ( j == columns-2 ) {
erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
}
if ( j == columns-1 ) {
erOverEzA(i,j) = 0;
erOverEzC(i,j) = 0;
ephiOverEzA(i,j) = 0;
ephiOverEzC(i,j) = 0;
deltaEzA(i,j) = 0;
deltaEzC(i,j) = 0;
}
}
}
}
AliInfo("Step 2: Interpolation to Standard grid");
// -------------------------------------------------------------------------------
// Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
Double_t phi = fgkPhiList[k] ;
// final lookup table
TMatrixD &erOverEzFinal = *fLookUpErOverEz[k] ;
TMatrixD &ephiOverEzFinal = *fLookUpEphiOverEz[k];
TMatrixD &deltaEzFinal = *fLookUpDeltaEz[k] ;
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
for ( Int_t i = 0 ; i < kNR ; i++ ) {
r = fgkRList[i] ;
// Interpolate Lookup tables onto standard grid
if (fgkZList[j]>0) {
erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist2, zedlist2, philist2, arrayofEroverEzA2 );
ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist2, zedlist2, philist2, arrayofEphioverEzA2);
deltaEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist2, zedlist2, philist2, arrayofDeltaEzA2 );
} else {
erOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist2, zedlist2, philist2, arrayofEroverEzC2 );
ephiOverEzFinal(i,j) += Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist2, zedlist2, philist2, arrayofEphioverEzC2);
deltaEzFinal(i,j) -= Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist2, zedlist2, philist2, arrayofDeltaEzC2 );
}
} // end r loop
} // end z loop
} // end phi loop
// clear the temporary arrays lists
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
delete arrayofErA2[k];
delete arrayofEphiA2[k];
delete arrayofdEzA2[k];
delete arrayofErC2[k];
delete arrayofEphiC2[k];
delete arrayofdEzC2[k];
delete arrayofEroverEzA2[k];
delete arrayofEphioverEzA2[k];
delete arrayofDeltaEzA2[k];
delete arrayofEroverEzC2[k];
delete arrayofEphioverEzC2[k];
delete arrayofDeltaEzC2[k];
}
fRPhi->Close();
// FINISHED
fInitLookUp = kTRUE;
}
void AliTPCSpaceCharge3D::InitSpaceCharge3DDistortionCourse() {
//
// Initialization of the Lookup table which contains the solutions of the
// "space charge" (poisson) problem
//
// The sum-up uses a look-up table which contains different discretized Space charge fields
// in order to calculate the corresponding field deviations due to a given (discretized)
// space charge distribution ....
//
// Method of calculation: Weighted sum-up of the different fields within the look up table
// Note: Full 3d version: Course and slow ...
if (fInitLookUp) {
AliInfo("Lookup table was already initialized!");
// return;
}
AliInfo("Preparation of the weighted look-up table");
TFile *f = new TFile(fSCLookUpPOCsFileName3D.Data(),"READ");
if ( !f ) {
AliError("Precalculated POC-looup-table could not be found");
return;
}
// units are in [m]
TVector *gridf = (TVector*) f->Get("constants");
TVector &grid = *gridf;
TMatrix *coordf = (TMatrix*) f->Get("coordinates");
TMatrix &coord = *coordf;
TMatrix *coordPOCf = (TMatrix*) f->Get("POCcoord");
TMatrix &coordPOC = *coordPOCf;
Bool_t flagRadSym = 0;
if (grid(1)==1 && grid(4)==1) {
// AliInfo("LOOK UP TABLE IS RADIAL SYMETTRIC - Field in Phi is ZERO");
flagRadSym=1;
}
Int_t rows = (Int_t)grid(0); // number of points in r direction
Int_t phiSlices = (Int_t)grid(1); // number of points in phi
Int_t columns = (Int_t)grid(2); // number of points in z direction
const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius)/(rows-1); // unit in [cm]
const Float_t gridSizePhi = TMath::TwoPi()/phiSlices; // unit in [rad]
const Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
// temporary matrices needed for the calculation
TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
arrayofErA[k] = new TMatrixD(rows,columns) ;
arrayofEphiA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
arrayofdEzA[k] = new TMatrixD(rows,columns) ;
arrayofErC[k] = new TMatrixD(rows,columns) ;
arrayofEphiC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
arrayofdEzC[k] = new TMatrixD(rows,columns) ;
arrayofEroverEzA[k] = new TMatrixD(rows,columns) ;
arrayofEphioverEzA[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
arrayofDeltaEzA[k] = new TMatrixD(rows,columns) ;
arrayofEroverEzC[k] = new TMatrixD(rows,columns) ;
arrayofEphioverEzC[k] = new TMatrixD(rows,columns) ; // zero if radial symmetric
arrayofDeltaEzC[k] = new TMatrixD(rows,columns) ;
// Set the values to zero the lookup tables
// not necessary, it is done in the constructor of TMatrix - code deleted
}
// list of points as used in the interpolation (during sum up)
Double_t rlist[kNRows], zedlist[kNColumns] , philist[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
philist[k] = gridSizePhi * k;
for ( Int_t i = 0 ; i < rows ; i++ ) {
rlist[i] = fgkIFCRadius + i*gridSizeR ;
for ( Int_t j = 0 ; j < columns ; j++ ) {
zedlist[j] = j * gridSizeZ ;
}
}
} // only done once
TTree *treePOC = (TTree*)f->Get("POCall");
TVector *bEr = 0; TVector *bEphi= 0; TVector *bEz = 0;
treePOC->SetBranchAddress("Er",&bEr);
if (!flagRadSym) treePOC->SetBranchAddress("Ephi",&bEphi);
treePOC->SetBranchAddress("Ez",&bEz);
// Read the complete tree and do a weighted sum-up over the POC configurations
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t treeNumPOC = (Int_t)treePOC->GetEntries(); // Number of POC conf. in the look-up table
Int_t ipC = 0; // POC Conf. counter (note: different to the POC number in the tree!)
for (Int_t itreepC=0; itreepCGetEntry(itreepC);
// center of the POC voxel in [meter]
Double_t r0 = coordPOC(ipC,0);
Double_t phi0 = coordPOC(ipC,1);
Double_t z0 = coordPOC(ipC,2);
ipC++; // POC configuration counter
// weights (charge density) at POC position on the A and C side (in C/m^3/e0)
// note: coordinates are in [cm]
Double_t weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
Double_t weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
// Summing up the vector components according to their weight
Int_t ip = 0;
for ( Int_t j = 0 ; j < columns ; j++ ) {
for ( Int_t i = 0 ; i < rows ; i++ ) {
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
// check wether the coordinates were screwed
if (TMath::Abs((coord(0,ip)*100-rlist[i]))>1 ||
TMath::Abs((coord(1,ip)-philist[k]))>1 ||
TMath::Abs((coord(2,ip)*100-zedlist[j]))>1) {
AliError("internal error: coordinate system was screwed during the sum-up");
printf("lookup: (r,phi,z)=(%f,%f,%f)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
printf("sum-up: (r,phi,z)=(%f,%f,%f)\n",rlist[i],philist[k],zedlist[j]);
AliError("Don't trust the results of the space charge calculation!");
}
// unfortunately, the lookup tables were produced to be faster for phi symmetric charges
// This will be the most frequent usage (hopefully)
// That's why we have to do this here ...
TMatrixD &erA = *arrayofErA[k] ;
TMatrixD &ephiA = *arrayofEphiA[k];
TMatrixD &dEzA = *arrayofdEzA[k] ;
TMatrixD &erC = *arrayofErC[k] ;
TMatrixD &ephiC = *arrayofEphiC[k];
TMatrixD &dEzC = *arrayofdEzC[k] ;
// Sum up - Efield values in [V/m] -> transition to [V/cm]
erA(i,j) += ((*bEr)(ip)) * weightA /100;
erC(i,j) += ((*bEr)(ip)) * weightC /100;
if (!flagRadSym) {
ephiA(i,j) += ((*bEphi)(ip)) * weightA/100; // [V/rad]
ephiC(i,j) += ((*bEphi)(ip)) * weightC/100; // [V/rad]
}
dEzA(i,j) += ((*bEz)(ip)) * weightA /100;
dEzC(i,j) += ((*bEz)(ip)) * weightC /100;
// increase the counter
ip++;
}
}
} // end coordinate loop
// Rotation and summation in the rest of the dPhiSteps
// which were not stored in the this tree due to storage & symmetry reasons
Int_t phiPoints = (Int_t) grid(1);
Int_t phiPOC = (Int_t) grid(4);
// printf("%d %d\n",phiPOC,flagRadSym);
for (Int_t phiiC = 1; phiiC=rotVal) {
ipR = ip-rotVal;
} else {
ipR = ip+(phiPoints-rotVal);
}
// unfortunately, the lookup tables were produced to be faster for phi symmetric charges
// This will be the most frequent usage
// That's why we have to do this here and not outside the loop ...
TMatrixD &erA = *arrayofErA[k] ;
TMatrixD &ephiA = *arrayofEphiA[k];
TMatrixD &dEzA = *arrayofdEzA[k] ;
TMatrixD &erC = *arrayofErC[k] ;
TMatrixD &ephiC = *arrayofEphiC[k];
TMatrixD &dEzC = *arrayofdEzC[k] ;
// Sum up - Efield values in [V/m] -> transition to [V/cm]
erA(i,j) += ((*bEr)(ipR)) * weightA /100;
erC(i,j) += ((*bEr)(ipR)) * weightC /100;
if (!flagRadSym) {
ephiA(i,j) += ((*bEphi)(ipR)) * weightA/100; // [V/rad]
ephiC(i,j) += ((*bEphi)(ipR)) * weightC/100; // [V/rad]
}
dEzA(i,j) += ((*bEz)(ipR)) * weightA /100;
dEzC(i,j) += ((*bEz)(ipR)) * weightC /100;
// increase the counter
ip++;
}
}
} // end coordinate loop
} // end phi-POC summation (phiiC)
// printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
}
// -------------------------------------------------------------------------------
// Division by the Ez (drift) field and integration along z
AliInfo("Division and integration");
Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm;
for ( Int_t k = 0 ; k < phiSlices ; k++ ) { // phi loop
// matrices holding the solution - summation of POC charges // see above
TMatrixD &erA = *arrayofErA[k] ;
TMatrixD &ephiA = *arrayofEphiA[k];
TMatrixD &dezA = *arrayofdEzA[k] ;
TMatrixD &erC = *arrayofErC[k] ;
TMatrixD &ephiC = *arrayofEphiC[k];
TMatrixD &dezC = *arrayofdEzC[k] ;
// matrices which will contain the integrated fields (divided by the drift field)
TMatrixD &erOverEzA = *arrayofEroverEzA[k] ;
TMatrixD &ephiOverEzA = *arrayofEphioverEzA[k];
TMatrixD &deltaEzA = *arrayofDeltaEzA[k];
TMatrixD &erOverEzC = *arrayofEroverEzC[k] ;
TMatrixD &ephiOverEzC = *arrayofEphioverEzC[k];
TMatrixD &deltaEzC = *arrayofDeltaEzC[k];
for ( Int_t i = 0 ; i < rows ; i++ ) { // r loop
for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {// z loop
// Count backwards to facilitate integration over Z
Int_t index = 1 ; // Simpsons rule if N=odd.If N!=odd then add extra point by trapezoidal rule.
erOverEzA(i,j) = 0; ephiOverEzA(i,j) = 0; deltaEzA(i,j) = 0;
erOverEzC(i,j) = 0; ephiOverEzC(i,j) = 0; deltaEzC(i,j) = 0;
for ( Int_t m = j ; m < columns ; m++ ) { // integration
erOverEzA(i,j) += index*(gridSizeZ/3.0)*erA(i,m)/(-1*ezField) ;
erOverEzC(i,j) += index*(gridSizeZ/3.0)*erC(i,m)/(-1*ezField) ;
if (!flagRadSym) {
ephiOverEzA(i,j) += index*(gridSizeZ/3.0)*ephiA(i,m)/(-1*ezField) ;
ephiOverEzC(i,j) += index*(gridSizeZ/3.0)*ephiC(i,m)/(-1*ezField) ;
}
deltaEzA(i,j) += index*(gridSizeZ/3.0)*dezA(i,m)/(-1) ;
deltaEzC(i,j) += index*(gridSizeZ/3.0)*dezC(i,m)/(-1) ;
if ( index != 4 ) index = 4; else index = 2 ;
}
if ( index == 4 ) {
erOverEzA(i,j) -= (gridSizeZ/3.0)*erA(i,columns-1)/(-1*ezField) ;
erOverEzC(i,j) -= (gridSizeZ/3.0)*erC(i,columns-1)/(-1*ezField) ;
if (!flagRadSym) {
ephiOverEzA(i,j) -= (gridSizeZ/3.0)*ephiA(i,columns-1)/(-1*ezField) ;
ephiOverEzC(i,j) -= (gridSizeZ/3.0)*ephiC(i,columns-1)/(-1*ezField) ;
}
deltaEzA(i,j) -= (gridSizeZ/3.0)*dezA(i,columns-1)/(-1) ;
deltaEzC(i,j) -= (gridSizeZ/3.0)*dezC(i,columns-1)/(-1) ;
}
if ( index == 2 ) {
erOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*erA(i,columns-2)-2.5*erA(i,columns-1))/(-1*ezField) ;
erOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*erC(i,columns-2)-2.5*erC(i,columns-1))/(-1*ezField) ;
if (!flagRadSym) {
ephiOverEzA(i,j) += (gridSizeZ/3.0)*(0.5*ephiA(i,columns-2)-2.5*ephiA(i,columns-1))/(-1*ezField) ;
ephiOverEzC(i,j) += (gridSizeZ/3.0)*(0.5*ephiC(i,columns-2)-2.5*ephiC(i,columns-1))/(-1*ezField) ;
}
deltaEzA(i,j) += (gridSizeZ/3.0)*(0.5*dezA(i,columns-2)-2.5*dezA(i,columns-1))/(-1) ;
deltaEzC(i,j) += (gridSizeZ/3.0)*(0.5*dezC(i,columns-2)-2.5*dezC(i,columns-1))/(-1) ;
}
if ( j == columns-2 ) {
erOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*erA(i,columns-2)+1.5*erA(i,columns-1))/(-1*ezField) ;
erOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*erC(i,columns-2)+1.5*erC(i,columns-1))/(-1*ezField) ;
if (!flagRadSym) {
ephiOverEzA(i,j) = (gridSizeZ/3.0)*(1.5*ephiA(i,columns-2)+1.5*ephiA(i,columns-1))/(-1*ezField) ;
ephiOverEzC(i,j) = (gridSizeZ/3.0)*(1.5*ephiC(i,columns-2)+1.5*ephiC(i,columns-1))/(-1*ezField) ;
}
deltaEzA(i,j) = (gridSizeZ/3.0)*(1.5*dezA(i,columns-2)+1.5*dezA(i,columns-1))/(-1) ;
deltaEzC(i,j) = (gridSizeZ/3.0)*(1.5*dezC(i,columns-2)+1.5*dezC(i,columns-1))/(-1) ;
}
if ( j == columns-1 ) {
erOverEzA(i,j) = 0;
erOverEzC(i,j) = 0;
if (!flagRadSym) {
ephiOverEzA(i,j) = 0;
ephiOverEzC(i,j) = 0;
}
deltaEzA(i,j) = 0;
deltaEzC(i,j) = 0;
}
}
}
}
AliInfo("Interpolation to Standard grid");
// -------------------------------------------------------------------------------
// Interpolate results onto the standard grid which is used for all AliTPCCorrections classes
const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
Double_t r, phi, z ;
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
phi = fgkPhiList[k] ;
TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
z = TMath::Abs(fgkZList[j]) ; // z position is symmetric
for ( Int_t i = 0 ; i < kNR ; i++ ) {
r = fgkRList[i] ;
// Interpolate Lookup tables onto standard grid
if (fgkZList[j]>0) {
erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist, zedlist, philist, arrayofEroverEzA );
ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist, zedlist, philist, arrayofEphioverEzA);
deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist, zedlist, philist, arrayofDeltaEzA );
} else {
erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist, zedlist, philist, arrayofEroverEzC );
ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist, zedlist, philist, arrayofEphioverEzC);
deltaEz(i,j) = - Interpolate3DTable(order, r, z, phi, rows, columns, phiSlices,
rlist, zedlist, philist, arrayofDeltaEzC );
// negative coordinate system on C side
}
} // end r loop
} // end z loop
} // end phi loop
// clear the temporary arrays lists
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
delete arrayofErA[k];
delete arrayofEphiA[k];
delete arrayofdEzA[k];
delete arrayofErC[k];
delete arrayofEphiC[k];
delete arrayofdEzC[k];
delete arrayofEroverEzA[k];
delete arrayofEphioverEzA[k];
delete arrayofDeltaEzA[k];
delete arrayofEroverEzC[k];
delete arrayofEphioverEzC[k];
delete arrayofDeltaEzC[k];
}
fInitLookUp = kTRUE;
}
void AliTPCSpaceCharge3D::SetSCDataFileName(TString fname) {
//
// Set & load the Space charge density distribution from a file
// (linear interpolation onto a standard grid)
//
fSCDataFileName = fname;
TFile *f = new TFile(fSCDataFileName.Data(),"READ");
if (!f) {
AliError(Form("File %s, which should contain the space charge distribution, could not be found",
fSCDataFileName.Data()));
return;
}
TH2F *densityRZ = (TH2F*) f->Get("SpaceChargeInRZ");
if (!densityRZ) {
AliError(Form("The indicated file (%s) does not contain a histogram called %s",
fSCDataFileName.Data(),"SpaceChargeInRZ"));
return;
}
TH3F *densityRPhi = (TH3F*) f->Get("SpaceChargeInRPhi");
if (!densityRPhi) {
AliError(Form("The indicated file (%s) does not contain a histogram called %s",
fSCDataFileName.Data(),"SpaceChargeInRPhi"));
return;
}
Double_t r, phi, z ;
TMatrixD &scDensityInRZ = *fSCdensityInRZ;
TMatrixD &scDensityInRPhiA = *fSCdensityInRPhiA;
TMatrixD &scDensityInRPhiC = *fSCdensityInRPhiC;
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
phi = fgkPhiList[k] ;
TMatrixD &scDensity = *fSCdensityDistribution[k] ;
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
z = fgkZList[j] ;
for ( Int_t i = 0 ; i < kNR ; i++ ) {
r = fgkRList[i] ;
// partial load in (r,z)
if (k==0) // do just once
scDensityInRZ(i,j) = densityRZ->Interpolate(r,z);
// partial load in (r,phi)
if ( j==0 || j == kNZ/2 ) {
if (z>0)
scDensityInRPhiA(i,k) = densityRPhi->Interpolate(r,phi,0.499); // A side
else
scDensityInRPhiC(i,k) = densityRPhi->Interpolate(r,phi,-0.499); // C side
}
// Full 3D configuration ...
if (z>0)
scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiA(i,k);
else
scDensity(i,j) = scDensityInRZ(i,j) + scDensityInRPhiC(i,k);
}
}
}
f->Close();
fInitLookUp = kFALSE;
}
Float_t AliTPCSpaceCharge3D::GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z, Int_t mode) {
//
// returns the (input) space charge density at a given point according
// Note: input in [cm], output in [C/m^3/e0] !!
//
while (phi<0) phi += TMath::TwoPi();
while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
// Float_t sc =fSCdensityDistribution->Interpolate(r0,phi0,z0);
Int_t order = 1; //
Float_t sc = 0;
if (mode == 0) { // return full load
sc = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
fgkRList, fgkZList, fgkPhiList, fSCdensityDistribution );
} else if (mode == 1) { // return partial load in (r,z)
TMatrixD &scDensityInRZ = *fSCdensityInRZ;
sc = Interpolate2DTable(order, r, z, kNR, kNZ, fgkRList, fgkZList, scDensityInRZ );
} else if (mode == 2) { // return partial load in (r,phi)
if (z>0) {
TMatrixD &scDensityInRPhi = *fSCdensityInRPhiA;
sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
} else {
TMatrixD &scDensityInRPhi = *fSCdensityInRPhiC;
sc = Interpolate2DTable(order, r, phi, kNR, kNPhi, fgkRList, fgkPhiList, scDensityInRPhi );
}
} else {
// should i give a warning?
sc = 0;
}
// printf("%f %f %f: %f\n",r,phi,z,sc);
return sc;
}
TH2F * AliTPCSpaceCharge3D::CreateHistoSCinXY(Float_t z, Int_t nx, Int_t ny, Int_t mode) {
//
// return a simple histogramm containing the space charge distribution (input for the calculation)
//
TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"x [cm]","y [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
nx,-250.,250.,ny,-250.,250.);
for (Int_t iy=1;iy<=ny;++iy) {
Double_t yp = h->GetYaxis()->GetBinCenter(iy);
for (Int_t ix=1;ix<=nx;++ix) {
Double_t xp = h->GetXaxis()->GetBinCenter(ix);
Float_t r = TMath::Sqrt(xp*xp+yp*yp);
Float_t phi = TMath::ATan2(yp,xp);
if (85.<=r && r<=250.) {
Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
h->SetBinContent(ix,iy,sc);
} else {
h->SetBinContent(ix,iy,0.);
}
}
}
return h;
}
TH2F * AliTPCSpaceCharge3D::CreateHistoSCinZR(Float_t phi, Int_t nz, Int_t nr,Int_t mode ) {
//
// return a simple histogramm containing the space charge distribution (input for the calculation)
//
TH2F *h=CreateTH2F("spaceCharge",GetTitle(),"z [cm]","r [cm]","#rho_{sc} [C/m^{3}/e_{0}]",
nz,-250.,250.,nr,85.,250.);
for (Int_t ir=1;ir<=nr;++ir) {
Float_t r = h->GetYaxis()->GetBinCenter(ir);
for (Int_t iz=1;iz<=nz;++iz) {
Float_t z = h->GetXaxis()->GetBinCenter(iz);
Float_t sc = GetSpaceChargeDensity(r,phi,z,mode)/fgke0; // in [C/m^3/e0]
h->SetBinContent(iz,ir,sc);
}
}
return h;
}
void AliTPCSpaceCharge3D::WriteChargeDistributionToFile(const char* fname) {
//
// Example on how to write a Space charge distribution into a File
// (see below: estimate from scaling STAR measurements to Alice)
// Charge distribution is splitted into two (RZ and RPHI) in order to speed up
// the needed calculation time
//
TFile *f = new TFile(fname,"RECREATE");
// some grid, not too course
Int_t nr = 350;
Int_t nphi = 180;
Int_t nz = 500;
Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1);
Double_t dphi = TMath::TwoPi()/(nphi+1);
Double_t dz = 500./(nz+1);
Double_t safty = 0.; // due to a root bug which does not interpolate the boundary (first and last bin) correctly
// Charge distribution in ZR (rotational symmetric) ------------------
TH2F *histoZR = new TH2F("chargeZR","chargeZR",
nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
nz,-250-dz-safty,250+dz+safty);
for (Int_t ir=1;ir<=nr;++ir) {
Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir);
for (Int_t iz=1;iz<=nz;++iz) {
Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz);
// recalculation to meter
Double_t lZ = 2.5; // approx. TPC drift length
Double_t rpM = rp/100.; // in [m]
Double_t zpM = TMath::Abs(zp/100.); // in [m]
// setting of mb multiplicity and Interaction rate
Double_t multiplicity = 950;
Double_t intRate = 7800;
// calculation of "scaled" parameters
Double_t a = multiplicity*intRate/79175;
Double_t b = a/lZ;
Double_t charge = (a - b*zpM)/(rpM*rpM); // charge in [C/m^3/e0]
charge = charge*fgke0; // [C/m^3]
if (zp<0) charge *= 0.9; // e.g. slightly less on C side due to front absorber
// charge = 0; // for tests
histoZR->SetBinContent(ir,iz,charge);
}
}
histoZR->Write("SpaceChargeInRZ");
// Charge distribution in RPhi (e.g. Floating GG wire) ------------
TH3F *histoRPhi = new TH3F("chargeRPhi","chargeRPhi",
nr,fgkIFCRadius-dr-safty,fgkOFCRadius+dr+safty,
nphi,0-dphi-safty,TMath::TwoPi()+dphi+safty,
2,-1,1); // z part - to allow A and C side differences
// some 'arbitrary' GG leaks
Int_t nGGleaks = 5;
Double_t secPosA[5] = {3,6,6,11,13}; // sector
Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm
Double_t secPosC[5] = {1,8,12,15,15}; // sector
Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm
for (Int_t ir=1;ir<=nr;++ir) {
Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir);
for (Int_t iphi=1;iphi<=nphi;++iphi) {
Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi);
for (Int_t iz=1;iz<=2;++iz) {
Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz);
Double_t charge = 0;
for (Int_t igg = 0; igg(TMath::Pi()/9*secPos) ) ) { // sector slice
if ( rp>(radialPos-2.5) && rp<(radialPos+2.5)) // 5 cm slice
charge = 300;
}
}
charge = charge*fgke0; // [C/m^3]
histoRPhi->SetBinContent(ir,iphi,iz,charge);
}
}
}
histoRPhi->Write("SpaceChargeInRPhi");
f->Close();
}
void AliTPCSpaceCharge3D::Print(const Option_t* option) const {
//
// Print function to check the settings of the boundary vectors
// option=="a" prints the C0 and C1 coefficents for calibration purposes
//
TString opt = option; opt.ToLower();
printf("%s\n",GetTitle());
printf(" - Space Charge effect with arbitrary 3D charge density (as input).\n");
printf(" SC correction factor: %f \n",fCorrectionFactor);
if (opt.Contains("a")) { // Print all details
printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
}
if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitSpaceCharge3DDistortion() ...");
}