AliTPCCorrection.cxx ... additional funtionality to deal with 3D problem (e.g. 3D-poisson solver) + calculation of dz distortions
AliTPCInverseCorrection.cxx ... updated desctructor
AliTPCComposedCorrection.cxx ... updated destructor and additional protectections
AliTPCBoundaryVoltError.h ... minor bugfix plus calculation of z distortions due to E field changes
AliTPCBoundaryVoltError.cxx ... minor bugfix plus calculation of z distortions due to E field changes
AliTPCGGVoltError.cxx
AliTPCFCVoltError3D.h ... new class to calculate distortions due to Field cage imperfections (e.g ROD movements)
AliTPCFCVoltError3D.cxx ... new class to calculate distortions due to Field cage imperfections (e.g ROD movements)
AliTPCROCVoltError3D.h ... new class to calculate distortions due to ReadOut chamber misalignments in z
AliTPCROCVoltError3D.cxx ... new class to calculate distortions due to ReadOut chamber misalignments in z
AliTPCSpaceCharge.h ... new class to calculate distortions due to a fixed 2D Space Charge distribution
AliTPCSpaceCharge.cxx ... new class to calculate distortions due to a fixed 2D Space Charge distribution
Calib/maps ... new folder for lookup tables (e.g. ROC z misalignments)
Calib/maps/TPCROCdzSurvey.root ... Look up table for ROC z misalignments (taken from survey in 2006)
TPCbaseLinkDef.h ... included the new classes (see above)
libTPCbase.pkg ... included the new classes (see above)
AliTPCcalibAlign.cxx ... corrected warning
// //
// The class allows "effective Omega Tau" corrections. //
// //
-// NOTE: This class is not capable of calculation z distortions due to //
+// NOTE: This class is capable of calculating z distortions due to //
// drift velocity change in dependence of the electric field!!! //
// //
// date: 01/06/2010 //
AliTPCBoundaryVoltError::AliTPCBoundaryVoltError()
: AliTPCCorrection("BoundaryVoltError","Boundary Voltage Error"),
fC0(0.),fC1(0.),
+ fROCdisplacement(kTRUE),
fInitLookUp(kFALSE)
{
//
// Correction Terms for effective omegaTau; obtained by a laser calibration run
SetOmegaTauT1T2(wt,fT1,fT2);
-
}
// Calculates the correction due e.g. residual voltage errors on the TPC boundaries
//
- if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitBoundaryVoltErrorDistortion() ...");
+ if (!fInitLookUp) {
+ AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
+ InitBoundaryVoltErrorDistortion();
+ }
Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
// note that the poisson solution was linearly mirroed on this grid!
- Double_t intEr, intEphi ;
+ Double_t intEr, intEphi, intdEz ;
Double_t r, phi, z ;
Int_t sign;
// Get the E field integral
Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr );
+ // Get DeltaEz field integral
+ Interpolate2DEdistortion( order, r, z, fLookUpDeltaEz, intdEz );
// Calculate distorted position
if ( r > 0.0 ) {
// Calculate correction in cartesian coordinates
dx[0] = r * TMath::Cos(phi) - x[0];
dx[1] = r * TMath::Sin(phi) - x[1];
- dx[2] = 0.; // z distortion not implemented (1st order distortions)
+ dx[2] = intdEz; // z distortion - (internally scaled with driftvelocity dependency
+ // on the Ez field plus the actual ROC misalignment (if set TRUE)
+
}
TMatrixD voltArrayA(kRows,kColumns), voltArrayC(kRows,kColumns); // boundary vectors
TMatrixD chargeDensity(kRows,kColumns); // dummy charge
TMatrixD arrayErOverEzA(kRows,kColumns), arrayErOverEzC(kRows,kColumns); // solution
+ TMatrixD arrayDeltaEzA(kRows,kColumns), arrayDeltaEzC(kRows,kColumns); // solution
Double_t rList[kRows], zedList[kColumns] ;
// A side boundary vectors
if ( i == 0 ) voltArrayA(i,j) += zed *((fBoundariesA[1]-fBoundariesA[0])/((kColumns-1)*gridSizeZ))
+ fBoundariesA[0] ; // IFC
- if ( j == kColumns-1 ) voltArrayA(i,j) += radius*((fBoundariesA[3]-fBoundariesA[2])/((kRows-1)*gridSizeR+fgkIFCRadius))
+ if ( j == kColumns-1 ) voltArrayA(i,j) += (radius-fgkIFCRadius)*((fBoundariesA[3]-fBoundariesA[2])/((kRows-1)*gridSizeR))
+ fBoundariesA[2] ; // ROC
if ( i == kRows-1 ) voltArrayA(i,j) += zed *((fBoundariesA[4]-fBoundariesA[5])/((kColumns-1)*gridSizeZ))
+ fBoundariesA[5] ; // OFC
- if ( j == 0 ) voltArrayA(i,j) += radius*((fBoundariesA[6]-fBoundariesA[7])/((kRows-1)*gridSizeR+fgkIFCRadius))
+ if ( j == 0 ) voltArrayA(i,j) += (radius-fgkIFCRadius)*((fBoundariesA[6]-fBoundariesA[7])/((kRows-1)*gridSizeR))
+ fBoundariesA[7] ; // CE
-
+
if (symmetry==0) {
// C side boundary vectors
if ( i == 0 ) voltArrayC(i,j) += zed *((fBoundariesC[1]-fBoundariesC[0])/((kColumns-1)*gridSizeZ))
+ fBoundariesC[0] ; // IFC
- if ( j == kColumns-1 ) voltArrayC(i,j) += radius*((fBoundariesC[3]-fBoundariesC[2])/((kRows-1)*gridSizeR+fgkIFCRadius))
+ if ( j == kColumns-1 ) voltArrayC(i,j) += (radius-fgkIFCRadius)*((fBoundariesC[3]-fBoundariesC[2])/((kRows-1)*gridSizeR))
+ fBoundariesC[2] ; // ROC
if ( i == kRows-1 ) voltArrayC(i,j) += zed *((fBoundariesC[4]-fBoundariesC[5])/((kColumns-1)*gridSizeZ))
+ fBoundariesC[5] ; // OFC
- if ( j == 0 ) voltArrayC(i,j) += radius*((fBoundariesC[6]-fBoundariesC[7])/((kRows-1)*gridSizeR+fgkIFCRadius))
+ if ( j == 0 ) voltArrayC(i,j) += (radius-fgkIFCRadius)*((fBoundariesC[6]-fBoundariesC[7])/((kRows-1)*gridSizeR))
+ fBoundariesC[7] ; // CE
-
}
+
}
}
// always solve the problem on the A side
- PoissonRelaxation2D( voltArrayA, chargeDensity, arrayErOverEzA, kRows, kColumns, kIterations ) ;
+ PoissonRelaxation2D( voltArrayA, chargeDensity, arrayErOverEzA, arrayDeltaEzA,
+ kRows, kColumns, kIterations, fROCdisplacement ) ;
if (symmetry!=0) { // A and C side are the same ("anti-symmetric" or "symmetric")
for ( Int_t j = 0 ; j < kColumns ; j++ ) {
for ( Int_t i = 0 ; i < kRows ; i++ ) {
arrayErOverEzC(i,j) = symmetry*arrayErOverEzA(i,j);
+ arrayDeltaEzC(i,j) = -symmetry*arrayDeltaEzA(i,j);
}
}
} else if (symmetry==0) { // A and C side are different - Solve the problem on the C side too
- PoissonRelaxation2D( voltArrayC, chargeDensity, arrayErOverEzC, kRows, kColumns, kIterations ) ;
+ PoissonRelaxation2D( voltArrayC, chargeDensity, arrayErOverEzC, arrayDeltaEzC,
+ kRows, kColumns, kIterations, fROCdisplacement ) ;
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) {
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ arrayDeltaEzC(i,j) = -arrayDeltaEzC(i,j); // negative z coordinate!
+ }
+ }
}
//Interpolate results onto standard grid for Electric Fields
Int_t ilow=0, jlow=0 ;
Double_t z,r;
Float_t saveEr[2] ;
+ Float_t saveEz[2] ;
for ( Int_t i = 0 ; i < kNZ ; ++i ) {
z = TMath::Abs( fgkZList[i] ) ;
for ( Int_t j = 0 ; j < kNR ; ++j ) {
// Linear interpolation !!
r = fgkRList[j] ;
- Search( kRows, rList, r, ilow ) ; // Note switch - R in rows and Z in columns
- Search( kColumns, zedList, z, jlow ) ;
- if ( ilow < 0 ) ilow = 0 ; // check if out of range
- if ( jlow < 0 ) jlow = 0 ;
- if ( ilow + 1 >= kRows - 1 ) ilow = kRows - 2 ;
- if ( jlow + 1 >= kColumns - 1 ) jlow = kColumns - 2 ;
-
- if (fgkZList[i]>0) { // A side solution
- saveEr[0] = arrayErOverEzA(ilow,jlow) +
- (arrayErOverEzA(ilow,jlow+1)-arrayErOverEzA(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
- saveEr[1] = arrayErOverEzA(ilow+1,jlow) +
- (arrayErOverEzA(ilow+1,jlow+1)-arrayErOverEzA(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
- } else if (fgkZList[i]<0) { // C side solution
- saveEr[0] = arrayErOverEzC(ilow,jlow) +
- (arrayErOverEzC(ilow,jlow+1)-arrayErOverEzC(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
- saveEr[1] = arrayErOverEzC(ilow+1,jlow) +
- (arrayErOverEzC(ilow+1,jlow+1)-arrayErOverEzC(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
- } else {
- AliWarning("Field calculation at z=0 (CE) is not allowed!");
- saveEr[0]=0; saveEr[1]=0;
- }
- fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ;
+ Search( kRows, rList, r, ilow ) ; // Note switch - R in rows and Z in columns
+ Search( kColumns, zedList, z, jlow ) ;
+ if ( ilow < 0 ) ilow = 0 ; // check if out of range
+ if ( jlow < 0 ) jlow = 0 ;
+ if ( ilow + 1 >= kRows - 1 ) ilow = kRows - 2 ;
+ if ( jlow + 1 >= kColumns - 1 ) jlow = kColumns - 2 ;
+
+ if (fgkZList[i]>0) { // A side solution
+ saveEr[0] = arrayErOverEzA(ilow,jlow) +
+ (arrayErOverEzA(ilow,jlow+1)-arrayErOverEzA(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEr[1] = arrayErOverEzA(ilow+1,jlow) +
+ (arrayErOverEzA(ilow+1,jlow+1)-arrayErOverEzA(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEz[0] = arrayDeltaEzA(ilow,jlow) +
+ (arrayDeltaEzA(ilow,jlow+1)-arrayDeltaEzA(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEz[1] = arrayDeltaEzA(ilow+1,jlow) +
+ (arrayDeltaEzA(ilow+1,jlow+1)-arrayDeltaEzA(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
+
+ } else if (fgkZList[i]<0) { // C side solution
+ saveEr[0] = arrayErOverEzC(ilow,jlow) +
+ (arrayErOverEzC(ilow,jlow+1)-arrayErOverEzC(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEr[1] = arrayErOverEzC(ilow+1,jlow) +
+ (arrayErOverEzC(ilow+1,jlow+1)-arrayErOverEzC(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEz[0] = arrayDeltaEzC(ilow,jlow) +
+ (arrayDeltaEzC(ilow,jlow+1)-arrayDeltaEzC(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEz[1] = arrayDeltaEzC(ilow+1,jlow) +
+ (arrayDeltaEzC(ilow+1,jlow+1)-arrayDeltaEzC(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
+
+ } else {
+ AliWarning("Field calculation at z=0 (CE) is not allowed!");
+ saveEr[0]=0; saveEr[1]=0;
+ saveEz[0]=0; saveEz[1]=0;
}
+ fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ;
+ fLookUpDeltaEz[i][j] = saveEz[0] + (saveEz[1]-saveEz[0])*(r-rList[ilow])/gridSizeR ;
+ }
}
- /* delete [] saveEr;
- delete [] sVec;
- delete [] rList;
- delete [] zedList;
- */
-
+ voltArrayA.Clear();
+ voltArrayC.Clear();
+ chargeDensity.Clear();
+ arrayErOverEzA.Clear();
+ arrayErOverEzC.Clear();
+ arrayDeltaEzA.Clear();
+ arrayDeltaEzC.Clear();
+
fInitLookUp = kTRUE;
}
printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
}
- if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitBoundaryVoltErrorDistortion() ...");
+ if (!fInitLookUp)
+ AliError("Lookup table was not initialized! You should do InitBoundaryVoltErrorDistortion() ...");
}
fBoundariesA[i]= boundariesA[i];
if (i>5) fBoundariesC[i]= -boundariesA[i]; // setting for the CE is passed to C side
}
-
+ fInitLookUp=kFALSE;
}
void AliTPCBoundaryVoltError::SetBoundariesC(Float_t boundariesC[6]){
//
- // set voltage errors on the TPC boundaries - A side
+ // set voltage errors on the TPC boundaries - C side
//
// Start at IFC at the Central electrode and work clockwise (for C side) through
// IFC, ROC and OFC. The boundary conditions are currently defined to be a linear
for (Int_t i=0; i<6; i++) {
fBoundariesC[i]= boundariesC[i];
}
-
+ fInitLookUp=kFALSE;
}
// setters and getters for conical
void SetBoundariesA(Float_t boundariesA[8]);
void SetBoundariesC(Float_t boundariesC[6]); // CE settings from the A side
+ Float_t GetBoundariesA(Int_t i) const {return fBoundariesA[i]; }
+ Float_t GetBoundariesC(Int_t i) const {return fBoundariesC[i]; }
- Float_t GetBoundariesA(Int_t i) const {return fBoundariesA[i];}
- Float_t GetBoundariesC(Int_t i) const {return fBoundariesC[i];}
+ void SetROCDisplacement(Bool_t flag) { fROCdisplacement = flag; fInitLookUp=kFALSE;}
+ Bool_t GetROCDisplacement() const { return fROCdisplacement; }
void InitBoundaryVoltErrorDistortion();
Float_t fBoundariesA[8]; // Boundary values on the A side (see Setter function)
Float_t fBoundariesC[8]; // Boundary values on the C side (see Setter function)
- Bool_t fInitLookUp; // flag to check it the Look Up table was created
+ Bool_t fROCdisplacement; // flag for ROC displacement (important for z distortions)
+ Bool_t fInitLookUp; // flag to check it the Look Up table was created
Double_t fLookUpErOverEz[kNZ][kNR]; // Array to store electric field integral (int Er/Ez)
+ Double_t fLookUpDeltaEz[kNZ][kNR]; // Array to store electric field integral (int Delta Ez)
// basic numbers for the poisson relaxation //can be set individually in each class
enum {kRows =257}; // grid size in r direction used in the poisson relaxation // ( 2**n + 1 ) eg. 65, 129, 257 etc.
#include <TCollection.h>
#include <TTimeStamp.h>
#include <TIterator.h>
-
+#include "AliLog.h"
#include "AliTPCComposedCorrection.h"
AliTPCComposedCorrection::~AliTPCComposedCorrection() {
//
- // virtual destructor
+ // destructor
//
+ if (!fCorrections) {
+ AliInfo("No Correction-models were set: can not delete them");
+ } else {
+ TIterator *i=fCorrections->MakeIterator();
+ AliTPCCorrection *c;
+ while (0!=(c=dynamic_cast<AliTPCCorrection*>(i->Next()))) {
+ delete c;
+ }
+ delete i;
+ }
}
// This applies all correction and the specified manner (see general
// class description for details).
//
- TIterator *i=fCorrections->MakeIterator();
+
+ if (!fCorrections) {
+ AliInfo("No Corrections-models were set: can not calculate distortions");
+ return;
+ }
+ TIterator *i=fCorrections->MakeIterator();
AliTPCCorrection *c;
switch (fMode) {
case kParallel:
// class descxiption for details).
//
+ if (!fCorrections) {
+ AliInfo("No Corrections-models were set: can not calculate distortions");
+ return;
+ }
TIterator *i=fCorrections->MakeReverseIterator();
AliTPCCorrection *c;
switch (fMode) {
printf("Composed TPC spacepoint correction \"%s\" -- composed of:\n",GetTitle());
TString opt = option; opt.ToLower();
Int_t in=1;
+ if (!fCorrections) {
+ printf(" - composed correction is empty!\n");
+ return;
+ }
TIterator *i=fCorrections->MakeIterator();
AliTPCCorrection *c;
while (0!=(c=dynamic_cast<AliTPCCorrection*>(i->Next()))) {
//
// Initialization funtion (not used at the moment)
//
-
+ if (!fCorrections) {
+ AliInfo("No Correction-models were set");
+ return;
+ }
TIterator *i=fCorrections->MakeIterator();
AliTPCCorrection *c;
while (0!=(c=dynamic_cast<AliTPCCorrection*>(i->Next())))
//
// Update function
//
+ if (!fCorrections) {
+ AliInfo("No Correction-models were set");
+ return;
+ }
TIterator *i=fCorrections->MakeIterator();
AliTPCCorrection *c;
// Note: overwrites previously set values!
//
+ if (!fCorrections) {
+ AliInfo("No Correction-models were set");
+ return;
+ }
+
TIterator *i=fCorrections->MakeIterator();
AliTPCCorrection *c;
while (0!=(c=dynamic_cast<AliTPCCorrection*>(i->Next()))) {
#include <AliCDBStorage.h>
#include <AliCDBId.h>
#include <AliCDBMetaData.h>
-#include "TVectorD.h"
+#include "TVectorD.h"
+#include "AliTPCParamSR.h"
+#include "AliTPCCorrection.h"
+#include "AliLog.h"
-#include "TRandom.h"
#include "AliExternalTrackParam.h"
#include "AliTrackPointArray.h"
#include "TDatabasePDG.h"
#include "AliTPCROC.h"
#include "THnSparse.h"
+#include "AliTPCLaserTrack.h"
+#include "AliESDVertex.h"
+#include "AliVertexerTracks.h"
+#include "TDatabasePDG.h"
+#include "TF1.h"
#include "TRandom.h"
+
+#include "TDatabasePDG.h"
+
#include "AliTPCTransform.h"
#include "AliTPCcalibDB.h"
#include "AliTPCExB.h"
-#include "AliTPCCorrection.h"
-#include "AliTPCRecoParam.h"
-#include "AliExternalTrackParam.h"
-#include "AliTrackPointArray.h"
-#include "TDatabasePDG.h"
-#include "AliTrackerBase.h"
-#include "AliTPCROC.h"
-#include "THnSparse.h"
+#include "AliTPCRecoParam.h"
-#include "AliTPCLaserTrack.h"
-#include "AliESDVertex.h"
-#include "AliVertexerTracks.h"
-#include "TDatabasePDG.h"
-#include "TF1.h"
-
-#include "AliTPCCorrection.h"
-#include "AliLog.h"
ClassImp(AliTPCCorrection)
// FIXME: the following values should come from the database
-const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
-const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
-const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
-const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
-const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
-const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
+const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
+const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
+const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
+const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
+const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
+const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
+const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
+const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
+const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
// FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
-84.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
-90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
-110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
-130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
-150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
-170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
-190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
-210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
-230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
-249.0, 249.5, 250.0, 251.5, 252.0 } ;
-
+ 83.06, 83.5, 84.0, 84.5, 85.0, 85.5, 86.0, 86.5,
+ 87.0, 87.5, 88.0, 88.5, 89.0, 89.5, 90.0, 90.5, 91.0, 92.0,
+ 93.0, 94.0, 95.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
+ 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
+ 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
+ 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
+ 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
+ 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
+ 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
+ 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 241.0, 242.0, 243.0, 244.0,
+ 245.0, 245.5, 246.0, 246.5, 247.0, 247.5, 248.0, 248.5, 249.0, 249.5,
+ 250.0, 250.5, 251.0, 251.5, 252.0, 252.5, 253.0, 253.5, 254.0, 254.5
+} ;
+
const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
-249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
-240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
AliTPCCorrection::AliTPCCorrection()
- : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
+ : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
{
//
// default constructor
//
if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
+
+ // Initialization of interpolation points
+ for (Int_t i = 0; i<kNPhi; i++) {
+ fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);
+ }
+
}
AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
-: TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
+: TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
{
//
// default constructor, that set the name and title
//
if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
+
+ // Initialization of interpolation points
+ for (Int_t i = 0; i<kNPhi; i++) {
+ fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);
+ }
+
}
AliTPCCorrection::~AliTPCCorrection() {
// in respect to position z within the XY plane.
// The histogramm has nx times ny entries.
//
-
+ AliTPCParam* tpcparam = new AliTPCParamSR;
+
TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
nx,-250.,250.,ny,-250.,250.);
Float_t x[3],dx[3];
x[0]=h->GetXaxis()->GetBinCenter(ix);
GetCorrection(x,roc,dx);
Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
- if (90.<=r0 && r0<=250.) {
+ if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
h->SetBinContent(ix,iy,r1-r0);
}
h->SetBinContent(ix,iy,0.);
}
}
+ delete tpcparam;
return h;
}
// The histogramm has nx times ny entries.
//
+ AliTPCParam* tpcparam = new AliTPCParamSR;
+
TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
nx,-250.,250.,ny,-250.,250.);
Float_t x[3],dx[3];
x[0]=h->GetXaxis()->GetBinCenter(ix);
GetCorrection(x,roc,dx);
Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
- if (90.<=r0 && r0<=250.) {
+ if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
Float_t phi0=TMath::ATan2(x[1] ,x[0] );
Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
h->SetBinContent(ix,iy,0.);
}
}
+ delete tpcparam;
+ return h;
+}
+
+TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
+ // in respect to position z within the XY plane.
+ // The histogramm has nx times ny entries.
+ //
+
+ AliTPCParam* tpcparam = new AliTPCParamSR;
+
+ TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
+ nx,-250.,250.,ny,-250.,250.);
+ Float_t x[3],dx[3];
+ x[2]=z;
+ Int_t roc=z>0.?0:18; // FIXME
+ for (Int_t iy=1;iy<=ny;++iy) {
+ x[1]=h->GetYaxis()->GetBinCenter(iy);
+ for (Int_t ix=1;ix<=nx;++ix) {
+ x[0]=h->GetXaxis()->GetBinCenter(ix);
+ GetCorrection(x,roc,dx);
+ Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
+ if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
+ h->SetBinContent(ix,iy,dx[2]);
+ }
+ else
+ h->SetBinContent(ix,iy,0.);
+ }
+ }
+ delete tpcparam;
return h;
}
return h;
}
+TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
+ // in respect to angle phi within the ZR plane.
+ // The histogramm has nx times ny entries.
+ //
+ TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
+ nz,-250.,250.,nr,85.,250.);
+ Float_t x[3],dx[3];
+ for (Int_t ir=1;ir<=nr;++ir) {
+ Float_t radius=h->GetYaxis()->GetBinCenter(ir);
+ x[0]=radius*TMath::Cos(phi);
+ x[1]=radius*TMath::Sin(phi);
+ for (Int_t iz=1;iz<=nz;++iz) {
+ x[2]=h->GetXaxis()->GetBinCenter(iz);
+ Int_t roc=x[2]>0.?0:18; // FIXME
+ GetCorrection(x,roc,dx);
+ h->SetBinContent(iz,ir,dx[2]);
+ }
+ }
+ return h;
+
+}
+
+
TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
const char *xlabel,const char *ylabel,const char *zlabel,
Int_t nbinsx,Double_t xlow,Double_t xup,
return h;
}
-
// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
}
+void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
+ const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
+ Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
+ //
+ // Interpolate table - 3D interpolation
+ //
+
+ Double_t saveEr[10], savedEr[10] ;
+ Double_t saveEphi[10], savedEphi[10] ;
+ Double_t saveEz[10], savedEz[10] ;
+
+ Search( kNZ, fgkZList, z, fILow ) ;
+ Search( kNPhi, fgkPhiList, z, fJLow ) ;
+ Search( kNR, fgkRList, r, fKLow ) ;
+
+ if ( fILow < 0 ) fILow = 0 ; // check if out of range
+ if ( fJLow < 0 ) fJLow = 0 ;
+ if ( fKLow < 0 ) fKLow = 0 ;
+
+ if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
+ if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
+ if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
+
+ for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
+ for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
+ saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
+ saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
+ saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
+ }
+ savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
+ savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
+ savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
+ }
+ erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
+ ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
+ ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
+
+}
+
+Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
+ const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
+ const TMatrixD &array ) {
+ //
+ // Interpolate table (TMatrix format) - 2D interpolation
+ //
+
+ static Int_t jlow = 0, klow = 0 ;
+ Double_t saveArray[10] ;
+
+ Search( nx, xv, x, jlow ) ;
+ Search( ny, yv, y, klow ) ;
+ if ( jlow < 0 ) jlow = 0 ; // check if out of range
+ if ( klow < 0 ) klow = 0 ;
+ if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
+ if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
+
+ for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
+ {
+ Double_t *ajkl = &((TMatrixD&)array)(j,klow);
+ saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
+ }
+
+ return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
+
+}
+
+Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
+ const Int_t nx, const Int_t ny, const Int_t nz,
+ const Double_t xv[], const Double_t yv[], const Double_t zv[],
+ TMatrixD **arrayofArrays ) {
+ //
+ // Interpolate table (TMatrix format) - 3D interpolation
+ //
+
+ static Int_t ilow = 0, jlow = 0, klow = 0 ;
+ Double_t saveArray[10], savedArray[10] ;
+
+ Search( nx, xv, x, ilow ) ;
+ Search( ny, yv, y, jlow ) ;
+ Search( nz, zv, z, klow ) ;
+
+ if ( ilow < 0 ) ilow = 0 ; // check if out of range
+ if ( jlow < 0 ) jlow = 0 ;
+ if ( klow < 0 ) klow = 0 ;
+
+ if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
+ if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
+ if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
+
+ for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
+ {
+ TMatrixD &table = *arrayofArrays[k] ;
+ for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
+ {
+ saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
+ }
+ savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
+ }
+ return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
+
+}
+
+
Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
const Int_t order, const Double_t x ) {
}
-void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
- TMatrixD &arrayErOverEz, const Int_t rows,
- const Int_t columns, const Int_t iterations ) {
+void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
+ TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
+ const Int_t rows, const Int_t columns, const Int_t iterations,
+ const Bool_t rocDisplacement ) {
//
// Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
//
// NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
// So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
//
+ // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
+ //
// Original code by Jim Thomas (STAR TPC Collaboration)
//
iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
+ sumChargeDensity.Clear();
}
// Differentiate V(r) and solve for E(r) using special equations for the first and last rows
// Integrate Er/Ez from Z to zero
for ( Int_t j = 0 ; j < columns ; j++ ) {
for ( Int_t i = 0 ; i < rows ; i++ ) {
+
Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
arrayErOverEz(i,j) = 0.0 ;
+ arrayDeltaEz(i,j) = 0.0 ;
+
for ( Int_t k = j ; k < columns ; k++ ) {
arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
+ arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
if ( index != 4 ) index = 4; else index = 2 ;
}
- if ( index == 4 ) arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
- if ( index == 2 ) arrayErOverEz(i,j) +=
- (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
- -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
- if ( j == columns-2 ) arrayErOverEz(i,j) =
- (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
- +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
- if ( j == columns-1 ) arrayErOverEz(i,j) = 0.0 ;
+ if ( index == 4 ) {
+ arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
+ arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
+ }
+ if ( index == 2 ) {
+ arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
+ -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
+ arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
+ -2.5*(arrayEz(i,columns-1)-ezField));
+ }
+ if ( j == columns-2 ) {
+ arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
+ +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
+ arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
+ +1.5*(arrayEz(i,columns-1)-ezField) ) ;
+ }
+ if ( j == columns-1 ) {
+ arrayErOverEz(i,j) = 0.0 ;
+ arrayDeltaEz(i,j) = 0.0 ;
+ }
}
}
+ // calculate z distortion from the integrated Delta Ez residuals
+ // and include the aquivalence (Volt to cm) of the ROC shift !!
+
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+
+ // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
+ arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
+
+ // ROC Potential in cm aquivalent
+ Double_t dzROCShift = arrayV(i, columns -1)/ezField;
+ if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
+
+ }
+ }
+
+ arrayEr.Clear();
+ arrayEz.Clear();
+
}
+void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
+ TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
+ const Int_t rows, const Int_t columns, const Int_t phislices,
+ const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
+ Bool_t rocDisplacement ) {
+ //
+ // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
+ //
+ // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
+ // The number of rows and COLUMNS can be different.
+ //
+ // ROWS == 2**M + 1
+ // COLUMNS == 2**N + 1
+ // PHISLICES == Arbitrary but greater than 3
+ //
+ // DeltaPhi in Radians
+ //
+ // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
+ // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
+ //
+ // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
+
+ const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
+
+ const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
+ const Float_t gridSizePhi = deltaphi ;
+ const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
+ const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
+ const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
+
+ TMatrixD arrayE(rows,columns) ;
+
+ // Check that the number of rows and columns is suitable for a binary expansion
+ if ( !IsPowerOfTwo((rows-1)) ) {
+ AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
+ return; }
+ if ( !IsPowerOfTwo((columns-1)) ) {
+ AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
+ return; }
+ if ( phislices <= 3 ) {
+ AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
+ return; }
+ if ( phislices > 1000 ) {
+ AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
+ return; }
+
+ // Solve Poisson's equation in cylindrical coordinates by relaxation technique
+ // Allow for different size grid spacing in R and Z directions
+ // Use a binary expansion of the matrix to speed up the solution of the problem
+
+ Int_t loops, mplus, mminus, signplus, signminus ;
+ Int_t ione = (rows-1)/4 ;
+ Int_t jone = (columns-1)/4 ;
+ loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
+ loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
+
+ TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
+
+ for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
+
+ for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
+
+ Float_t tempgridSizeR = gridSizeR * ione ;
+ Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
+ Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
+
+ std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ coef1[i] = 1.0 + tempgridSizeR/(2*radius);
+ coef2[i] = 1.0 - tempgridSizeR/(2*radius);
+ coef3[i] = tempratioPhi/(radius*radius);
+ coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
+ }
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+ TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
+ TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
+ for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
+ if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
+ else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
+ Float_t weight = 0.0 ;
+ Float_t sum = 0.0 ;
+ sumChargeDensity(i,j) = 0.0 ;
+ for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
+ for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
+ if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
+ else
+ weight = 1.0 ;
+ sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
+ sum += weight*radius ;
+ }
+ }
+ sumChargeDensity(i,j) /= sum ;
+ }
+ sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
+ }
+ }
+ }
+
+ for ( Int_t k = 1 ; k <= iterations; k++ ) {
+
+ // over-relaxation index, >= 1 but < 2
+ Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
+ Float_t overRelaxM1 = overRelax - 1.0 ;
+
+ std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
+ std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
+
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ overRelaxcoef4[i] = overRelax * coef4[i] ;
+ overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
+ }
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+
+ mplus = m + 1; signplus = 1 ;
+ mminus = m - 1 ; signminus = 1 ;
+ if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
+ if ( mplus > phislices-1 ) mplus = phislices - 2 ;
+ if ( mminus < 0 ) mminus = 1 ;
+ }
+ else if (symmetry==-1) { // Anti-symmetry in phi
+ if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
+ if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
+ }
+ else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
+ if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
+ if ( mminus < 0 ) mminus = m - 1 + phislices ;
+ }
+ TMatrixD& arrayV = *arrayofArrayV[m] ;
+ TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
+ TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
+ TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
+
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+
+ arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
+ + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
+ - overRelaxcoef5[i] * arrayV(i,j)
+ + coef1[i] * arrayV(i+ione,j)
+ + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
+ + sumChargeDensity(i,j)
+ ) * overRelaxcoef4[i] ;
+ // Note: over-relax the solution at each step. This speeds up the convergance.
+
+ }
+ }
+
+ if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+
+ if ( ione > 1 ) {
+ arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
+ if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
+ }
+ if ( jone > 1 ) {
+ arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
+ if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
+ }
+ if ( ione > 1 && jone > 1 ) {
+ arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
+ if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
+ if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
+ // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
+ }
+ }
+ }
+ }
+
+ }
+ }
+
+ ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
+ jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
+
+ }
+
+ //Differentiate V(r) and solve for E(r) using special equations for the first and last row
+ //Integrate E(r)/E(z) from point of origin to pad plane
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+ TMatrixD& arrayV = *arrayofArrayV[m] ;
+ TMatrixD& eroverEz = *arrayofEroverEz[m] ;
+
+ for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+
+ // Differentiate in R
+ for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
+ arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
+ arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
+ // Integrate over Z
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ eroverEz(i,j) = 0.0 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+
+ eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
+ if ( index == 2 ) eroverEz(i,j) +=
+ (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-2 ) eroverEz(i,j) =
+ (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
+ }
+ }
+ // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
+ // eroverEz.Draw("surf") ; } // JT test
+ }
+
+ //Differentiate V(r) and solve for E(phi)
+ //Integrate E(phi)/E(z) from point of origin to pad plane
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+
+ mplus = m + 1; signplus = 1 ;
+ mminus = m - 1 ; signminus = 1 ;
+ if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
+ if ( mplus > phislices-1 ) mplus = phislices - 2 ;
+ if ( mminus < 0 ) mminus = 1 ;
+ }
+ else if (symmetry==-1) { // Anti-symmetry in phi
+ if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
+ if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
+ }
+ else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
+ if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
+ if ( mminus < 0 ) mminus = m - 1 + phislices ;
+ }
+ TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
+ TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
+ TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
+ for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+ // Differentiate in Phi
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
+ }
+ // Integrate over Z
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ ePhioverEz(i,j) = 0.0 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+
+ ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
+ if ( index == 2 ) ePhioverEz(i,j) +=
+ (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-2 ) ePhioverEz(i,j) =
+ (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
+ }
+ }
+ // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
+ // arrayE.Draw("surf") ; } // JT test
+ }
+
+
+ // Differentiate V(r) and solve for E(z) using special equations for the first and last row
+ // Integrate (E(z)-Ezstd) from point of origin to pad plane
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+ TMatrixD& arrayV = *arrayofArrayV[m] ;
+ TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
+
+ // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
+ for ( Int_t i = 0 ; i < rows ; i++) {
+ for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
+ arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
+ arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
+ }
+
+ for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+ // Integrate over Z
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ deltaEz(i,j) = 0.0 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+ deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
+ if ( index == 2 ) deltaEz(i,j) +=
+ (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
+ if ( j == columns-2 ) deltaEz(i,j) =
+ (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
+ if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
+ }
+ }
+ // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
+ // eroverEz.Draw("surf") ; } // JT test
+
+ // calculate z distortion from the integrated Delta Ez residuals
+ // and include the aquivalence (Volt to cm) of the ROC shift !!
+
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+
+ // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
+ deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
+
+ // ROC Potential in cm aquivalent
+ Double_t dzROCShift = arrayV(i, columns -1)/ezField;
+ if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
+
+ }
+ }
+
+ } // end loop over phi
+
+
+
+ for ( Int_t k = 0 ; k < phislices ; k++ )
+ {
+ arrayofSumChargeDensities[k]->Delete() ;
+ }
+
+
+
+ arrayE.Clear();
+}
Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
// corrArray - array with partial corrections
// step - skipe entries - if 1 all entries processed - it is slow
// debug 0 if debug on also space points dumped - it is slow
+
const Double_t kMaxSnp = 0.85;
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
// const Double_t kB2C=-0.299792458e-3;
void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
+ //
+ // Fast method to simulate the influence of the given distortion on the vertex reconstruction
+ //
- AliVertexerTracks *vertexer = new AliVertexerTracks(5);// 5kGaus
+ AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+ if (!magF) AliError("Magneticd field - not initialized");
+ Double_t bz = magF->SolenoidField(); //field in kGauss
+ printf("bz: %lf\n",bz);
+ AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
- TObjArray ATrk; // Original Track array of Aside
- TObjArray dATrk; // Distorted Track array of A side
- UShort_t *AId = new UShort_t[nTracks]; // A side Track ID
- TObjArray CTrk;
- TObjArray dCTrk;
- UShort_t *CId = new UShort_t [nTracks];
- Int_t ID=0;
+ TObjArray aTrk; // Original Track array of Aside
+ TObjArray daTrk; // Distorted Track array of A side
+ UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
+ TObjArray cTrk;
+ TObjArray dcTrk;
+ UShort_t *cId = new UShort_t [nTracks];
+ Int_t id=0;
Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
fpt.SetParameters(7.24,0.120);
for(Int_t nt=0; nt<nTracks; nt++){
Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
- Double_t pt = fpt.GetRandom();// momentum for f1
+ Double_t pt = fpt.GetRandom(); // momentum for f1
+ // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
Short_t sign=1;
if(gRandom->Rndm() < 0.5){
sign =1;
"\n";
if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
if (td){
- dATrk.AddLast(td);
- ATrk.AddLast(t);
- Int_t nn=ATrk.GetEntriesFast();
- AId[nn]=ID;
+ daTrk.AddLast(td);
+ aTrk.AddLast(t);
+ Int_t nn=aTrk.GetEntriesFast();
+ aId[nn]=id;
}
}else if(( eta<-0.07 )&&( eta>-etaCuts )){
if (td){
- dCTrk.AddLast(td);
- CTrk.AddLast(t);
- Int_t nn=CTrk.GetEntriesFast();
- CId[nn]=ID;
+ dcTrk.AddLast(td);
+ cTrk.AddLast(t);
+ Int_t nn=cTrk.GetEntriesFast();
+ cId[nn]=id;
}
}
- ID++;
+ id++;
}// end of track loop
vertexer->SetTPCMode();
vertexer->SetConstraintOff();
- aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dATrk,AId));
- avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&ATrk,AId));
- cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dCTrk,CId));
- cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&CTrk,CId));
+ aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
+ avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
+ cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
+ cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
if (pcstream) (*pcstream)<<"vertex"<<
"x="<<orgVertex[0]<<
"y="<<orgVertex[1]<<
"avO.="<<&avOrg<< // original vertex A side
"cvO.="<<&cvOrg<<
"\n";
- delete []AId;
- delete []CId;
+ delete []aId;
+ delete []cId;
}
void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
TH2F* CreateHistoDRinXY (Float_t z=10.,Int_t nx=100,Int_t ny=100);
TH2F* CreateHistoDRPhiinXY(Float_t z=10.,Int_t nx=100,Int_t nphi=100);
+ TH2F* CreateHistoDZinXY (Float_t z=10.,Int_t nx=100,Int_t ny=100);
+
TH2F* CreateHistoDRinZR (Float_t phi=0.,Int_t nZ=100,Int_t nR=100);
TH2F* CreateHistoDRPhiinZR(Float_t phi=0.,Int_t nZ=100,Int_t nR=100);
+ TH2F* CreateHistoDZinZR (Float_t phi=0.,Int_t nZ=100,Int_t nR=100);
+
+
TTree* CreateDistortionTree(Double_t step=5);
static void MakeDistortionMap(THnSparse * his0, TTreeSRedirector *pcstream, const char* hname, Int_t run);
// normally called directly in the correction classes which inherit from this class
Int_t nbinsx,Double_t xlow,Double_t xup,
Int_t nbinsy,Double_t ylow,Double_t yup);
- static const Double_t fgkTPCZ0; // nominal gating grid position
+ static const Double_t fgkTPCZ0; // nominal gating grid position
static const Double_t fgkIFCRadius; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
static const Double_t fgkOFCRadius; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
static const Double_t fgkZOffSet; // Offset from CE: calculate all distortions closer to CE as if at this point
static const Double_t fgkCathodeV; // Cathode Voltage (volts)
static const Double_t fgkGG; // Gating Grid voltage (volts)
+ static const Double_t fgkdvdE; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
+ static const Double_t fgkEM; // charge/mass in [C/kg]
+ static const Double_t fgke0; // vacuum permittivity [A·s/(V·m)]
- enum {kNR= 92}; // Number of R points in the table for interpolating distortion data
+
+ enum {kNR= 118}; // Number of R points in the table for interpolating distortion data
+ enum {kNPhi= 18*10+1}; // Number of Phi points in the table for interpolating distortion data ( plus one extra for 360 == 0 )
enum {kNZ= 270}; // Number of Z points in the table for interpolating distortion data
+
static const Double_t fgkRList[kNR]; // points in the radial direction (for the lookup table)
+ Double_t fgkPhiList[kNPhi]; // points in the phi direction (for the lookup table)
static const Double_t fgkZList[kNZ]; // points in the z direction (for the lookup table)
// Simple Interpolation functions: e.g. with tricubic interpolation (not yet in TH3)
- Int_t fJLow; // variable to help in the interpolation
- Int_t fKLow; // variable to help in the interpolation
+ Int_t fILow, fJLow, fKLow; // variable to help in the interpolation
void Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
const Double_t er[kNZ][kNR], Double_t &erValue );
+ void Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
+ const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
+ Double_t &erValue, Double_t &ephiValue, Double_t &ezValue);
+ Double_t Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
+ const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
+ const TMatrixD &array );
+ Double_t Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
+ const Int_t nx, const Int_t ny, const Int_t nz,
+ const Double_t xv[], const Double_t yv[], const Double_t zv[],
+ TMatrixD **arrayofArrays );
Double_t Interpolate( const Double_t xArray[], const Double_t yArray[],
const Int_t order, const Double_t x );
void Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low );
virtual Int_t IsPowerOfTwo ( Int_t i ) const ;
-
- // Algorithms to solve the laplace or possion equation
- void PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity, TMatrixD &arrayErOverEz, const Int_t rows, const Int_t columns, const Int_t iterations );
+
+
+ // Algorithms to solve the laplace or poisson equation
+ void PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
+ TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
+ const Int_t rows, const Int_t columns, const Int_t iterations,
+ const Bool_t rocDisplacement = kTRUE);
+
+ void PoissonRelaxation3D( TMatrixD **arrayofArrayV, TMatrixD **arrayofChargeDensities,
+ TMatrixD **arrayofEroverEz, TMatrixD **arrayofEPhioverEz, TMatrixD **arrayofEz,
+ const Int_t rows, const Int_t columns, const Int_t phislices,
+ const Float_t deltaphi, const Int_t iterations, const Int_t summetry,
+ const Bool_t rocDisplacement = kTRUE);
+
protected:
Double_t fT1; // tensor term of wt - T1
Double_t fT2; // tensor term of wt - T2
static TObjArray *fgVisualCorrection; // array of orrection for visualization
private:
- ClassDef(AliTPCCorrection,2);
+ AliTPCCorrection(const AliTPCCorrection &); // not implemented
+ AliTPCCorrection &operator=(const AliTPCCorrection &); // not implemented
+
+ ClassDef(AliTPCCorrection,3);
};
#endif
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCFCVoltError3D class //
+// The class calculates the space point distortions due to residual voltage //
+// errors on the Field Cage boundaries (rods) of the TPC in 3D. //
+// //
+// The class allows "effective Omega Tau" corrections. //
+// //
+// NOTE: This class is capable of calculating z distortions due to //
+// drift velocity changes in dependence of the electric field!!! //
+// //
+// date: 08/08/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+// //
+// Example usage : //
+// AliTPCFCVoltError3D fcerror; //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "AliMagF.h"
+#include "TGeoGlobalMagField.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
+#include "AliLog.h"
+#include "TMatrixD.h"
+
+#include "TMath.h"
+#include "AliTPCROC.h"
+#include "AliTPCFCVoltError3D.h"
+
+ClassImp(AliTPCFCVoltError3D)
+
+AliTPCFCVoltError3D::AliTPCFCVoltError3D()
+ : AliTPCCorrection("FieldCageVoltErrors","FieldCage (Rods) Voltage Errors"),
+ fC0(0.),fC1(0.),
+ fInitLookUp(kFALSE)
+{
+ //
+ // default constructor
+ //
+
+ // flags for filled 'basic lookup tables'
+ for (Int_t i=0; i<5; i++){
+ fInitLookUpBasic[i]= kFALSE;
+ }
+
+ // Boundary settings
+ for (Int_t i=0; i<36; i++){
+ fRodVoltShiftA[i] = 0;
+ fRodVoltShiftC[i] = 0;
+ }
+ for (Int_t i=0; i<2; i++){
+ fRotatedClipVoltA[i] = 0;
+ fRotatedClipVoltC[i] = 0;
+ }
+ //
+ for (Int_t i=0; i<18; i++){
+ fOFCRodShiftA[i] = 0;
+ fOFCRodShiftC[i] = 0;
+ }
+
+ // Array which will contain the solution according to the setted boundary conditions
+ // it represents a sum up of the 4 basic look up tables (created individually)
+ // see InitFCVoltError3D() 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);
+ }
+
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ fLookUpBasic1ErOverEz[k] = 0;
+ fLookUpBasic1EphiOverEz[k] = 0;
+ fLookUpBasic1DeltaEz[k] = 0;
+
+ fLookUpBasic2ErOverEz[k] = 0;
+ fLookUpBasic2EphiOverEz[k] = 0;
+ fLookUpBasic2DeltaEz[k] = 0;
+
+ fLookUpBasic3ErOverEz[k] = 0;
+ fLookUpBasic3EphiOverEz[k] = 0;
+ fLookUpBasic3DeltaEz[k] = 0;
+
+ fLookUpBasic4ErOverEz[k] = 0;
+ fLookUpBasic4EphiOverEz[k] = 0;
+ fLookUpBasic4DeltaEz[k] = 0;
+
+ fLookUpBasic5ErOverEz[k] = 0;
+ fLookUpBasic5EphiOverEz[k] = 0;
+ fLookUpBasic5DeltaEz[k] = 0;
+ }
+
+}
+
+AliTPCFCVoltError3D::~AliTPCFCVoltError3D() {
+ //
+ // destructor
+ //
+
+ for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
+ delete fLookUpErOverEz[k];
+ delete fLookUpEphiOverEz[k];
+ delete fLookUpDeltaEz[k];
+ }
+
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ delete fLookUpBasic1ErOverEz[k]; // does nothing if pointer is zero!
+ delete fLookUpBasic1EphiOverEz[k];
+ delete fLookUpBasic1DeltaEz[k];
+
+ delete fLookUpBasic2ErOverEz[k]; // does nothing if pointer is zero!
+ delete fLookUpBasic2EphiOverEz[k];
+ delete fLookUpBasic2DeltaEz[k];
+
+ delete fLookUpBasic3ErOverEz[k]; // does nothing if pointer is zero!
+ delete fLookUpBasic3EphiOverEz[k];
+ delete fLookUpBasic3DeltaEz[k];
+
+ delete fLookUpBasic4ErOverEz[k]; // does nothing if pointer is zero!
+ delete fLookUpBasic4EphiOverEz[k];
+ delete fLookUpBasic4DeltaEz[k];
+
+ delete fLookUpBasic5ErOverEz[k]; // does nothing if pointer is zero!
+ delete fLookUpBasic5EphiOverEz[k];
+ delete fLookUpBasic5DeltaEz[k];
+ }
+}
+
+void AliTPCFCVoltError3D::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);
+
+ InitFCVoltError3D();
+}
+
+void AliTPCFCVoltError3D::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);
+
+
+}
+
+
+
+void AliTPCFCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
+ //
+ // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
+ //
+
+ if (!fInitLookUp) {
+ AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
+ InitFCVoltError3D();
+ return;
+ }
+
+ Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
+ // note that the poisson solution was linearly mirroed on this grid!
+ Double_t intEr, intEphi, intDeltaEz;
+ 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 );
+ intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
+ fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
+
+ // printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
+
+ // Calculate distorted position
+ if ( r > 0.0 ) {
+ phi = phi + ( fC0*intEphi - fC1*intEr ) / r;
+ r = r + ( fC0*intEr + fC1*intEphi );
+ }
+
+ // Calculate correction in cartesian coordinates
+ dx[0] = r * TMath::Cos(phi) - x[0];
+ dx[1] = r * TMath::Sin(phi) - x[1];
+ dx[2] = intDeltaEz; // z distortion - (internally scaled with driftvelocity dependency
+ // on the Ez field plus the actual ROC misalignment (if set TRUE)
+
+}
+
+void AliTPCFCVoltError3D::InitFCVoltError3D() {
+ //
+ // Initialization of the Lookup table which contains the solutions of the
+ // Dirichlet boundary problem
+ // Calculation of the single 3D-Poisson solver is done just if needed
+ // (see basic lookup tables in header file)
+ //
+
+ const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
+ const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
+ const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
+ const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
+
+ // temporary arrays to create the boundary conditions
+ TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
+ arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
+ }
+ Double_t innerList[kPhiSlices], outerList[kPhiSlices] ;
+
+ // list of point as used in the poisson relation and the interpolation (during sum up)
+ Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ philist[k] = gridSizePhi * k;
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ rlist[i] = fgkIFCRadius + i*gridSizeR ;
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
+ zedlist[j] = j * gridSizeZ ;
+ }
+ }
+ }
+
+ // ==========================================================================
+ // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
+ // Allow for different size grid spacing in R and Z directions
+
+ const Int_t symmetry[5] = {1,1,-1,-1,1}; // shifted rod (2x), rotated clip (2x), only rod shift on OFC (1x)
+
+ // check which lookup tables are needed ---------------------------------
+
+ Bool_t needTable[5] ={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
+
+ // Shifted rods & strips
+ for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
+ if (fRodVoltShiftA[rod]!=0) needTable[0]=kTRUE;
+ if (fRodVoltShiftC[rod]!=0) needTable[0]=kTRUE;
+ }
+ for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
+ if (fRodVoltShiftA[rod]!=0) needTable[1]=kTRUE;
+ if (fRodVoltShiftC[rod]!=0) needTable[1]=kTRUE;
+ }
+ // Rotated clips
+ if (fRotatedClipVoltA[0]!=0 || fRotatedClipVoltC[0]!=0) needTable[2]=kTRUE;
+ if (fRotatedClipVoltA[1]!=0 || fRotatedClipVoltC[1]!=0) needTable[3]=kTRUE;
+
+ // shifted Copper rods on OFC
+ for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
+ if (fOFCRodShiftA[rod]!=0) needTable[4]=kTRUE;
+ if (fOFCRodShiftC[rod]!=0) needTable[4]=kTRUE;
+ }
+
+
+ // reserve the arrays for the basic lookup tables ----------------------
+ if (needTable[0] && !fInitLookUpBasic[0]) {
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
+ fLookUpBasic1ErOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic1EphiOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic1DeltaEz[k] = new TMatrixD(kRows,kColumns);
+ // will be deleted in destructor
+ }
+ }
+ if (needTable[1] && !fInitLookUpBasic[1]) {
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
+ fLookUpBasic2ErOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic2EphiOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic2DeltaEz[k] = new TMatrixD(kRows,kColumns);
+ // will be deleted in destructor
+ }
+ }
+ if (needTable[2] && !fInitLookUpBasic[2]) {
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
+ fLookUpBasic3ErOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic3EphiOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic3DeltaEz[k] = new TMatrixD(kRows,kColumns);
+ // will be deleted in destructor
+ }
+ }
+ if (needTable[3] && !fInitLookUpBasic[3]) {
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
+ fLookUpBasic4ErOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic4EphiOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic4DeltaEz[k] = new TMatrixD(kRows,kColumns);
+ // will be deleted in destructor
+ }
+ }
+ if (needTable[4] && !fInitLookUpBasic[4]) {
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) { // Possibly make an extra table to be used for 0 == 360
+ fLookUpBasic5ErOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic5EphiOverEz[k] = new TMatrixD(kRows,kColumns);
+ fLookUpBasic5DeltaEz[k] = new TMatrixD(kRows,kColumns);
+ // will be deleted in destructor
+ }
+ }
+
+ // Set bondaries and solve Poisson's equation --------------------------
+
+ for (Int_t look=0; look<5; look++) {
+
+ Float_t inner18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
+ Float_t outer18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
+
+ if (needTable[look] && !fInitLookUpBasic[look]) {
+
+ // Specify which rod is the reference; other distortions calculated by summing over multiple rotations of refrence
+ // Note: the interpolation later on depends on it!! Do not change if not really needed!
+ if (look==0) {
+ AliInfo(Form("IFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
+ inner18[0] = 1;
+ } else if (look==1) {
+ AliInfo(Form("OFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
+ outer18[0] = 1;
+ } else if (look==2) {
+ AliInfo(Form("IFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
+ inner18[0] = 1;
+ } else if (look==3) {
+ AliInfo(Form("OFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
+ outer18[0] = 1;
+ } else if (look==4) {
+ AliInfo(Form("OFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
+ outer18[0] = 1;
+ }
+ // Build potentials on boundary stripes for a rotated clip (SYMMETRY==-1) or a shifted rod (SYMMETRY==1)
+ if (look!=4) {
+ // in these cases, the strips follow the actual rod shift (linear interpolation between the rods)
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ Int_t slices = kPhiSlicesPerSector ;
+ Int_t ipoint = k/slices ;
+ innerList[k] = (((ipoint+1)*slices-k)*-1*inner18[ipoint]+(k-ipoint*slices)*inner18[ipoint+1])/slices ;
+ outerList[k] = (((ipoint+1)*slices-k)*-1*outer18[ipoint]+(k-ipoint*slices)*outer18[ipoint+1])/slices ;
+ if ( (k%slices) == 0 && symmetry[look] == -1 ) { innerList[k] = 0.0 ; outerList[k] = 0.0 ; }
+ // above, force Zero if Anti-Sym
+ }
+ } else {
+ // in this case, we assume that the strips stay at the correct position, but the copper plated OFC-rods move
+ // the distortion is then much more localized around the rod (indicated by real data)
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
+ innerList[k] = outerList[k] = 0;
+ outerList[0] = outer18[0]; // point at rod
+ outerList[1] = outer18[0]/4; // point close to rod (educated-`guessed` scaling)
+ }
+
+ // Fill arrays with initial conditions. V on the boundary and Charge in the volume.
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ TMatrixD &arrayV = *arrayofArrayV[k] ;
+ TMatrixD &charge = *arrayofCharge[k] ;
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
+ arrayV(i,j) = 0.0 ;
+ charge(i,j) = 0.0 ;
+ if ( i == 0 ) arrayV(i,j) = innerList[k] ;
+ if ( i == kRows-1 ) arrayV(i,j) = outerList[k] ;
+ }
+ }
+ // no charge in the volume
+ for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
+ for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
+ charge(i,j) = 0.0 ;
+ }
+ }
+ }
+
+ // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
+ // Allow for different size grid spacing in R and Z directions
+ if (look==0) {
+ PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
+ fLookUpBasic1ErOverEz, fLookUpBasic1EphiOverEz, fLookUpBasic1DeltaEz,
+ kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[0]) ;
+ AliInfo("IFC - ROD&Strip shift : done ");
+ } else if (look==1) {
+ PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
+ fLookUpBasic2ErOverEz, fLookUpBasic2EphiOverEz, fLookUpBasic2DeltaEz,
+ kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[1]) ;
+
+ AliInfo("OFC - ROD&Strip shift : done ");
+ } else if (look==2) {
+ PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
+ fLookUpBasic3ErOverEz, fLookUpBasic3EphiOverEz, fLookUpBasic3DeltaEz,
+ kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[2]) ;
+ AliInfo("IFC - Clip rot. : done ");
+ } else if (look==3) {
+ PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
+ fLookUpBasic4ErOverEz, fLookUpBasic4EphiOverEz, fLookUpBasic4DeltaEz,
+ kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[3]) ;
+ AliInfo("OFC - Clip rot. : done ");
+ } else if (look==4) {
+ PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
+ fLookUpBasic5ErOverEz, fLookUpBasic5EphiOverEz, fLookUpBasic5DeltaEz,
+ kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[4]) ;
+ AliInfo("OFC - CopperRod shift : done ");
+ }
+
+ fInitLookUpBasic[look] = kTRUE;
+ }
+ }
+
+
+ // =============================================================================
+ // Create the final lookup table with corresponding ROD shifts or clip rotations
+
+ Float_t erIntegralSum = 0.0 ;
+ Float_t ephiIntegralSum = 0.0 ;
+ Float_t ezIntegralSum = 0.0 ;
+
+ Double_t phiPrime = 0. ;
+ Double_t erIntegral = 0. ;
+ Double_t ephiIntegral = 0. ;
+ Double_t ezIntegral = 0. ;
+
+
+ AliInfo("Calculation of combined Look-up Table");
+
+ // Interpolate and sum the Basic lookup tables onto the standard grid
+ 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 i = 0 ; i < kNZ ; i++ ) {
+ z = TMath::Abs(fgkZList[i]) ; // Symmetric solution in Z that depends only on ABS(Z)
+ for ( Int_t j = 0 ; j < kNR ; j++ ) {
+ r = fgkRList[j] ;
+ // Interpolate basicLookup tables; once for each rod, then sum the results
+
+ erIntegralSum = 0.0 ;
+ ephiIntegralSum = 0.0 ;
+ ezIntegralSum = 0.0 ;
+
+ // SHIFTED RODS =========================================================
+
+ // IFC ROD SHIFTS +++++++++++++++++++++++++++++
+ for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
+
+ erIntegral = ephiIntegral = ezIntegral = 0.0 ;
+
+ if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
+ if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
+
+ // Rotate to a position where we have maps and prepare to sum
+ phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
+
+ if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
+
+ if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic1ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic1DeltaEz );
+
+ } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
+
+ phiPrime = TMath::TwoPi() - phiPrime ;
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic1ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic1DeltaEz );
+
+ // Flip symmetry of phi integral for symmetric boundary conditions
+ if ( symmetry[0] == 1 ) ephiIntegral *= -1 ;
+ // Flip symmetry of r integral if anti-symmetric boundary conditions
+ if ( symmetry[0] == -1 ) erIntegral *= -1 ;
+
+ }
+
+ if ( fgkZList[i] > 0 ) {
+ erIntegralSum += fRodVoltShiftA[rod]*erIntegral ;
+ ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
+ ezIntegralSum += fRodVoltShiftA[rod]*ezIntegral;
+ } else if ( fgkZList[i] < 0 ) {
+ erIntegralSum += fRodVoltShiftC[rod]*erIntegral ;
+ ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
+ ezIntegralSum -= fRodVoltShiftC[rod]*ezIntegral;
+ }
+ }
+
+ // OFC ROD SHIFTS +++++++++++++++++++++++++++++
+ for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
+
+ erIntegral = ephiIntegral = ezIntegral = 0.0 ;
+
+ if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
+ if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
+
+ // Rotate to a position where we have maps and prepare to sum
+ phiPrime = phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;
+
+
+ if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
+
+ if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic2ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic2DeltaEz );
+
+ } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
+
+ phiPrime = TMath::TwoPi() - phiPrime ;
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic2ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic2DeltaEz );
+
+ // Flip symmetry of phi integral for symmetric boundary conditions
+ if ( symmetry[1] == 1 ) ephiIntegral *= -1 ;
+ // Flip symmetry of r integral if anti-symmetric boundary conditions
+ if ( symmetry[1] == -1 ) erIntegral *= -1 ;
+
+ }
+
+ if ( fgkZList[i] > 0 ) {
+ erIntegralSum += fRodVoltShiftA[rod]*erIntegral ;
+ ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
+ ezIntegralSum += fRodVoltShiftA[rod]*ezIntegral;
+ } else if ( fgkZList[i] < 0 ) {
+ erIntegralSum += fRodVoltShiftC[rod]*erIntegral ;
+ ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
+ ezIntegralSum -= fRodVoltShiftC[rod]*ezIntegral;
+ }
+
+ } // rod loop - shited ROD
+
+
+ // Rotated clips =========================================================
+
+ Int_t rodIFC = 11; // resistor rod positions, IFC
+ Int_t rodOFC = 3; // resistor rod positions, OFC
+ // just one rod on IFC and OFC
+
+ // IFC ROTATED CLIP +++++++++++++++++++++++++++++
+ for ( Int_t rod = rodIFC ; rod < rodIFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
+
+ erIntegral = ephiIntegral = ezIntegral = 0.0 ;
+ if ( fRotatedClipVoltA[0] == 0 && fgkZList[i] > 0) continue ;
+ if ( fRotatedClipVoltC[0] == 0 && fgkZList[i] < 0) continue ;
+
+ // Rotate to a position where we have maps and prepare to sum
+ phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
+
+ if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
+
+ if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic3ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic3DeltaEz );
+
+ } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
+
+ phiPrime = TMath::TwoPi() - phiPrime ;
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic3ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic3DeltaEz );
+
+ // Flip symmetry of phi integral for symmetric boundary conditions
+ if ( symmetry[2] == 1 ) ephiIntegral *= -1 ;
+ // Flip symmetry of r integral if anti-symmetric boundary conditions
+ if ( symmetry[2] == -1 ) erIntegral *= -1 ;
+
+ }
+
+ if ( fgkZList[i] > 0 ) {
+ erIntegralSum += fRotatedClipVoltA[0]*erIntegral ;
+ ephiIntegralSum += fRotatedClipVoltA[0]*ephiIntegral ;
+ ezIntegralSum += fRotatedClipVoltA[0]*ezIntegral;
+ } else if ( fgkZList[i] < 0 ) {
+ erIntegralSum += fRotatedClipVoltC[0]*erIntegral ;
+ ephiIntegralSum += fRotatedClipVoltC[0]*ephiIntegral ;
+ ezIntegralSum -= fRotatedClipVoltC[0]*ezIntegral;
+ }
+ }
+
+ // OFC: ROTATED CLIP +++++++++++++++++++++++++++++
+ for ( Int_t rod = rodOFC ; rod < rodOFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
+
+ erIntegral = ephiIntegral = ezIntegral = 0.0 ;
+
+ if ( fRotatedClipVoltA[1] == 0 && fgkZList[i] > 0) continue ;
+ if ( fRotatedClipVoltC[1] == 0 && fgkZList[i] < 0) continue ;
+
+ // Rotate to a position where we have maps and prepare to sum
+ phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
+
+
+ if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
+
+ if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic4ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic4DeltaEz );
+
+ } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
+
+ phiPrime = TMath::TwoPi() - phiPrime ;
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic4ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic4DeltaEz );
+
+ // Flip symmetry of phi integral for symmetric boundary conditions
+ if ( symmetry[3] == 1 ) ephiIntegral *= -1 ;
+ // Flip symmetry of r integral if anti-symmetric boundary conditions
+ if ( symmetry[3] == -1 ) erIntegral *= -1 ;
+
+ }
+
+ if ( fgkZList[i] > 0 ) {
+ erIntegralSum += fRotatedClipVoltA[1]*erIntegral ;
+ ephiIntegralSum += fRotatedClipVoltA[1]*ephiIntegral ;
+ ezIntegralSum += fRotatedClipVoltA[1]*ezIntegral;
+ } else if ( fgkZList[i] < 0 ) {
+ erIntegralSum += fRotatedClipVoltC[1]*erIntegral ;
+ ephiIntegralSum += fRotatedClipVoltC[1]*ephiIntegral ;
+ ezIntegralSum -= fRotatedClipVoltC[1]*ezIntegral;
+ }
+ }
+
+ // OFC Cooper rod shift +++++++++++++++++++++++++++++
+ for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
+
+ erIntegral = ephiIntegral = ezIntegral = 0.0 ;
+
+ if ( fOFCRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
+ if ( fOFCRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
+
+ // Rotate to a position where we have maps and prepare to sum
+ phiPrime = phi - rod*kPhiSlicesPerSector*gridSizePhi ;
+
+ if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ; // Stay in range from 0 to TwoPi
+
+ if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic5ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic5DeltaEz );
+
+ } else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
+
+ phiPrime = TMath::TwoPi() - phiPrime ;
+
+ erIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic5ErOverEz );
+ ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
+ ezIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, fLookUpBasic5DeltaEz );
+
+ // Flip symmetry of phi integral for symmetric boundary conditions
+ if ( symmetry[0] == 1 ) ephiIntegral *= -1 ;
+ // Flip symmetry of r integral if anti-symmetric boundary conditions
+ if ( symmetry[0] == -1 ) erIntegral *= -1 ;
+
+ }
+
+ if ( fgkZList[i] > 0 ) {
+ erIntegralSum += fOFCRodShiftA[rod]*erIntegral ;
+ ephiIntegralSum += fOFCRodShiftA[rod]*ephiIntegral ;
+ ezIntegralSum += fOFCRodShiftA[rod]*ezIntegral;
+ } else if ( fgkZList[i] < 0 ) {
+ erIntegralSum += fOFCRodShiftC[rod]*erIntegral ;
+ ephiIntegralSum += fOFCRodShiftC[rod]*ephiIntegral ;
+ ezIntegralSum -= fOFCRodShiftC[rod]*ezIntegral;
+ }
+ }
+
+ // put the sum into the final lookup table
+ erOverEz(j,i) = erIntegralSum;
+ ephiOverEz(j,i) = ephiIntegralSum;
+ deltaEz(j,i) = ezIntegralSum;
+
+ // if (j==1) printf("%lf %lf %lf: %lf %lf %lf\n",r, phi, z, erIntegralSum,ephiIntegralSum,ezIntegralSum);
+
+ } // endo r loop
+ } // end of z loop
+ } // end of phi loop
+
+
+ // clear the temporary arrays lists
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ delete arrayofArrayV[k];
+ delete arrayofCharge[k];
+ }
+
+ AliInfo(" done");
+ fInitLookUp = kTRUE;
+
+}
+
+void AliTPCFCVoltError3D::Print(const Option_t* option) const {
+ //
+ // Print function to check the settings of the Rod shifts and the rotated clips
+ // option=="a" prints the C0 and C1 coefficents for calibration purposes
+ //
+
+ TString opt = option; opt.ToLower();
+ printf("%s\n",GetTitle());
+ printf(" - ROD shifts (residual voltage settings): 40V correspond to 1mm of shift\n");
+ Int_t count = 0;
+ printf(" : A-side - IFC (Rod & Strips) \n ");
+ for (Int_t i=0; i<18;i++) {
+ if (fRodVoltShiftA[i]!=0) {
+ printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftA[i]);
+ count++;
+ }
+ if (count==0&&i==17)
+ printf("-> all at 0 V\n");
+ else if (i==17){
+ printf("\n");
+ count=0;
+ }
+ }
+ printf(" : C-side - IFC (Rod & Strips) \n ");
+ for (Int_t i=0; i<18;i++) {
+ if (fRodVoltShiftC[i]!=0) {
+ printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftC[i]);
+ count++;
+ }
+ if (count==0&&i==17)
+ printf("-> all at 0 V\n");
+ else if (i==17){
+ printf("\n");
+ count=0;
+ }
+ }
+ printf(" : A-side - OFC (only Rod) \n ");
+ for (Int_t i=18; i<36;i++) {
+ if (fRodVoltShiftA[i]!=0) {
+ printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftA[i]);
+ count++;
+ }
+ if (count==0&&i==35)
+ printf("-> all at 0 V\n");
+ else if (i==35) {
+ printf("\n");
+ count=0;
+ }
+ }
+ printf(" : C-side - OFC (only Rod) \n ");
+ for (Int_t i=18; i<36;i++) {
+ if (fRodVoltShiftC[i]!=0) {
+ printf("Rod %2d: %3.1f V \t",i,fRodVoltShiftC[i]);
+ count++;
+ }
+ if (count==0&&i==35)
+ printf("-> all at 0 V\n");
+ else if (i==35){
+ printf("\n");
+ count=0;
+ }
+ }
+
+ printf(" - Rotated clips (residual voltage settings): 40V correspond to 1mm of 'offset'\n");
+ if (fRotatedClipVoltA[0]!=0) { printf(" A side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[0]); count++; }
+ if (fRotatedClipVoltA[1]!=0) { printf(" A side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[1]); count++; }
+ if (fRotatedClipVoltC[0]!=0) { printf(" C side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[0]); count++; }
+ if (fRotatedClipVoltC[1]!=0) { printf(" C side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[1]); count++; }
+ if (count==0)
+ printf(" -> no rotated clips \n");
+
+ count=0;
+ printf(" - Copper ROD shifts (without strips):\n");
+ printf(" : A-side - OFC (Copper Rod shift) \n ");
+ for (Int_t i=0; i<18;i++) {
+ if (fOFCRodShiftA[i]!=0) {
+ printf("Rod %2d: %3.1f V \t",i,fOFCRodShiftA[i]);
+ count++;
+ }
+ if (count==0&&i==17)
+ printf("-> all at 0 V\n");
+ else if (i==17){
+ printf("\n");
+ count=0;
+ }
+ }
+ printf(" : C-side - OFC (Copper Rod shift) \n ");
+ for (Int_t i=0; i<18;i++) {
+ if (fOFCRodShiftC[i]!=0) {
+ printf("Rod %2d: %3.1f V \t",i,fOFCRodShiftC[i]);
+ count++;
+ }
+ if (count==0&&i==17)
+ printf("-> all at 0 V\n");
+ else if (i==17){
+ printf("\n");
+ count=0;
+ }
+ }
+
+ 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 InitFCVoltError3D() ...");
+
+}
--- /dev/null
+#ifndef ALITPCFCVOLTERROR3D_H
+#define ALITPCFCVOLTERROR3D_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCFCVoltError3D class //
+// date: 01/06/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliTPCCorrection.h"
+
+
+class AliTPCFCVoltError3D : public AliTPCCorrection {
+public:
+ AliTPCFCVoltError3D();
+ virtual ~AliTPCFCVoltError3D();
+
+ // initialization and update functions
+ virtual void Init();
+ virtual void Update(const TTimeStamp &timeStamp);
+
+ // common setters and getters for tangled ExB effect
+ virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2) {
+ fT1=t1; fT2=t2;
+ const Double_t wt0=t2*omegaTau; fC0=1./(1.+wt0*wt0);
+ const Double_t wt1=t1*omegaTau; fC1=wt1/(1.+wt1*wt1);
+ };
+ void SetC0C1(Float_t c0,Float_t c1) {fC0=c0;fC1=c1;} // CAUTION: USE WITH CARE
+ Float_t GetC0() const {return fC0;}
+ Float_t GetC1() const {return fC1;}
+
+ // setters and getters
+
+ // Set rod shift in Voltage equivalents (40V ~ 1mm)
+ // rod numbers: 0-17 (IFC), 18-35 (OFC)
+ // note: strips move accordingly
+ void SetRodVoltShiftA(Int_t rod, Float_t voltOffset) {fRodVoltShiftA[rod]=voltOffset; fInitLookUp=kFALSE;}
+ void SetRodVoltShiftC(Int_t rod, Float_t voltOffset) {fRodVoltShiftC[rod]=voltOffset; fInitLookUp=kFALSE;}
+ Float_t GetRodVoltShiftA(Int_t i) const {return fRodVoltShiftA[i];}// 0-17: IFC, 18-35; OFC
+ Float_t GetRodVoltShiftC(Int_t i) const {return fRodVoltShiftC[i];}// 0-17: IFC, 18-35; OFC
+
+ // Set rotated clip (just at High Voltage RODs) in Voltage equivalents (40V ~ 1mm)
+ // rod number: 0 (IFC), 1 (OFC)
+ void SetRotatedClipVoltA(Int_t rod, Float_t voltOffset) {fRotatedClipVoltA[rod]=voltOffset; fInitLookUp=kFALSE;}
+ void SetRotatedClipVoltC(Int_t rod, Float_t voltOffset) {fRotatedClipVoltC[rod]=voltOffset; fInitLookUp=kFALSE;}
+ Float_t GetRotatedClipVoltA(Int_t i) const {return fRotatedClipVoltA[i];}// (0,1):(IFC,OFC)
+ Float_t GetRotatedClipVoltC(Int_t i) const {return fRotatedClipVoltC[i];}// (0,1):(IFC,OFC)
+
+ // Set rod shift in Voltage equivalents (40V ~ 1mm)
+ // rod numbers: 0-17 (OFC)
+ // note: strips DO NOT move, only the copper rods do ...
+ void SetOFCRodShiftA(Int_t rod, Float_t voltOffset) {fOFCRodShiftA[rod]=voltOffset; fInitLookUp=kFALSE;}
+ void SetOFCRodShiftC(Int_t rod, Float_t voltOffset) {fOFCRodShiftC[rod]=voltOffset; fInitLookUp=kFALSE;}
+ Float_t GetOFCRodShiftA(Int_t i) const {return fOFCRodShiftA[i];}// 0-17: OFC
+ Float_t GetOFCRodShiftC(Int_t i) const {return fOFCRodShiftC[i];}// 0-17: OFC
+
+
+ void InitFCVoltError3D(); // Fill the lookup tables
+
+ virtual void Print(const Option_t* option="") const;
+
+protected:
+ virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
+
+private:
+
+ AliTPCFCVoltError3D(const AliTPCFCVoltError3D &); // not implemented
+ AliTPCFCVoltError3D &operator=(const AliTPCFCVoltError3D &); // not implemented
+
+ Float_t fC0; // coefficient C0 (compare Jim Thomas's notes for definitions)
+ Float_t fC1; // coefficient C1 (compare Jim Thomas's notes for definitions)
+ Float_t fRodVoltShiftA[36]; // Rod (&strips) shift in Volt (40V~1mm)
+ Float_t fRodVoltShiftC[36]; // Rod (&strips) shift in Volt (40V~1mm)
+ Float_t fRotatedClipVoltA[2]; // rotated clips at HV rod
+ Float_t fRotatedClipVoltC[2]; // rotated clips at HV rod
+ Float_t fOFCRodShiftA[18]; // only Rod shift on OFC
+ Float_t fOFCRodShiftC[18]; // only Rod shift on OFC
+
+ Bool_t fInitLookUp; // flag to check it the Look Up table was created (SUM)
+ Bool_t fInitLookUpBasic[5]; // flag if the basic lookup was created (shifted Rod (IFC,OFC) or rotated clip (IFC,OFC))
+
+
+ TMatrixD *fLookUpErOverEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+ TMatrixD *fLookUpEphiOverEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+ TMatrixD *fLookUpDeltaEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+
+ // basic numbers for the poisson relaxation //can be set individually in each class
+ enum {kRows =257}; // grid size in r direction used in the poisson relaxation // ( 2**n + 1 ) eg. 65, 129, 257 etc.
+ enum {kColumns=129}; // grid size in z direction used in the poisson relaxation // ( 2**m + 1 ) eg. 65, 129, 257 etc.
+ enum {kPhiSlicesPerSector = 10 }; // number of points in phi slices
+ enum {kPhiSlices = 1+kPhiSlicesPerSector*3 }; // number of points in phi for the basic lookup tables
+ enum {kIterations=100}; // Number of iterations within the poisson relaxation
+
+ // ugly way to store "partial" look up tables
+ // needed for the faster calculation of the final distortion table
+
+ // for Rod and Strip shift
+ TMatrixD *fLookUpBasic1ErOverEz[kPhiSlices]; // Array to store electric field integral (int Er/Ez)
+ TMatrixD *fLookUpBasic1EphiOverEz[kPhiSlices]; // Array to store electric field integral (int Ephi/Ez)
+ TMatrixD *fLookUpBasic1DeltaEz[kPhiSlices]; // Array to store electric field integral (int Ez)
+
+ TMatrixD *fLookUpBasic2ErOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic2EphiOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic2DeltaEz[kPhiSlices]; // Array to store electric field integral
+
+ // for rotated clips
+ TMatrixD *fLookUpBasic3ErOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic3EphiOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic3DeltaEz[kPhiSlices]; // Array to store electric field integral
+
+ TMatrixD *fLookUpBasic4ErOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic4EphiOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic4DeltaEz[kPhiSlices]; // Array to store electric field integral
+
+ // for (only rod) shift (just OFC since they are copper plated)
+ TMatrixD *fLookUpBasic5ErOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic5EphiOverEz[kPhiSlices]; // Array to store electric field integral
+ TMatrixD *fLookUpBasic5DeltaEz[kPhiSlices]; // Array to store electric field integral
+
+
+ ClassDef(AliTPCFCVoltError3D,0);
+};
+
+#endif
// Calculate correction in cartesian coordinates
dx[0] = r * TMath::Cos(phi) - x[0];
dx[1] = r * TMath::Sin(phi) - x[1];
- dx[2] = 0.; // z distortion not implemented (1st order distortions)
+ dx[2] = 0.; // z distortion not implemented (1st order distortions) - see e.g. AliTPCBoundaryVoltError-class
//
// virtual destructor
//
+ if (fCorrection) delete fCorrection;
}
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCROCVoltError3D class //
+// The class calculates the space point distortions due to residual voltage //
+// errors on Read Out Chambers of the TPC in 3D. //
+// //
+// The class allows "effective Omega Tau" corrections. //
+// //
+// NOTE: This class is capable of calculating z distortions due to //
+// misalignment and the vd dependency on the residual drift field //
+// //
+// date: 08/08/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+// //
+// Example usage : //
+// AliTPCROCVoltError3D ROCerror; //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "AliMagF.h"
+#include "TGeoGlobalMagField.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
+#include "AliLog.h"
+#include "TMatrixD.h"
+#include "TFile.h"
+
+#include "TMath.h"
+#include "AliTPCROC.h"
+#include "AliTPCROCVoltError3D.h"
+
+ClassImp(AliTPCROCVoltError3D)
+
+AliTPCROCVoltError3D::AliTPCROCVoltError3D()
+ : AliTPCCorrection("ROCVoltErrors","ROC z alignment Errors"),
+ fC0(0.),fC1(0.),
+ fROCdisplacement(kTRUE),
+ fInitLookUp(kFALSE),
+ fROCDataFileName("$(ALICE_ROOT)/TPC/Calib/maps/TPCROCdzSurvey.root"), // standard file name of ROC survey
+ fdzDataLinFit(0)
+{
+ //
+ // default constructor
+ //
+
+ // Array which will contain the solution according to the setted boundary conditions
+ // main input: z alignment of the Read Out chambers
+ // see InitROCVoltError3D() 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);
+ }
+
+ SetROCDataFileName(fROCDataFileName); // initialization of fdzDataLinFit is included
+
+}
+
+AliTPCROCVoltError3D::~AliTPCROCVoltError3D() {
+ //
+ // destructor
+ //
+
+ for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
+ delete fLookUpErOverEz[k];
+ delete fLookUpEphiOverEz[k];
+ delete fLookUpDeltaEz[k];
+ }
+
+ delete fdzDataLinFit;
+}
+
+void AliTPCROCVoltError3D::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);
+
+ InitROCVoltError3D();
+}
+
+void AliTPCROCVoltError3D::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);
+
+}
+
+void AliTPCROCVoltError3D::SetROCDataFileName(char *const fname) {
+ //
+ // Set / load the ROC data (linear fit of ROC misalignments)
+ //
+
+ fROCDataFileName = fname;
+
+ TFile f(fROCDataFileName,"READ");
+ TMatrixD *m = (TMatrixD*) f.Get("dzSurveyLinFitData");
+ TMatrixD &mf = *m;
+
+ // prepare some space
+
+ if (fdzDataLinFit) delete fdzDataLinFit;
+ fdzDataLinFit = new TMatrixD(72,3);
+ TMatrixD &dataIntern = *fdzDataLinFit;
+
+ for (Int_t iroc=0;iroc<72;iroc++) {
+ dataIntern(iroc,0) = mf(iroc,0); // z0 offset
+ dataIntern(iroc,1) = mf(iroc,1); // slope in x
+ dataIntern(iroc,2) = mf(iroc,2); // slope in y
+ }
+
+ f.Close();
+
+ fInitLookUp = kFALSE;
+
+}
+
+void AliTPCROCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
+ //
+ // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
+ //
+
+ if (!fInitLookUp) {
+ AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
+ InitROCVoltError3D();
+ return;
+ }
+
+ Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
+
+ Double_t intEr, intEphi, intDeltaEz;
+ 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);
+ intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
+ fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz );
+
+ // printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
+
+ // Calculate distorted position
+ if ( r > 0.0 ) {
+ phi = phi + ( fC0*intEphi - fC1*intEr ) / r;
+ r = r + ( fC0*intEr + fC1*intEphi );
+ }
+
+ // Calculate correction in cartesian coordinates
+ dx[0] = r * TMath::Cos(phi) - x[0];
+ dx[1] = r * TMath::Sin(phi) - x[1];
+ dx[2] = intDeltaEz; // z distortion - (internally scaled with driftvelocity dependency
+ // on the Ez field plus the actual ROC misalignment (if set TRUE)
+
+}
+
+void AliTPCROCVoltError3D::InitROCVoltError3D() {
+ //
+ // Initialization of the Lookup table which contains the solutions of the
+ // Dirichlet boundary problem
+ // Calculation of the single 3D-Poisson solver is done just if needed
+ // (see basic lookup tables in header file)
+ //
+
+ const Int_t order = 1 ; // Linear interpolation = 1, Quadratic = 2
+ const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
+ const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
+ const Float_t gridSizePhi = TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
+
+ // temporary arrays to create the boundary conditions
+ TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
+ TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
+
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ arrayofArrayV[k] = new TMatrixD(kRows,kColumns) ;
+ arrayofCharge[k] = new TMatrixD(kRows,kColumns) ;
+ arrayofEroverEz[k] = new TMatrixD(kRows,kColumns) ;
+ arrayofEphioverEz[k] = new TMatrixD(kRows,kColumns) ;
+ arrayofDeltaEz[k] = new TMatrixD(kRows,kColumns) ;
+ }
+
+ // list of point as used in the poisson relation and the interpolation (during sum up)
+ Double_t rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ philist[k] = gridSizePhi * k;
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ rlist[i] = fgkIFCRadius + i*gridSizeR ;
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
+ zedlist[j] = j * gridSizeZ ;
+ }
+ }
+ }
+
+ // ==========================================================================
+ // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
+ // Allow for different size grid spacing in R and Z directions
+
+ const Int_t symmetry = 0;
+
+ // Set bondaries and solve Poisson's equation --------------------------
+
+ if ( !fInitLookUp ) {
+
+ AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
+
+ for ( Int_t side = 0 ; side < 2 ; side++ ) { // Solve Poisson3D twice; once for +Z and once for -Z
+
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ TMatrixD &arrayV = *arrayofArrayV[k] ;
+ TMatrixD &charge = *arrayofCharge[k] ;
+
+ //Fill arrays with initial conditions. V on the boundary and Charge in the volume.
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
+ arrayV(i,j) = 0.0 ;
+ charge(i,j) = 0.0 ;
+
+ Float_t radius0 = rlist[i] ;
+ Float_t phi0 = gridSizePhi * k ;
+
+ // To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
+ if ( j == (kColumns-1) )
+ arrayV(i,j) = 0.5* ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
+
+ }
+ }
+
+ for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
+ for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
+ charge(i,j) = 0.0 ;
+ }
+ }
+ }
+
+ // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
+ // Allow for different size grid spacing in R and Z directions
+
+ PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
+ arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
+ kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
+ symmetry, fROCdisplacement) ;
+
+
+ //Interpolate results onto a custom grid which is used just for these calculations.
+ 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]) ; // Symmetric solution in Z that depends only on ABS(Z)
+
+ if ( side == 0 && fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
+ if ( side == 1 && fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
+
+ for ( Int_t i = 0 ; i < kNR ; i++ ) {
+ r = fgkRList[i] ;
+
+ // Interpolate basicLookup tables; once for each rod, then sum the results
+ erOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, arrayofEroverEz );
+ ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, arrayofEphioverEz);
+ deltaEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
+ rlist, zedlist, philist, arrayofDeltaEz );
+
+ if (side == 1) deltaEz(i,j) = - deltaEz(i,j); // negative coordinate system on C side
+
+ } // end r loop
+ }// end z loop
+ }// end phi loop
+
+ if ( side == 0 ) AliInfo(" A side done");
+ if ( side == 1 ) AliInfo(" C side done");
+ } // end side loop
+ }
+
+ // clear the temporary arrays lists
+ for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
+ delete arrayofArrayV[k];
+ delete arrayofCharge[k];
+ delete arrayofEroverEz[k];
+ delete arrayofEphioverEz[k];
+ delete arrayofDeltaEz[k];
+ }
+
+
+ fInitLookUp = kTRUE;
+
+}
+
+
+Float_t AliTPCROCVoltError3D::GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0) {
+ //
+ // Returns the dz alignment data (in voltage equivalents) at
+ // the given position
+ //
+
+ Float_t xp = r0*TMath::Cos(phi0);
+ Float_t yp = r0*TMath::Sin(phi0);
+
+ // phi0 should be between 0 and 2pi
+ if (phi0<0) phi0+=TMath::TwoPi();
+ Int_t roc = (Int_t)TMath::Floor((TMath::RadToDeg()*phi0)/20);
+ if (side==1) roc+=18; // C side
+ if (r0>132) roc+=36; // OROC
+
+ // linear-plane data: z = z0 + kx*x + ky*y
+ TMatrixD &fitData = *fdzDataLinFit;
+ Float_t dz = fitData(roc,0)+fitData(roc,1)*xp + fitData(roc,2)*yp; // value in cm
+
+ // aproximated Voltage-offset-aquivalent to the z misalignment
+ // (linearly scaled with the z position)
+ Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
+ Float_t voltOff = dz*ezField; // values in "Volt equivalents"
+
+ return voltOff;
+}
+
+TH2F * AliTPCROCVoltError3D::CreateHistoOfZSurvey(Int_t side, Int_t nx, Int_t ny) {
+ //
+ // return a simple histogramm containing the input to the poisson solver
+ // (z positions of the Read-out chambers, linearly interpolated)
+
+ char hname[100];
+ if (side==0) sprintf(hname,"survey_dz_Aside");
+ if (side==1) sprintf(hname,"survey_dz_Cside");
+
+ TH2F *h = new TH2F(hname,hname,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 r0 = TMath::Sqrt(xp*xp+yp*yp);
+ Float_t phi0 = TMath::ATan2(yp,xp);
+
+ Float_t dz = GetROCVoltOffset(side,r0,phi0); // in [volt]
+
+ Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
+ dz = dz/ezField; // in [cm]
+
+ if (85.<=r0 && r0<=245.) {
+ h->SetBinContent(ix,iy,dz);
+ } else {
+ h->SetBinContent(ix,iy,0.);
+ }
+ }
+ }
+
+ h->GetXaxis()->SetTitle("x [cm]");
+ h->GetYaxis()->SetTitle("y [cm]");
+ h->GetZaxis()->SetTitle("dz [cm]");
+ h->SetStats(0);
+ // h->DrawCopy("colz");
+
+ return h;
+}
+
+void AliTPCROCVoltError3D::Print(const Option_t* option) const {
+ //
+ // Print function to check the settings of the Rod shifts and the rotated clips
+ // option=="a" prints the C0 and C1 coefficents for calibration purposes
+ //
+
+ TString opt = option; opt.ToLower();
+ printf("%s\n",GetTitle());
+ printf(" - Voltage settings on the TPC Read-Out chambers - linearly interpolated\n");
+ printf(" info: Check the following data-file for more details: %s \n",fROCDataFileName);
+
+ 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 InitROCVoltError3D() ...");
+
+}
--- /dev/null
+#ifndef ALITPCROCVOLTERROR3D_H
+#define ALITPCROCVOLTERROR3D_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCROCVoltError3D class //
+// date: 01/06/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliTPCCorrection.h"
+#include "TH2F.h"
+
+
+class AliTPCROCVoltError3D : public AliTPCCorrection {
+public:
+ AliTPCROCVoltError3D();
+ virtual ~AliTPCROCVoltError3D();
+
+ // initialization and update functions
+ virtual void Init();
+ virtual void Update(const TTimeStamp &timeStamp);
+
+ // common setters and getters for tangled ExB effect
+ virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2) {
+ fT1=t1; fT2=t2;
+ const Double_t wt0=t2*omegaTau; fC0=1./(1.+wt0*wt0);
+ const Double_t wt1=t1*omegaTau; fC1=wt1/(1.+wt1*wt1);
+ };
+ void SetC0C1(Float_t c0,Float_t c1) {fC0=c0;fC1=c1;} // CAUTION: USE WITH CARE
+ Float_t GetC0() const {return fC0;}
+ Float_t GetC1() const {return fC1;}
+
+ // setters and getters
+ void SetROCDataFileName(char *const fname);
+ char* GetROCDataFileName() const {return fROCDataFileName;}
+
+ // flag to wheter or not include the z aligment in the dz calculation
+ // if FALSE, the dz offset is purely due to the electric field change
+ void SetROCDisplacement(Bool_t flag) { fROCdisplacement = flag; fInitLookUp=kFALSE;}
+ Bool_t GetROCDisplacement() const { return fROCdisplacement; }
+
+
+ void InitROCVoltError3D(); // Fill the lookup tables
+
+ Float_t GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0);
+ TH2F* CreateHistoOfZSurvey(Int_t side, Int_t nx=250, Int_t ny=250);
+
+ virtual void Print(const Option_t* option="") const;
+
+protected:
+ virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
+
+private:
+
+ AliTPCROCVoltError3D(const AliTPCROCVoltError3D &); // not implemented
+ AliTPCROCVoltError3D &operator=(const AliTPCROCVoltError3D &); // not implemented
+
+ Float_t fC0; // coefficient C0 (compare Jim Thomas's notes for definitions)
+ Float_t fC1; // coefficient C1 (compare Jim Thomas's notes for definitions)
+
+ Bool_t fROCdisplacement; // flag on wheter to consider the ROC displacement
+ // when calculating the z distortions
+ Bool_t fInitLookUp; // flag to check it the Look Up table was created (SUM)
+
+ TMatrixD *fLookUpErOverEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+ TMatrixD *fLookUpEphiOverEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+ TMatrixD *fLookUpDeltaEz[kNPhi]; // Array to store electric field integral (int Er/Ez)
+
+ char *fROCDataFileName; // filename of the survey data containing the lin Fit values
+ TMatrixD *fdzDataLinFit; // Linear fits of dz survey points (each sector=72) (z0,slopeX,slopeY)
+
+ // basic numbers for the poisson relaxation //can be set individually in each class
+ enum {kRows =257}; // grid size in r direction used in the poisson relaxation // ( 2**n + 1 ) eg. 65, 129, 257 etc.
+ enum {kColumns=129}; // grid size in z direction used in the poisson relaxation // ( 2**m + 1 ) eg. 65, 129, 257 etc.
+ enum {kPhiSlicesPerSector=6}; // phi slices per sector
+ enum {kPhiSlices = 18*kPhiSlicesPerSector }; // number of points in phi for the basic lookup tables
+ enum {kIterations=100}; // Number of iterations within the poisson relaxation
+
+ ClassDef(AliTPCROCVoltError3D,0);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCSpaceCharge class //
+// The class calculates the space point distortions due to a space charge //
+// effect .... //
+// The class allows "effective Omega Tau" corrections. //
+// //
+// NOTE: This class is capable of calculating z distortions due to //
+// drift velocity change in dependence of the electric field!!! //
+// //
+// date: 23/08/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+// //
+// Example usage: //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "AliMagF.h"
+#include "TGeoGlobalMagField.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
+#include "AliLog.h"
+#include "TMatrixD.h"
+
+#include "TMath.h"
+#include "AliTPCROC.h"
+#include "AliTPCSpaceCharge.h"
+
+ClassImp(AliTPCSpaceCharge)
+
+AliTPCSpaceCharge::AliTPCSpaceCharge()
+ : AliTPCCorrection("SpaceCharge2D","Space Charge 2D"),
+ fC0(0.),fC1(0.),fCorrectionFactor(0.001),
+ fInitLookUp(kFALSE)
+{
+ //
+ // default constructor
+ //
+
+}
+
+AliTPCSpaceCharge::~AliTPCSpaceCharge() {
+ //
+ // default destructor
+ //
+}
+
+
+
+void AliTPCSpaceCharge::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);
+
+ InitSpaceChargeDistortion(); // fill the look up table
+}
+
+void AliTPCSpaceCharge::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 AliTPCSpaceCharge::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! Perform the inizialisation now ...");
+ InitSpaceChargeDistortion();
+ }
+ 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!");
+
+ // Efield is symmetric in phi - 2D calculation
+ intEphi = 0.0;
+ // Get the E field integrals
+ Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr );
+ // Get DeltaEz field integral
+ Interpolate2DEdistortion( order, r, z, fLookUpDeltaEz, intdEz );
+
+
+ // 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;
+
+ // 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 - (internally scaled with driftvelocity dependency
+ // on the Ez field
+
+}
+
+void AliTPCSpaceCharge::InitSpaceChargeDistortion() {
+ //
+ // Initialization of the Lookup table which contains the solutions of the
+ // poisson problem
+ //
+
+ const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
+ const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
+
+ TMatrixD voltArray(kRows,kColumns); // dummy boundary vectors
+ TMatrixD chargeDensity(kRows,kColumns); // charge
+ TMatrixD arrayErOverEz(kRows,kColumns); // solution in Er
+ TMatrixD arrayDeltaEz(kRows,kColumns); // solution in Ez
+
+ Double_t rList[kRows], zedList[kColumns] ;
+
+ // Fill arrays with initial conditions. V on the boundary and ChargeDensity in the volume.
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) {
+ Double_t zed = j*gridSizeZ ;
+ zedList[j] = zed ;
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ Double_t radius = fgkIFCRadius + i*gridSizeR ;
+ rList[i] = radius ;
+ voltArray(i,j) = 0; // Initialize voltArray to zero - not used in this class
+ chargeDensity(i,j) = 0; // Initialize ChargeDensity to zero
+ }
+ }
+
+ // Fill the initial conditions
+ for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
+ Double_t zed = j*gridSizeZ ;
+ for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
+ Double_t radius = fgkIFCRadius + i*gridSizeR ;
+
+ Double_t zterm = (fgkTPCZ0-zed) * (fgkOFCRadius*fgkOFCRadius - fgkIFCRadius*fgkIFCRadius) / fgkTPCZ0 ;
+ // for 1/R**2 charge density in the TPC; then integrated in Z due to drifting ions
+ chargeDensity(i,j) = zterm / ( TMath::Log(fgkOFCRadius/fgkIFCRadius) * ( radius*radius ) ) ;
+ }
+ }
+
+
+ // Solve the electrosatic problem in 2D
+
+ PoissonRelaxation2D( voltArray, chargeDensity, arrayErOverEz, arrayDeltaEz, kRows, kColumns, kIterations ) ;
+
+ //Interpolate results onto standard grid for Electric Fields
+ Int_t ilow=0, jlow=0 ;
+ Double_t z,r;
+ Float_t saveEr[2], saveEz[2] ;
+ for ( Int_t i = 0 ; i < kNZ ; ++i ) {
+ z = TMath::Abs( fgkZList[i] ) ; // assume symmetric behaviour on A and C side
+ for ( Int_t j = 0 ; j < kNR ; ++j ) {
+
+ // Linear interpolation !!
+ r = fgkRList[j] ;
+ Search( kRows, rList, r, ilow ) ; // Note switch - R in rows and Z in columns
+ Search( kColumns, zedList, z, jlow ) ;
+ if ( ilow < 0 ) ilow = 0 ; // check if out of range
+ if ( jlow < 0 ) jlow = 0 ;
+ if ( ilow + 1 >= kRows - 1 ) ilow = kRows - 2 ;
+ if ( jlow + 1 >= kColumns - 1 ) jlow = kColumns - 2 ;
+
+ saveEr[0] = arrayErOverEz(ilow,jlow) +
+ (arrayErOverEz(ilow,jlow+1)-arrayErOverEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEr[1] = arrayErOverEz(ilow+1,jlow) +
+ (arrayErOverEz(ilow+1,jlow+1)-arrayErOverEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEz[0] = arrayDeltaEz(ilow,jlow) +
+ (arrayDeltaEz(ilow,jlow+1)-arrayDeltaEz(ilow,jlow))*(z-zedList[jlow])/gridSizeZ ;
+ saveEz[1] = arrayDeltaEz(ilow+1,jlow) +
+ (arrayDeltaEz(ilow+1,jlow+1)-arrayDeltaEz(ilow+1,jlow))*(z-zedList[jlow])/gridSizeZ ;
+
+
+ fLookUpErOverEz[i][j] = saveEr[0] + (saveEr[1]-saveEr[0])*(r-rList[ilow])/gridSizeR ;
+ fLookUpDeltaEz[i][j] = saveEz[0] + (saveEz[1]-saveEz[0])*(r-rList[ilow])/gridSizeR ;
+ }
+ }
+
+ fInitLookUp = kTRUE;
+
+}
+
+void AliTPCSpaceCharge::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 effects assuming a radial symmetric z over r^2 SC-distribution.\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 InitSpaceChargeDistortion() ...");
+
+}
--- /dev/null
+#ifndef ALITPCSPACECHARGE_H
+#define ALITPCSPACECHARGE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCSpaceCharge class //
+// date: 19/06/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliTPCCorrection.h"
+
+
+class AliTPCSpaceCharge : public AliTPCCorrection {
+public:
+ AliTPCSpaceCharge();
+ virtual ~AliTPCSpaceCharge();
+
+ // initialization and update functions
+ virtual void Init();
+ virtual void Update(const TTimeStamp &timeStamp);
+
+
+ // common setters and getters for tangled ExB effect
+ virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2) {
+ fT1=t1; fT2=t2;
+ const Double_t wt0=t2*omegaTau; fC0=1./(1.+wt0*wt0);
+ const Double_t wt1=t1*omegaTau; fC1=wt1/(1.+wt1*wt1);
+ };
+ void SetC0C1(Float_t c0,Float_t c1) {fC0=c0;fC1=c1;} // CAUTION: USE WITH CARE
+ Float_t GetC0() const {return fC0;}
+ Float_t GetC1() const {return fC1;}
+
+ // setters and getters for conical
+ void SetCorrectionFactor(Float_t correctionFactor) {fCorrectionFactor=correctionFactor;}
+ Float_t GetCorrectionFactor() const {return fCorrectionFactor;}
+
+ void InitSpaceChargeDistortion();
+
+ virtual void Print(const Option_t* option="") const;
+
+protected:
+ virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
+
+private:
+ Float_t fC0; // coefficient C0 (compare Jim Thomas's notes for definitions)
+ Float_t fC1; // coefficient C1 (compare Jim Thomas's notes for definitions)
+ Float_t fCorrectionFactor; // Space Charge Correction factor in comparison to initialized
+ // look up table which was created for M_mb = 900 and IR = 3000
+ // compare Internal Note Nr: ???
+
+ Bool_t fInitLookUp; // flag to check it the Look Up table was created
+
+ Double_t fLookUpErOverEz[kNZ][kNR]; // Array to store electric field integral (int Er/Ez)
+ Double_t fLookUpDeltaEz[kNZ][kNR]; // Array to store electric field integral (int Delta Ez)
+
+ // basic numbers for the poisson relaxation //can be set individually in each class
+ enum {kRows =257}; // grid size in r direction used in the poisson relaxation // ( 2**n + 1 ) eg. 65, 129, 257 etc.
+ enum {kColumns=129}; // grid size in z direction used in the poisson relaxation // ( 2**m + 1 ) eg. 65, 129, 257 etc.
+ enum {kIterations=100}; // Number of iterations within the poisson relaxation
+
+ ClassDef(AliTPCSpaceCharge,0);
+};
+
+#endif
TVectorD vPosG(3); //vertex position
TVectorD vPosL(3); // vertex position in the TPC local system
TVectorF vImpact(2); //track impact parameter
- Double_t tofSignal= fCurrentTrack->GetTOFsignal(); // tof signal
+ // Double_t tofSignal= fCurrentTrack->GetTOFsignal(); // tof signal
TVectorF tpcPosG(3); // global position of track at the middle of fXmiddle
Double_t lphi = TMath::ATan2(pyf[0],fXmiddle); // expected local phi angle - if vertex at 0
Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi; // expected global phi if vertex at 0
#pragma link C++ class AliTPCComposedCorrection+;
#pragma link C++ class AliTPCExBBShape+;
#pragma link C++ class AliTPCExBTwist+;
-#pragma link C++ class AliTPCExBConical+;
#pragma link C++ class AliTPCGGVoltError+;
+#pragma link C++ class AliTPCFCVoltError3D+;
+#pragma link C++ class AliTPCROCVoltError3D+;
#pragma link C++ class AliTPCBoundaryVoltError+;
#pragma link C++ class AliTPCCalibGlobalMisalignment+;
+#pragma link C++ class AliTPCSpaceCharge+;
+
#pragma link C++ class AliTPCExBEffective+;
#pragma link C++ class AliTPCExBEffectiveSector+;
AliTPCkalmanTime.cxx AliESDcosmic.cxx AliTPCPointCorrection.cxx AliTPCTransformation.cxx \
AliTPCkalmanFit.cxx AliTPCLaserTrack.cxx AliTPCcalibBase.cxx \
AliTPCCorrection.cxx AliTPCInverseCorrection.cxx AliTPCComposedCorrection.cxx \
- AliTPCExBBShape.cxx AliTPCExBTwist.cxx AliTPCGGVoltError.cxx AliTPCBoundaryVoltError.cxx AliTPCExBConical.cxx AliXRDPROOFtoolkit.cxx AliTPCCalibGlobalMisalignment.cxx AliTPCExBEffective.cxx AliTPCExBEffectiveSector.cxx
+ AliTPCExBBShape.cxx AliTPCExBTwist.cxx AliTPCGGVoltError.cxx AliTPCFCVoltError3D.cxx AliTPCROCVoltError3D.cxx AliTPCBoundaryVoltError.cxx AliTPCSpaceCharge.cxx AliXRDPROOFtoolkit.cxx AliTPCCalibGlobalMisalignment.cxx AliTPCExBEffective.cxx AliTPCExBEffectiveSector.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff