--- /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. *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCBoundaryVoltError class //
+// The class calculates the space point distortions due to residual voltage //
+// errors on the main boundaries of the TPC. For example, the inner vessel //
+// of the TPC is shifted by a certain amount, whereas the ROCs on the A side//
+// and the ROCs on the C side follow this mechanical shift (at the inner //
+// vessel) in z direction (see example below). This example is commonly //
+// named "conical deformation" of the TPC field cage. //
+// //
+// The class allows "effective Omega Tau" corrections. //
+// //
+// NOTE: This class is not capable of calculation z distortions due to //
+// drift velocity change in dependence of the electric field!!! //
+// //
+// date: 01/06/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+// //
+// Example usage (e.g +1mm shift of "conical deformation") //
+// AliTPCBoundaryVoltError bve; //
+// Float_t boundA[8] = {-40,-40,-40,0,0,0,0,-40}; // voltages A-side //
+// Float_t boundC[6] = { 40, 40, 40,0,0,0}; // voltages C-side //
+// bve.SetBoundariesA(boundA); //
+// bve.SetBoundariesC(boundC); //
+// bve.SetOmegaTauT1T2(0.32,1.,1.); // values ideally from OCDB //
+// // initialization of the look up //
+// bve.InitBoundaryVoltErrorDistortion(); //
+// // plot dRPhi distortions ... //
+// bve.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2"); //
+//////////////////////////////////////////////////////////////////////////////
+
+#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 "AliTPCBoundaryVoltError.h"
+
+ClassImp(AliTPCBoundaryVoltError)
+
+AliTPCBoundaryVoltError::AliTPCBoundaryVoltError()
+ : AliTPCCorrection("BoundaryVoltError","Boundary Voltage Error"),
+ fC0(0.),fC1(0.),
+ fInitLookUp(kFALSE)
+{
+ //
+ // default constructor
+ //
+ for (Int_t i=0; i<8; i++){
+ fBoundariesA[i]= 0;
+ if (i<6) fBoundariesC[i]= 0;
+ }
+}
+
+AliTPCBoundaryVoltError::~AliTPCBoundaryVoltError() {
+ //
+ // default destructor
+ //
+}
+
+
+
+void AliTPCBoundaryVoltError::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);
+
+ InitBoundaryVoltErrorDistortion();
+}
+
+void AliTPCBoundaryVoltError::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 AliTPCBoundaryVoltError::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) AliError("Lookup table was not initialized! You should do 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 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
+
+
+ intEphi = 0.0; // Efield is symmetric in phi - 2D calculation
+
+ 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 E field integral
+ Interpolate2DEdistortion( order, r, z, fLookUpErOverEz, intEr );
+
+ // 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] = 0.; // z distortion not implemented (1st order distortions)
+
+}
+
+void AliTPCBoundaryVoltError::InitBoundaryVoltErrorDistortion() {
+ //
+ // Initialization of the Lookup table which contains the solutions of the
+ // Dirichlet boundary problem
+ //
+
+ const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
+ const Float_t gridSizeZ = fgkTPCZ0 / (kColumns-1) ;
+
+ TMatrixD voltArrayA(kRows,kColumns), voltArrayC(kRows,kColumns); // boundary vectors
+ TMatrixD chargeDensity(kRows,kColumns); // dummy charge
+ TMatrixD arrayErOverEzA(kRows,kColumns), arrayErOverEzC(kRows,kColumns); // solution
+
+ 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 ;
+ voltArrayA(i,j) = 0; // Initialize voltArrayA to zero
+ voltArrayC(i,j) = 0; // Initialize voltArrayC to zero
+ chargeDensity(i,j) = 0; // Initialize ChargeDensity to zero - not used in this class
+ }
+ }
+
+
+ // check if boundary values are the same for the C side (for later, saving some calculation time)
+
+ Int_t symmetry = -1; // assume that A and C side are identical (but anti-symmetric!) // e.g conical deformation
+ Int_t sVec[8];
+
+ // check if boundaries are different (regardless the sign)
+ for (Int_t i=0; i<8; i++) {
+ if ((TMath::Abs(fBoundariesA[i]) - TMath::Abs(fBoundariesC[i])) > 1e-5) symmetry = 0;
+ sVec[i] = (Int_t) (TMath::Sign((Float_t)1.,fBoundariesA[i])*TMath::Sign((Float_t)1.,fBoundariesC[i])); // == -1 for anti-symmetry
+ }
+ if (symmetry==-1) { // still the same values?
+ // check the kind of symmetry , if even ...
+ if (sVec[0]==1 && sVec[1]==1 && sVec[2]==1 && sVec[3]==1 && sVec[4]==1 && sVec[5]==1 && sVec[6]==1 && sVec[7]==1 )
+ symmetry = 1;
+ else if (sVec[0]==-1 && sVec[1]==-1 && sVec[2]==-1 && sVec[3]==-1 && sVec[4]==-1 && sVec[5]==-1 && sVec[6]==-1 && sVec[7]==-1 )
+ symmetry = -1;
+ else
+ symmetry = 0; // some of the values differ in the sign -> neither symmetric nor antisymmetric
+ }
+
+
+
+ // Solve the electrosatic problem in 2D
+
+ // Fill the complete Boundary vectors
+ // Start at IFC at CE and work anti-clockwise through IFC, ROC, OFC, and CE (clockwise for C side)
+ for ( Int_t j = 0 ; j < kColumns ; j++ ) {
+ Double_t zed = j*gridSizeZ ;
+ for ( Int_t i = 0 ; i < kRows ; i++ ) {
+ Double_t radius = fgkIFCRadius + i*gridSizeR ;
+
+ // 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))
+ + 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))
+ + 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))
+ + 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))
+ + fBoundariesC[7] ; // CE
+
+ }
+ }
+ }
+
+ voltArrayA(0,0) *= 0.5 ; // Use average boundary condition at corner
+ voltArrayA(kRows-1,0) *= 0.5 ; // Use average boundary condition at corner
+ voltArrayA(0,kColumns-1) *= 0.5 ; // Use average boundary condition at corner
+ voltArrayA(kRows-1,kColumns-1)*= 0.5 ; // Use average boundary condition at corner
+
+ if (symmetry==0) {
+ voltArrayC(0,0) *= 0.5 ; // Use average boundary condition at corner
+ voltArrayC(kRows-1,0) *= 0.5 ; // Use average boundary condition at corner
+ voltArrayC(0,kColumns-1) *= 0.5 ; // Use average boundary condition at corner
+ voltArrayC(kRows-1,kColumns-1)*= 0.5 ; // Use average boundary condition at corner
+ }
+
+
+ // always solve the problem on the A side
+ PoissonRelaxation2D( voltArrayA, chargeDensity, arrayErOverEzA, kRows, kColumns, kIterations ) ;
+
+ 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);
+ }
+ }
+ } 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 ) ;
+ }
+
+ //Interpolate results onto standard grid for Electric Fields
+ Int_t ilow=0, jlow=0 ;
+ Double_t z,r;
+ Float_t saveEr[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 ;
+ }
+ }
+
+ /* delete [] saveEr;
+ delete [] sVec;
+ delete [] rList;
+ delete [] zedList;
+ */
+
+ fInitLookUp = kTRUE;
+
+}
+
+void AliTPCBoundaryVoltError::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(" - Voltage settings (on the TPC boundaries) - linearly interpolated\n");
+ printf(" : A-side (anti-clockwise)\n");
+ printf(" (0,1):\t IFC (CE) : %3.1f V \t IFC (ROC): %3.1f V \n",fBoundariesA[0],fBoundariesA[1]);
+ printf(" (2,3):\t ROC (IFC): %3.1f V \t ROC (OFC): %3.1f V \n",fBoundariesA[2],fBoundariesA[3]);
+ printf(" (4,5):\t OFC (ROC): %3.1f V \t OFC (CE) : %3.1f V \n",fBoundariesA[4],fBoundariesA[5]);
+ printf(" (6,7):\t CE (OFC): %3.1f V \t CE (IFC): %3.1f V \n",fBoundariesA[6],fBoundariesA[7]);
+ printf(" : C-side (clockwise)\n");
+ printf(" (0,1):\t IFC (CE) : %3.1f V \t IFC (ROC): %3.1f V \n",fBoundariesC[0],fBoundariesC[1]);
+ printf(" (2,3):\t ROC (IFC): %3.1f V \t ROC (OFC): %3.1f V \n",fBoundariesC[2],fBoundariesC[3]);
+ printf(" (4,5):\t OFC (ROC): %3.1f V \t OFC (CE) : %3.1f V \n",fBoundariesC[4],fBoundariesC[5]);
+ printf(" (6,7):\t CE (OFC): %3.1f V \t CE (IFC): %3.1f V \n",fBoundariesC[6],fBoundariesC[7]);
+
+ // Check wether the settings of the Central Electrode agree (on the A and C side)
+ // Note: they have to be antisymmetric!
+ if (( TMath::Abs(fBoundariesA[6]+fBoundariesC[6])>1e-5) || ( TMath::Abs(fBoundariesA[7]+fBoundariesC[7])>1e-5 ) ){
+ AliWarning("Boundary parameters for the Central Electrode (CE) are not anti-symmetric! HOW DID YOU MANAGE THAT?");
+ AliWarning("Congratulations, you just splitted the Central Electrode of the TPC!");
+ AliWarning("Non-physical settings of the boundary parameter at the Central Electrode");
+ }
+
+ 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 InitBoundaryVoltErrorDistortion() ...");
+
+}
+
+
+void AliTPCBoundaryVoltError::SetBoundariesA(Float_t boundariesA[8]){
+ //
+ // set voltage errors on the TPC boundaries - A side
+ //
+ // Start at IFC at the Central electrode and work anti-clockwise (clockwise for C side) through
+ // IFC, ROC, OFC, and CE. The boundary conditions are currently defined to be a linear
+ // interpolation between pairs of numbers in the Boundary (e.g. fBoundariesA) vector.
+ // The first pair of numbers represent the beginning and end of the Inner Field cage, etc.
+ // The unit of the error potential vector is [Volt], whereas 1mm shift of the IFC would
+ // correspond to ~ 40 V
+ //
+ // Note: The setting for the CE will be passed to the C side!
+
+ for (Int_t i=0; i<8; i++) {
+ fBoundariesA[i]= boundariesA[i];
+ if (i>5) fBoundariesC[i]= -boundariesA[i]; // setting for the CE is passed to C side
+ }
+
+}
+void AliTPCBoundaryVoltError::SetBoundariesC(Float_t boundariesC[6]){
+ //
+ // set voltage errors on the TPC boundaries - A 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
+ // interpolation between pairs of numbers in the Boundary (e.g. fBoundariesC) vector.
+ // The first pair of numbers represent the beginning and end of the Inner Field cage, etc.
+ // The unit of the error potential vector is [Volt], whereas 1mm shift of the IFC would
+ // correspond to ~ 40 V
+ //
+ // Note: The setting for the CE will be taken from the A side (pos 6 and 7)!
+
+ for (Int_t i=0; i<6; i++) {
+ fBoundariesC[i]= boundariesC[i];
+ }
+
+}
--- /dev/null
+#ifndef ALITPCBOUNDARYVOLTERROR_H
+#define ALITPCBOUNDARYVOLTERROR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// AliTPCBoundaryVoltError class //
+// date: 01/06/2010 //
+// Authors: Jim Thomas, Stefan Rossegger //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliTPCCorrection.h"
+
+
+class AliTPCBoundaryVoltError : public AliTPCCorrection {
+public:
+ AliTPCBoundaryVoltError();
+ virtual ~AliTPCBoundaryVoltError();
+
+ // 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 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];}
+
+ void InitBoundaryVoltErrorDistortion();
+
+ 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 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
+
+ Double_t fLookUpErOverEz[kNZ][kNR]; // 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=257}; // grid size in r 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(AliTPCBoundaryVoltError,0);
+};
+
+#endif
#include <AliCDBMetaData.h>
#include "TVectorD.h"
+
+#include "TRandom.h"
+#include "AliExternalTrackParam.h"
+#include "AliTrackPointArray.h"
+#include "TDatabasePDG.h"
+#include "AliTrackerBase.h"
+#include "AliTPCROC.h"
+#include "THnSparse.h"
+
#include "TRandom.h"
#include "AliTPCTransform.h"
#include "AliTPCcalibDB.h"
#include "AliTrackerBase.h"
#include "AliTPCROC.h"
#include "THnSparse.h"
+
#include "AliTPCLaserTrack.h"
#include "AliTPCCorrection.h"
+#include "AliLog.h"
ClassImp(AliTPCCorrection)
h->SetBinContent(iz,ir,r1-r0);
}
}
- printf("SDF\n");
return h;
}
}
+void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
+ TMatrixD &arrayErOverEz, const Int_t rows,
+ const Int_t columns, const Int_t iterations ) {
+ //
+ // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
+ //
+ // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
+ // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
+ // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
+ // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
+ // you wish to solve Laplaces equation however it should not contain random numbers or you will get
+ // random numbers back as a solution.
+ // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
+ // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
+ // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
+ // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
+ // number of points and solves the problem again. This happens several times until the maximum number
+ // of points has been included in the array.
+ //
+ // 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.
+ //
+ // Original code by Jim Thomas (STAR TPC Collaboration)
+ //
+
+ 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 gridSizeZ = fgkTPCZ0 / (columns-1) ;
+ const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
+
+ TMatrixD arrayEr(rows,columns) ;
+ TMatrixD arrayEz(rows,columns) ;
+
+ //Check that number of rows and columns is suitable for a binary expansion
+
+ if ( !IsPowerOfTwo(rows-1) ) {
+ AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
+ return;
+ }
+ if ( !IsPowerOfTwo(columns-1) ) {
+ AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
+ 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 size of the matrix to speed up the solution of the problem
+
+ Int_t iOne = (rows-1)/4 ;
+ Int_t jOne = (columns-1)/4 ;
+ // Solve for N in 2**N, add one.
+ Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
+
+ for ( Int_t count = 0 ; count < loops ; count++ ) {
+ // Loop while the matrix expands & the resolution increases.
+
+ Float_t tempGridSizeR = gridSizeR * iOne ;
+ Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
+ Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
+
+ // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef1(rows) ;
+ std::vector<float> coef2(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);
+ }
+
+ TMatrixD sumChargeDensity(rows,columns) ;
+
+ for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+ Float_t radius = fgkIFCRadius + iOne*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 ;
+ // Note that this is cylindrical geometry
+ 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++ ) {
+ // Solve Poisson's Equation
+ // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
+ // as interations increase.
+ Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
+ Float_t overRelaxM1 = overRelax - 1.0 ;
+ Float_t overRelaxtempFourth, overRelaxcoef5 ;
+ overRelaxtempFourth = overRelax * tempFourth ;
+ overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
+
+ 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)
+ + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
+ - overRelaxcoef5 * arrayV(i,j)
+ + coef1[i] * arrayV(i+iOne,j)
+ + sumChargeDensity(i,j)
+ ) * overRelaxtempFourth;
+ }
+ }
+
+ 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 rows
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
+ arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
+ arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
+ }
+
+ // 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++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
+ arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
+ arrayEz(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 i = 0 ; i < rows ; i++) {
+ // Note: go back and compare to old version of this code. See notes below.
+ // JT Test ... attempt to divide by real Ez not Ez to first order
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ arrayEz(i,j) += ezField;
+ // This adds back the overall Z gradient of the field (main E field component)
+ }
+ // Warning: (-=) assumes you are using an error potetial without the overall Field included
+ }
+
+ // 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 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+ arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
+ 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 ;
+ }
+ }
+
+}
+
+
+
+const Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) {
+ //
+ // Helperfunction: Check if integer is a power of 2
+ //
+ Int_t j = 0;
+ while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
+ if ( j == 1 ) return(1) ; // True
+ return(0) ; // False
+}
+
AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
//
// track1.fP[2] - sinus of local inclination angle
// track1.fP[3] - tangent of deep angle
// track1.fP[4] - 1/pt
+
AliTPCROC * roc = AliTPCROC::Instance();
const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
if (r0>80){
corr->DistortPoint(distPoint, sector);
}
- Double_t value=distPoint[2]-gz;
+ // Double_t value=distPoint[2]-gz;
if (itype==0){
Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
#include <TNamed.h>
+#include "TMatrixD.h"
class TH2F;
class TTimeStamp;
class TCollection;
class TTree;
class THnSparse;
+
+
class AliTPCCorrection : public TNamed {
public:
enum CompositionType {kParallel,kQueue};
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 const Int_t IsPowerOfTwo ( Int_t i ) ;
+
+ // 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 );
+
+
protected:
Double_t fT1; // tensor term of wt - T1
Double_t fT2; // tensor term of wt - T2
AliTPCGGVoltError::AliTPCGGVoltError()
: AliTPCCorrection("GGVoltError","GatingGrid (GG) Voltage Error"),
fC0(0.),fC1(0.),
- fDeltaVGGA(0.),fDeltaVGGC(0.)
+ fDeltaVGGA(0.),fDeltaVGGC(0.),
+ fInitLookUp(kFALSE)
{
//
// default constructor
Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
SetOmegaTauT1T2(wt,fT1,fT2);
- InitGGVoltErrorDistortion();
+ // InitGGVoltErrorDistortion(); // not necessary in here since the Voltage should not change!
}
// Electrostatic Equations from StarNote SN0253 by Howard Wieman.
//
+ if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
+
Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
Double_t intEr, intEphi ;
dx[1] = r * TMath::Sin(phi) - x[1];
dx[2] = 0.; // z distortion not implemented (1st order distortions)
+
+
}
// This function is purely for calibration purposes
// Calculates the integral (int Er/Ez dz) for the setted GG voltage offset
//
-
+
+ if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
+
Int_t order = 1 ; // FIXME: so far hardcoded? Linear interpolation = 1, Quadratic = 2
Double_t intEr;
}
}
+
+ fInitLookUp = kTRUE;
}
printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
}
+ if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitGGVoltErrorDistortion() ...");
}
Double_t fDeltaVGGC; // Missmatch of gating grid voltage on C-side [V]
Double_t fGGVoltErrorER[kNZ][kNR]; // Array to store electric field for GGVoltError calculation
+ Bool_t fInitLookUp; // flag to check it the Look Up table was created
+
ClassDef(AliTPCGGVoltError,1);
};
}
+void AliTPCInverseCorrection::Print(Option_t* option) const {
+ //
+ // Print function to check which correction classes are used
+ // option=="d" prints details regarding the setted magnitude
+ // option=="a" prints the C0 and C1 coefficents for calibration purposes
+ //
+
+ printf("Inverse of ");
+ if (fCorrection) fCorrection->Print(option);
+}
+
void AliTPCInverseCorrection::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
//
// This is just calling the CalculateInverseCorrection of the wrapped
for (Int_t j=0;j<3;++j) dx[j]=0.;
}
+void AliTPCInverseCorrection:: SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2) {
+ //
+ // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
+ // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
+ // calibration run
+ //
+ if (fCorrection) fCorrection->SetOmegaTauT1T2(omegaTau, t1, t2);
+}
+
void AliTPCInverseCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
//
// This is just calling the CalculateCorrection of the wrapped
virtual void Init();
virtual void Update(const TTimeStamp &timeStamp);
-
+ virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2);
+
+ // convenience functions
+ virtual void Print(Option_t* option="") const;
+
private:
AliTPCCorrection *fCorrection; // The correction to be inverted.
#pragma link C++ class AliTPCExBTwist+;
#pragma link C++ class AliTPCExBConical+;
#pragma link C++ class AliTPCGGVoltError+;
+#pragma link C++ class AliTPCBoundaryVoltError+;
#pragma link C++ class AliTPCCalibGlobalMisalignment+;
#pragma link C++ class AliTPCExBEffective+;
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 AliTPCExBConical.cxx AliXRDPROOFtoolkit.cxx AliTPCCalibGlobalMisalignment.cxx AliTPCExBEffective.cxx
+ AliTPCExBBShape.cxx AliTPCExBTwist.cxx AliTPCGGVoltError.cxx AliTPCBoundaryVoltError.cxx AliTPCExBConical.cxx AliXRDPROOFtoolkit.cxx AliTPCCalibGlobalMisalignment.cxx AliTPCExBEffective.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff