* provided "as is" without express or implied warranty. *
**************************************************************************/
-//////////////////////////////////////////////////////////////////////////////
-// //
-// AliTPCSpaceCharge3D class //
-// The class calculates the space point distortions due to a space charge //
-// effect .... //
-// Method of calculation: //
-// The analytical solution for the poisson problem in 3D (cylindrical coord)//
-// is used in form of look up tables. PieceOfCake (POC) voxel were pre- //
-// calculated and can be sumed up (weighted) according to the what is needed//
-// //
-// The class allows "effective Omega Tau" corrections. //
-// //
-// NOTE: This class is not capable of calculation z distortions due to //
-// drift velocity change in dependence of the electric field!!! //
-// //
-// date: 19/06/2010 //
-// Authors: Stefan Rossegger //
-// //
-// Example usage: //
-//////////////////////////////////////////////////////////////////////////////
+// _________________________________________________________________
+//
+// Begin_Html
+// <h2> AliTPCSpaceCharge3D class </h2>
+// The class calculates the space point distortions due to an arbitrary space
+// charge distribution in 3D.
+// <p>
+// 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}.
+// <p>
+// 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.
+// <p>
+// 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("cAliTPCSpaceCharge3D","cAliTPCSpaceCharge3D",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
+// <p>
+// Date: 19/06/2010 <br>
+// Authors: Stefan Rossegger
+// End_Html
+// _________________________________________________________________
+
#include "AliMagF.h"
#include "TGeoGlobalMagField.h"
fSCLookUpPOCsFileNameRZ(""),
fSCLookUpPOCsFileNameRPhi(""),
fSCdensityInRZ(0),
- fSCdensityInRPhiA(0),
+ fSCdensityInRPhiA(0),
fSCdensityInRPhiC(0)
{
//
// 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);
+ fLookUpErOverEz[k] = new TMatrixF(kNR,kNZ);
+ fLookUpEphiOverEz[k] = new TMatrixF(kNR,kNZ);
+ fLookUpDeltaEz[k] = new TMatrixF(kNR,kNZ);
+ fSCdensityDistribution[k] = new TMatrixF(kNR,kNZ);
}
fSCdensityInRZ = new TMatrixD(kNR,kNZ);
fSCdensityInRPhiA = new TMatrixD(kNR,kNPhi);
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-36-26_34p-18p-01p-MN40.root";
- // fSCLookUpPOCsFileNameRPhi="$(ALICE_ROOT)/TPC/Calib/maps/sc_cChInZ_35-18-26_34p-18p-01p-MN40.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
Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
- Double_t intEr, intEphi, intdEz ;
+ Float_t intEr, intEphi, intdEz ;
Double_t r, phi, z ;
Int_t sign;
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)
+ 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)
}
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 *arrayofErA[phiSlices], *arrayofdEzA[phiSlices];
- TMatrixD *arrayofErC[phiSlices], *arrayofdEzC[phiSlices];
-
- TMatrixD *arrayofEroverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
- TMatrixD *arrayofEroverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
+ TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
+ TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
}
// list of points as used during sum up
- Double_t rlist1[rows], zedlist1[columns];// , philist1[phiSlices];
+ 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++ ) {
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)=(%lf,%lf)\n",rlist1[i],zedlist1[j]);
- printf("lookup: (r,z)=(%lf,%lf)\n",coord1(0,ip)*100,coord1(2,ip)*100);
+ 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!");
}
// phi = fgkPhiList[k] ;
// final lookup table
- TMatrixD &erOverEzFinal = *fLookUpErOverEz[k] ;
- TMatrixD &deltaEzFinal = *fLookUpDeltaEz[k] ;
+ TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
+ TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
// calculated and integrated tables - just one phi slice
TMatrixD &erOverEzA = *arrayofEroverEzA[0] ;
gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
// list of points as used during sum up
- Double_t rlist2[rows], philist2[phiSlices], zedlist2[columns];
+ 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++ ) {
// temporary matrices needed for the calculation
- TMatrixD *arrayofErA2[phiSlices], *arrayofEphiA2[phiSlices], *arrayofdEzA2[phiSlices];
- TMatrixD *arrayofErC2[phiSlices], *arrayofEphiC2[phiSlices], *arrayofdEzC2[phiSlices];
+ TMatrixD *arrayofErA2[kNPhiSlices], *arrayofEphiA2[kNPhiSlices], *arrayofdEzA2[kNPhiSlices];
+ TMatrixD *arrayofErC2[kNPhiSlices], *arrayofEphiC2[kNPhiSlices], *arrayofdEzC2[kNPhiSlices];
- TMatrixD *arrayofEroverEzA2[phiSlices], *arrayofEphioverEzA2[phiSlices], *arrayofDeltaEzA2[phiSlices];
- TMatrixD *arrayofEroverEzC2[phiSlices], *arrayofEphioverEzC2[phiSlices], *arrayofDeltaEzC2[phiSlices];
+ TMatrixD *arrayofEroverEzA2[kNPhiSlices], *arrayofEphioverEzA2[kNPhiSlices], *arrayofDeltaEzA2[kNPhiSlices];
+ TMatrixD *arrayofEroverEzC2[kNPhiSlices], *arrayofEphioverEzC2[kNPhiSlices], *arrayofDeltaEzC2[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
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%lf %lf : %e %e\n",r0,phi0,weightA,weightC);
+ // printf("-----\n%f %f : %e %e\n",r0,phi0,weightA,weightC);
// Summing up the vector components according to their weight
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)=(%lf,%lf,%lf)\n",coord2(0,ip)*100,coord2(1,ip),coord2(2,ip)*100);
- printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist2[i],philist2[k],zedlist2[j]);
+ 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!");
}
weightA = GetSpaceChargeDensity(r0*100,phi0R, 0.499, 2); // partial load in r,phi
weightC = GetSpaceChargeDensity(r0*100,phi0R,-0.499, 2); // partial load in r,phi
- // printf("%lf %lf : %e %e\n",r0,phi0R,weightA,weightC);
+ // printf("%f %f : %e %e\n",r0,phi0R,weightA,weightC);
// Summing up the vector components according to their weight
ip = 0;
ipC++; // POC configuration counter
- // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
+ // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
}
Double_t phi = fgkPhiList[k] ;
// final lookup table
- TMatrixD &erOverEzFinal = *fLookUpErOverEz[k] ;
- TMatrixD &ephiOverEzFinal = *fLookUpEphiOverEz[k];
- TMatrixD &deltaEzFinal = *fLookUpDeltaEz[k] ;
+ TMatrixF &erOverEzFinal = *fLookUpErOverEz[k] ;
+ TMatrixF &ephiOverEzFinal = *fLookUpEphiOverEz[k];
+ TMatrixF &deltaEzFinal = *fLookUpDeltaEz[k] ;
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
const Float_t gridSizeZ = fgkTPCZ0/(columns-1); // unit in [cm]
// temporary matrices needed for the calculation
- TMatrixD *arrayofErA[phiSlices], *arrayofEphiA[phiSlices], *arrayofdEzA[phiSlices];
- TMatrixD *arrayofErC[phiSlices], *arrayofEphiC[phiSlices], *arrayofdEzC[phiSlices];
+ TMatrixD *arrayofErA[kNPhiSlices], *arrayofEphiA[kNPhiSlices], *arrayofdEzA[kNPhiSlices];
+ TMatrixD *arrayofErC[kNPhiSlices], *arrayofEphiC[kNPhiSlices], *arrayofdEzC[kNPhiSlices];
- TMatrixD *arrayofEroverEzA[phiSlices], *arrayofEphioverEzA[phiSlices], *arrayofDeltaEzA[phiSlices];
- TMatrixD *arrayofEroverEzC[phiSlices], *arrayofEphioverEzC[phiSlices], *arrayofDeltaEzC[phiSlices];
+ TMatrixD *arrayofEroverEzA[kNPhiSlices], *arrayofEphioverEzA[kNPhiSlices], *arrayofDeltaEzA[kNPhiSlices];
+ TMatrixD *arrayofEroverEzC[kNPhiSlices], *arrayofEphioverEzC[kNPhiSlices], *arrayofDeltaEzC[kNPhiSlices];
for ( Int_t k = 0 ; k < phiSlices ; k++ ) {
}
// list of points as used in the interpolation (during sum up)
- Double_t rlist[rows], zedlist[columns] , philist[phiSlices];
+ 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++ ) {
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)=(%lf,%lf,%lf)\n",coord(0,ip)*100,coord(1,ip),coord(2,ip)*100);
- printf("sum-up: (r,phi,z)=(%lf,%lf,%lf)\n",rlist[i],philist[k],zedlist[j]);
+ 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!");
}
weightA = GetSpaceChargeDensity(r0*100,phi0, z0*100);
weightC = GetSpaceChargeDensity(r0*100,phi0,-z0*100);
- // printf("%lf %lf %lf: %e %e\n",r0,phi0,z0,weightA,weightC);
+ // printf("%f %f %f: %e %e\n",r0,phi0,z0,weightA,weightC);
// Summing up the vector components according to their weight
ip = 0;
} // end phi-POC summation (phiiC)
- // printf("POC: (r,phi,z) = (%lf %lf %lf) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
+ // printf("POC: (r,phi,z) = (%f %f %f) | weight(A,C): %03.1lf %03.1lf\n",r0,phi0,z0, weightA, weightC);
}
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
phi = fgkPhiList[k] ;
- TMatrixD &erOverEz = *fLookUpErOverEz[k] ;
- TMatrixD &ephiOverEz = *fLookUpEphiOverEz[k];
- TMatrixD &deltaEz = *fLookUpDeltaEz[k] ;
+ TMatrixF &erOverEz = *fLookUpErOverEz[k] ;
+ TMatrixF &ephiOverEz = *fLookUpEphiOverEz[k];
+ TMatrixF &deltaEz = *fLookUpDeltaEz[k] ;
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
TMatrixD &scDensityInRPhiC = *fSCdensityInRPhiC;
for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
phi = fgkPhiList[k] ;
- TMatrixD &scDensity = *fSCdensityDistribution[k] ;
+ TMatrixF &scDensity = *fSCdensityDistribution[k] ;
for ( Int_t j = 0 ; j < kNZ ; j++ ) {
z = fgkZList[j] ;
for ( Int_t i = 0 ; i < kNR ; i++ ) {
// Note: input in [cm], output in [C/m^3/e0] !!
//
- if (!fSCdensityDistribution || !fSCdensityInRZ || !fSCdensityInRPhiA || !fSCdensityInRPhiC ) {
- printf("Irgend a schaaas is nuul - argg\n");
- return 0.;
- }
-
while (phi<0) phi += TMath::TwoPi();
while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
sc = 0;
}
- // printf("%lf %lf %lf: %lf\n",r,phi,z,sc);
+ // printf("%f %f %f: %f\n",r,phi,z,sc);
return sc;
}
// 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,190}; // radius in cm
- Double_t secPosC[5] = {5,8,12,15,15}; // sector
- Double_t radialPosC[5] = {210,170,140,120,190}; // radius in cm
+ 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);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff