diff --git a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.cxx b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.cxx index 4a179ff2..66de4534 100644 --- a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.cxx +++ b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.cxx @@ -107,7 +107,7 @@ Double_t AliTPC3DCylindricalInterpolator::InterpolateCylindrical(Double_t r, Dou m = (kLow + k) % fNPhi; // interpolate for (Int_t i = iLow; i < iLow + fOrder + 1; i++) { - if (fOrder <= 2) { + if (fOrder < 3) { if (jLow >= 0) { index = m * (fNZ * fNR) + i * (fNZ) + jLow; saveArray[i - iLow] = Interpolate(&fZList[jLow], &fValue[index], z); @@ -176,19 +176,24 @@ Double_t AliTPC3DCylindricalInterpolator::InterpolatePhi( Int_t i2 = (iLow + 2) % lenX; Double_t xi2 = xArray[i2]; - if (xi1 < xi0) xi1 = TMath::TwoPi() + xi1; - if (xi2 < xi1) xi2 = TMath::TwoPi() + xi2; - if (x < xi0) x = TMath::TwoPi() + x; + if (fOrder <= 2) { + if (xi1 < xi0) xi1 = TMath::TwoPi() + xi1; + if (xi2 < xi1) xi2 = TMath::TwoPi() + xi2; + if (x < xi0) x = TMath::TwoPi() + x; + } Double_t y; if (fOrder > 2) { Double_t y2Array[fOrder + 1]; Double_t xArrayTemp[fOrder + 1]; Double_t dPhi = xArray[1] - xArray[0]; + // make list phi ascending order for (Int_t i = 0; i < fOrder + 1; i++) { xArrayTemp[i] = xArray[iLow] + (dPhi * i); } - if ((xArrayTemp[0] - x) > 1e-10) x = x + TMath::TwoPi(); + if (x < xArrayTemp[0]) x = TMath::TwoPi() + x; + if (x < xArrayTemp[0] || x > xArrayTemp[fOrder]) printf("x (%f) is outside of interpolation box (%f,%f)\n",x,xArrayTemp[0],xArrayTemp[fOrder]); + InitCubicSpline(xArrayTemp, yArray, fOrder + 1, y2Array, 1); y = InterpolateCubicSpline(xArrayTemp, yArray, y2Array, fOrder + 1, fOrder + 1, fOrder + 1, x, 1); } else if (fOrder == 2) { // Quadratic Interpolation = 2 @@ -422,6 +427,27 @@ void AliTPC3DCylindricalInterpolator::SetValue(TMatrixD **matricesVal) { } } +/// Set the value as interpolation point +/// +/// \param matricesVal TMatrixD** reference value for each point +void AliTPC3DCylindricalInterpolator::SetValue(TMatrixD **matricesVal,Int_t iZ) { + Int_t indexVal1D; + Int_t index1D; + TMatrixD *mat; + if (!fIsAllocatingLookUp) { + fValue = new Double_t[fNPhi * fNR * fNZ]; + fIsAllocatingLookUp = kTRUE; + } + for (Int_t m = 0; m < fNPhi; m++) { + indexVal1D = m * fNR * fNZ; + mat = matricesVal[m]; + for (Int_t i = 0; i < fNR; i++) { + index1D = indexVal1D + i * fNZ; + fValue[index1D + iZ] = (*mat)(i, iZ); + } + } +} + /// set the position of R /// /// \param rList diff --git a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.h b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.h index 8076cf9d..98ff6bcc 100644 --- a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.h +++ b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolator.h @@ -31,12 +31,14 @@ class AliTPC3DCylindricalInterpolator { void SetZList(Double_t *zList); void SetValue(Double_t *vList); void SetValue(TMatrixD **vList); + void SetValue(TMatrixD **vList,Int_t iZ); Int_t GetNR() { return fNR; } Int_t GetNPhi() { return fNPhi; } Int_t GetNZ() { return fNZ; } Int_t GetOrder() { return fOrder; } + Double_t * GetSecondDerZ() {return fSecondDerZ;} private: Int_t fOrder; ///< Order of interpolation, 1 - linear, 2 - quadratic, 3 >= - cubic, Int_t fNR; ///< Grid size in direction of R diff --git a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.cxx b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.cxx index 5a556676..f6f9aaf5 100644 --- a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.cxx +++ b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.cxx @@ -28,6 +28,7 @@ #include "TDecompSVD.h" #include "AliTPCPoissonSolver.h" #include "AliTPC3DCylindricalInterpolatorIrregular.h" +#include /// \cond CLASSIMP3 ClassImp(AliTPC3DCylindricalInterpolatorIrregular) @@ -87,6 +88,9 @@ AliTPC3DCylindricalInterpolatorIrregular::~AliTPC3DCylindricalInterpolatorIrregu delete fZList; } + if (fKDTreeIrregularPoints) { + delete[] fKDTreeIrregularPoints; + } delete[] fRBFWeightLookUp; delete[] fRBFWeight; } @@ -284,6 +288,7 @@ AliTPC3DCylindricalInterpolatorIrregular::Interpolate3DTableCylRBF( } } + index = new_phiIndex * (fNZ * fNR) + new_rIndex * fNZ + new_zIndex; phiStep = fStepPhi; rStep = fStepR; zStep = fStepZ; @@ -320,7 +325,7 @@ AliTPC3DCylindricalInterpolatorIrregular::Interpolate3DTableCylRBF( GetRBFWeightHalf(new_rIndex, new_zIndex, new_phiIndex, rStep, phiStep, zStep, radiusRBF0, 0, w); val = InterpRBFHalf(r, phi, z, startR, startPhi, startZ, rStep, phiStep, zStep, radiusRBF0, 0, w); } - delete w; + delete[] w; return val; } @@ -424,6 +429,20 @@ Double_t AliTPC3DCylindricalInterpolatorIrregular::GetValue( return Interpolate3DTableCylRBF(r, z, phi, rIndex, zIndex, phiIndex, rStep, phiStep, zStep, 0.0); } +// GetValue using searching at KDTree +Double_t AliTPC3DCylindricalInterpolatorIrregular::GetValue( + Double_t r, Double_t phi, Double_t z) { + + KDTreeNode n; + n.pR = &r; n.pPhi = φ n.pZ = &z; + KDTreeNode *nearest; + Double_t dist; + dist = 100000000.0; + Int_t startIndex = 0; // Z + Int_t dim = 3; // dimenstion + KDTreeNearest(fKDTreeIrregularRoot,&n, startIndex, dim, &nearest,&dist); + return Interpolate3DTableCylRBF(r,z,phi,nearest); +} /// Set value and distorted point for irregular grid interpolation /// /// \param matrixRicesValue @@ -467,6 +486,9 @@ void AliTPC3DCylindricalInterpolatorIrregular::SetValue( } } } + // KD Tree is used for look-up a point to irregular grid to find + // closest neughboor point + InitKDTree(); InitRBFWeight(); } @@ -488,9 +510,7 @@ void AliTPC3DCylindricalInterpolatorIrregular::InitRBFWeight() { Double_t radiusRBF0, minTemp, minTemp2; nd = fStepR * fStepPhi * fStepZ; - // init RBF weights for (Int_t m = 0; m < fNPhi; m++) { - indexInner = m * fNR * fNZ; for (Int_t i = 0; i < fNR; i++) { rIndex = indexInner + i * fNZ; @@ -831,7 +851,7 @@ void AliTPC3DCylindricalInterpolatorIrregular::GetRBFWeight( Int_t index = phiIndex * fNR * fNZ + rIndex * fNZ + zIndex; if (fRBFWeightLookUp[index] == 0) { RBFWeight(rIndex, zIndex, phiIndex, rStep, phiStep, zStep, radius0, kernelType, w); - \ + fRBFWeightLookUp[index] = 1; Int_t nd = rStep * zStep * phiStep; @@ -1162,3 +1182,204 @@ AliTPC3DCylindricalInterpolatorIrregular::GetRadius0RBF(const Int_t rIndex, cons return Distance(r0, 0.0, 0.0, r0 + gridSizeR, gridSizePhi, gridSizeR); } + + +// make kdtree for irregular look-up +void AliTPC3DCylindricalInterpolatorIrregular::InitKDTree() { + Int_t count = fNR * fNZ * fNPhi; + + fKDTreeIrregularPoints = new KDTreeNode[count]; + + for (Int_t i=0;ileft = MakeKDTree(t, (n - t), index, dim); + n->right = MakeKDTree(n + 1, (t + count) - (n + 1), index, dim); + } + return n; + +} + +// find median +AliTPC3DCylindricalInterpolatorIrregular::KDTreeNode * AliTPC3DCylindricalInterpolatorIrregular::FindMedian(KDTreeNode *start,KDTreeNode *end, Int_t index) { + if (end <= start) return NULL; + if (end == start + 1) + return start; + + + KDTreeNode *p, *store, *md = start + (end - start) / 2; + Double_t pivot; + + + while (1) { + if (index == 0) + pivot = *(md->pZ); + else if (index == 1) + pivot = *(md->pR); + else + pivot = *(md->pPhi); + + + Swap(md, end - 1); + + for (store = p = start; p < end; p++) { + + + if (((index ==0) && (*(p->pZ) < pivot)) || + ((index ==1) && (*(p->pR) < pivot)) || + ((index ==2) && (*(p->pPhi) < pivot))) + + { + if (p != store) + Swap(p, store); + store++; + } + + + } + Swap(store, end - 1); + + if ((index == 0) && (*(store->pZ) == *(md->pZ))) return md; + if ((index == 1) && (*(store->pR) == *(md->pR))) return md; + if ((index == 2) && (*(store->pPhi) == *(md->pPhi))) return md; + +// if (md->index == store->index) return md; + + if (store > md) end = store; + else start = store; + } + +} + +//swap +void AliTPC3DCylindricalInterpolatorIrregular::Swap(KDTreeNode *x, KDTreeNode *y) { + KDTreeNode *tmp = new KDTreeNode; + tmp->pR = x->pR; + tmp->pZ = x->pZ; + tmp->pPhi = x->pPhi; + tmp->index = x->index; + + x->pR = y->pR; + x->pZ = y->pZ; + x->pPhi = y->pPhi; + x->index = y->index; + + y->pR = tmp->pR; + y->pZ = tmp->pZ; + y->pPhi = tmp->pPhi; + y->index = tmp->index; + + delete tmp; +} + +// look for nearest point +void AliTPC3DCylindricalInterpolatorIrregular::KDTreeNearest(KDTreeNode *root, KDTreeNode *nd, Int_t index, Int_t dim, + KDTreeNode **best, Double_t *best_dist) +{ + Double_t d, dx2, dx; + + if (!root) return; + d = Distance(*(root->pR),*(root->pPhi),*(root->pZ),*(nd->pR),*(nd->pPhi),*(nd->pZ)); + if (index ==0) { + dx = *(root->pZ) - *(nd->pZ); + dx2 = Distance(*(nd->pR),*(nd->pPhi),*(root->pZ),*(nd->pR),*(nd->pPhi),*(nd->pZ)); + } else if (index == 1) { + dx = *(root->pR) - *(nd->pR); + dx2 = Distance(*(root->pR),*(nd->pPhi),*(nd->pZ),*(nd->pR),*(nd->pPhi),*(nd->pZ)); + } else { + dx = *(root->pPhi) - *(nd->pPhi); + dx2 = Distance(*(nd->pR),*(root->pPhi),*(nd->pZ),*(nd->pR),*(nd->pPhi),*(nd->pZ)); + } + + if (!*best || (d < *best_dist)) { + *best_dist = d; + *best = root; + } + + if (!*best_dist) return; + + if (++index >= dim) index = 0; + + KDTreeNearest(dx > 0 ? root->left : root->right, nd, index, dim, best, best_dist); + if (dx2 >= *best_dist) return; + KDTreeNearest(dx > 0 ? root->right : root->left, nd, index, dim, best, best_dist); + +} + +// interpolate on the nearest neighbor of irregular grid +Double_t +AliTPC3DCylindricalInterpolatorIrregular::Interpolate3DTableCylRBF( + Double_t r, Double_t z, Double_t phi, KDTreeNode * nearestNode) +{ + Double_t val = 0.0; + Int_t startPhi,startR,startZ; + Int_t phiIndex,rIndex,zIndex; + + + phiIndex = nearestNode->index / (fNR * fNZ); + rIndex = (nearestNode->index - (phiIndex * (fNR * fNZ)))/fNZ; + zIndex = nearestNode->index - (phiIndex * (fNR * fNZ)+ rIndex * fNZ); + + + startPhi =phiIndex - fStepPhi / 2; + startR = rIndex - fStepR / 2; + startZ = zIndex - fStepZ / 2; + + if (startPhi < 0) startPhi = fNPhi + startPhi; + + if (startR < 0) startR = 0; + if (startR + fStepR >= fNR) startR = fNR - fStepR; + + if (startZ < 0) startZ = 0; + if (startZ + fStepZ >= fNZ) startZ = fNZ - fStepZ; + Int_t indexPhi; + + Int_t index; + Double_t r0,z0,phi0; + + Int_t rStep = fStepR; + Int_t zStep = fStepZ; + Int_t phiStep = fStepPhi; + + Double_t *w; + + //Int_t nd = (phiStep-1) + (rStep-1) + (zStep-1) + 1; + Int_t nd = fStepPhi * fStepR * fStepZ; + + + + + + w = new Double_t[nd]; + + Float_t minTemp, minTemp2; + + Double_t radiusRBF0 = GetRadius0RBF(rIndex, phiIndex, zIndex); + + if (fType == 1) { + + for (Int_t i = 0; i < nd; i++) w[i] = 0.0; + GetRBFWeight(rIndex,zIndex,phiIndex, rStep, phiStep, zStep, radiusRBF0, 0, w); + val = InterpRBF(r, phi, z, startR, startPhi, startZ, rStep, phiStep, zStep, radiusRBF0, 0, w); + } else { + GetRBFWeightHalf(rIndex,zIndex, phiIndex, rStep, phiStep, zStep, radiusRBF0, 0, w); + val = InterpRBFHalf(r, phi, z, startR, startPhi, startZ, rStep, phiStep, zStep, radiusRBF0, 0, w); + } + delete[] w; + return val; +} + diff --git a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.h b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.h index 7cf3242f..c763e9cf 100644 --- a/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.h +++ b/TPCSpaceChargeBase/AliTPC3DCylindricalInterpolatorIrregular.h @@ -30,6 +30,7 @@ class AliTPC3DCylindricalInterpolatorIrregular { Double_t GetValue(Double_t r, Double_t phi, Double_t z, Int_t rIndex, Int_t phiIndex, Int_t zIndex, Int_t stepR, Int_t stepPhi, Int_t stepZ, Int_t minZColumnIndex); + Double_t GetValue(Double_t r, Double_t phi, Double_t z); void SetOrder(Int_t order) { fOrder = order; } void InitRBFWeight(); @@ -66,6 +67,14 @@ class AliTPC3DCylindricalInterpolatorIrregular { Int_t jy); private: + struct KDTreeNode { + Double_t *pR; + Double_t *pZ; + Double_t *pPhi; + Int_t index; + struct KDTreeNode *left, *right; + }; + Int_t fOrder; ///< Order of interpolation, 1 - linear, 2 - quadratic, 3 - cubic Int_t fType; ///< 0 INVERSE WEIGHT, 1 RBF FULL, 2 RBF Half Int_t fKernelType; ///< type kernel RBF 1--5 @@ -91,6 +100,7 @@ class AliTPC3DCylindricalInterpolatorIrregular { Int_t stepR, Int_t stepZ, Int_t stepPhi); Double_t Interpolate3DTableCylRBF(Double_t r, Double_t z, Double_t phi, Int_t rIndex, Int_t zIndex, Int_t phiIndex, Int_t stepR, Int_t stepZ, Int_t stepPhi, Double_t radiusRBF0); + Double_t Interpolate3DTableCylRBF(Double_t r, Double_t z, Double_t phi, KDTreeNode * nearestNode); void Search(Int_t n, const Double_t xArray[], Double_t x, Int_t &low); void Search(Int_t n, Double_t *xArray, Int_t offset, Double_t x, Int_t &low); @@ -116,6 +126,18 @@ class AliTPC3DCylindricalInterpolatorIrregular { Double_t radius0, Int_t kernelType, Double_t *weight); Double_t GetRadius0RBF(const Int_t rIndex, const Int_t phiIndex, const Int_t zIndex); + + KDTreeNode * fKDTreeIrregularPoints; // to save tree as list + KDTreeNode * fKDTreeIrregularRoot; // kdtree root + + void InitKDTree(); + KDTreeNode * MakeKDTree(KDTreeNode *tree,Int_t count, Int_t index, Int_t dimention); + + KDTreeNode * FindMedian(KDTreeNode *startTree,KDTreeNode *endTree, Int_t index); + void Swap(KDTreeNode *x, KDTreeNode *y); + + void KDTreeNearest(KDTreeNode *root, KDTreeNode *nd, Int_t index, Int_t dim, + KDTreeNode **best, Double_t *best_dist); /// \cond CLASSIMP ClassDef(AliTPC3DCylindricalInterpolatorIrregular,1); /// \endcond diff --git a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.cxx b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.cxx index 9da1ac46..0c2fa52b 100644 --- a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.cxx +++ b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.cxx @@ -44,6 +44,8 @@ AliTPCLookUpTable3DInterpolatorD::AliTPCLookUpTable3DInterpolatorD() { /// \param nZColumn Int_t size of grid Z direction /// \param zMin Double_t minimal value of Z /// \param zMax Double_t maximal value of Z + +/** AliTPCLookUpTable3DInterpolatorD::AliTPCLookUpTable3DInterpolatorD(Int_t nRRow, Double_t rMin, Double_t rMax, Int_t nPhiSlice, Double_t phiMin, Double_t phiMax, Int_t nZColumn, @@ -77,6 +79,8 @@ AliTPCLookUpTable3DInterpolatorD::AliTPCLookUpTable3DInterpolatorD(Int_t nRRow, for (Int_t m = 0; m < fNR; m++) fRList[m] = rMin + dR * m; for (Int_t m = 0; m < fNZ; m++) fZList[m] = zMin + dZ * m; } +**/ + /// Constructor /// @@ -105,12 +109,12 @@ AliTPCLookUpTable3DInterpolatorD::AliTPCLookUpTable3DInterpolatorD( SetNZ(nZColumn); SetLookUpZ(matricesZValue); SetZList(zList); - SetOrder(order); fInterpolatorR = new AliTPC3DCylindricalInterpolator(); fInterpolatorZ = new AliTPC3DCylindricalInterpolator(); fInterpolatorPhi = new AliTPC3DCylindricalInterpolator(); + SetOrder(order); fInterpolatorR->SetNR(nRRow); fInterpolatorR->SetNZ(nZColumn); fInterpolatorR->SetNPhi(nPhiSlice); @@ -173,6 +177,16 @@ void AliTPCLookUpTable3DInterpolatorD::CopyFromMatricesToInterpolator() { } +/// copy from matrices to 1D array for interpolation algorithm +void AliTPCLookUpTable3DInterpolatorD::CopyFromMatricesToInterpolator(Int_t iZ) { + fInterpolatorR->SetValue(fLookUpR,iZ); + fInterpolatorZ->SetValue(fLookUpZ,iZ); + fInterpolatorPhi->SetValue(fLookUpPhi,iZ); + + // no implementation for cubic spline interpolation +} + + /// get value of 3-components at a P(r,phi,z) /// /// \param r Double_t r position @@ -205,3 +219,13 @@ void AliTPCLookUpTable3DInterpolatorD::GetValue( zValue = fInterpolatorZ->GetValue(r, phi, z); } + +// Set Order of interpolation +// +void AliTPCLookUpTable3DInterpolatorD::SetOrder(Int_t order) { + fOrder = order; + fInterpolatorR->SetOrder(order); + fInterpolatorZ->SetOrder(order); + fInterpolatorPhi->SetOrder(order); +} + diff --git a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.h b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.h index 2f043475..95a42771 100644 --- a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.h +++ b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorD.h @@ -19,7 +19,7 @@ class AliTPCLookUpTable3DInterpolatorD { public: AliTPCLookUpTable3DInterpolatorD(); - AliTPCLookUpTable3DInterpolatorD(Int_t nRRow, Double_t rMin, Double_t rMax, Int_t nPhiSlice, Double_t phiMin, Double_t phiMax, Int_t nZColumn , Double_t zMin, Double_t zMax ); + //AliTPCLookUpTable3DInterpolatorD(Int_t nRRow, Double_t rMin, Double_t rMax, Int_t nPhiSlice, Double_t phiMin, Double_t phiMax, Int_t nZColumn , Double_t zMin, Double_t zMax ); AliTPCLookUpTable3DInterpolatorD(Int_t nRRow, TMatrixD**matricesRValue, Double_t *rList, Int_t nPhiSlice, TMatrixD**matricesPhiValue, Double_t *phiList, Int_t nZColumn, TMatrixD**matricesZValue, Double_t *zList, Int_t order); virtual ~AliTPCLookUpTable3DInterpolatorD(); @@ -38,11 +38,24 @@ class AliTPCLookUpTable3DInterpolatorD { void SetLookUpR(TMatrixD **matricesRValue) {fLookUpR = matricesRValue;} void SetLookUpPhi(TMatrixD **matricesPhiValue) {fLookUpPhi = matricesPhiValue;} void SetLookUpZ(TMatrixD **matricesZValue) {fLookUpZ = matricesZValue;} - void SetOrder(Int_t order) { fOrder = order; } + void SetOrder(Int_t order); void GetValue(Double_t r, Double_t phi, Double_t z, Double_t &rValue, Double_t &phiValue, Double_t &zValue); void GetValue(Double_t r, Double_t phi, Double_t z, Float_t &rValue, Float_t &phiValue, Float_t &zValue); void CopyFromMatricesToInterpolator(); + void CopyFromMatricesToInterpolator(Int_t iZ); // copy only iZ + TMatrixD ** GetLookUpR() { return fLookUpR;} + TMatrixD ** GetLookUpPhi() {return fLookUpPhi;} + TMatrixD ** GetLookUpZ() {return fLookUpZ;} + Double_t * GetRList() {return fRList;} + Double_t * GetPhiList() {return fPhiList;} + Double_t * GetZList() {return fZList;} + + + + AliTPC3DCylindricalInterpolator * GetInterpolatorR() {return fInterpolatorR;} + AliTPC3DCylindricalInterpolator * GetInterpolatorPhi() {return fInterpolatorPhi;} + AliTPC3DCylindricalInterpolator * GetInterpolatorZ() {return fInterpolatorZ;} private: Int_t fOrder; ///< order of interpolation Int_t fNR; ///< number of grid in R diff --git a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.cxx b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.cxx index 1d3f9b26..eda411fe 100644 --- a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.cxx +++ b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.cxx @@ -116,6 +116,7 @@ AliTPCLookUpTable3DInterpolatorIrregularD::~AliTPCLookUpTable3DInterpolatorIrreg /// copy from matrices to the interpolator void AliTPCLookUpTable3DInterpolatorIrregularD::CopyFromMatricesToInterpolator() { + fInterpolatorR->SetValue(fMatricesRValue, fMatricesRPoint, fMatricesPhiPoint, fMatricesZPoint); fInterpolatorZ->SetValue(fMatricesZValue, fMatricesRPoint, fMatricesPhiPoint, fMatricesZPoint); fInterpolatorPhi->SetValue(fMatricesPhiValue, fMatricesRPoint, fMatricesPhiPoint, fMatricesZPoint); @@ -196,3 +197,11 @@ void AliTPCLookUpTable3DInterpolatorIrregularD::GetValue( zValue = fInterpolatorZ->GetValue(r, phi, z, rIndex, phiIndex, zIndex, startR, startPhi, startZ); } +// using kdtree +void AliTPCLookUpTable3DInterpolatorIrregularD::GetValue( + Double_t r, Double_t phi, Double_t z, Double_t &rValue, Double_t &phiValue, Double_t &zValue) { + rValue = fInterpolatorR->GetValue(r, phi, z); + phiValue = fInterpolatorPhi->GetValue(r, phi, z); + zValue = fInterpolatorZ->GetValue(r, phi, z); +} + diff --git a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.h b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.h index 76161822..7fd59c18 100644 --- a/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.h +++ b/TPCSpaceChargeBase/AliTPCLookUpTable3DInterpolatorIrregularD.h @@ -46,6 +46,7 @@ class AliTPCLookUpTable3DInterpolatorIrregularD { Int_t phiIndex, Int_t zIndex, Int_t stepR, Int_t stepPhi, Int_t stepZ, Int_t minZColumnIndex); void GetValue(Double_t r, Double_t phi, Double_t z, Float_t &rValue, Float_t &phiValue, Float_t &zValue, Int_t rIndex, Int_t phiIndex, Int_t zIndex, Int_t stepR, Int_t stepPhi, Int_t stepZ); + void GetValue(Double_t r, Double_t phi, Double_t z, Double_t &rValue,Double_t &phiValue, Double_t &zValue); void SetOrder(Int_t order) { fOrder = order; } void CopyFromMatricesToInterpolator(); void CopyFromMatricesToInterpolator(Int_t j); diff --git a/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.cxx b/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.cxx index 3a5c97ee..223dedd8 100644 --- a/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.cxx +++ b/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.cxx @@ -17,7 +17,7 @@ /* $Id$ */ /// \class AliTPCSpaceCharge3DCalc -/// \brief This class provides distortion and correction map calculation with integration following electron drift +/// \brief This class provides distortion and correction map with integration following electron drift /// /// \author Rifki Sadikin , Indonesian Institute of Sciences /// \date Nov 20, 2017 @@ -41,20 +41,21 @@ ClassImp(AliTPCSpaceCharge3DCalc) AliTPCSpaceCharge3DCalc::AliTPCSpaceCharge3DCalc() : fC0(0.), fC1(0.), fCorrectionFactor(1.), fInitLookUp(kFALSE), fInterpolationOrder(5), fIrregularGridSize(3), fRBFKernelType(0), fNRRows(129), fNZColumns(129), fNPhiSlices(180), - fCorrectionType(0) { + fCorrectionType(1), fIntegrationStrategy(0) { InitAllocateMemory(); } -/// Construction for AliTPCSpaceCharge3DCalc class + /// Member values from params /// /// \param nRRow Int_t number of grid in r direction /// \param nZColumn Int_t number of grid in z direction /// \param nPhiSlice Int_t number of grid in \f$ \phi \f$ direction /// -AliTPCSpaceCharge3DCalc::AliTPCSpaceCharge3DCalc(Int_t nRRow, Int_t nZColumn, Int_t nPhiSlice) - : fC0(0.), fC1(0.), fCorrectionFactor(1.), fInitLookUp(kFALSE), - fInterpolationOrder(2), - fIrregularGridSize(3), fRBFKernelType(0), fCorrectionType(0) { +AliTPCSpaceCharge3DCalc::AliTPCSpaceCharge3DCalc( Int_t nRRow, + Int_t nZColumn, Int_t nPhiSlice) : + fC0(0.), fC1(0.), fCorrectionFactor(1.), fInitLookUp(kFALSE), + fInterpolationOrder(3), + fIrregularGridSize(3), fRBFKernelType(0), fCorrectionType(1), fIntegrationStrategy(0) { fNRRows = nRRow; fNPhiSlices = nPhiSlice; // the maximum of phi-slices so far = (8 per sector) fNZColumns = nZColumn; // the maximum on column-slices so ~ 2cm slicing @@ -75,7 +76,7 @@ AliTPCSpaceCharge3DCalc::AliTPCSpaceCharge3DCalc( Int_t nRRow, Int_t nZColumn, Int_t nPhiSlice, Int_t interpolationOrder, Int_t irregularGridSize, Int_t rbfKernelType) : fC0(0.), fC1(0.), fCorrectionFactor(1.), fInitLookUp(kFALSE), - fCorrectionType(0) { + fCorrectionType(1),fIntegrationStrategy(0) { fInterpolationOrder = interpolationOrder; fIrregularGridSize = irregularGridSize; @@ -89,7 +90,7 @@ AliTPCSpaceCharge3DCalc::AliTPCSpaceCharge3DCalc( /// Memory allocation for working/output memory /// void AliTPCSpaceCharge3DCalc::InitAllocateMemory() { - fPoissonSolver = new AliTPCPoissonSolver();; + fPoissonSolver = new AliTPCPoissonSolver(); fListR = new Double_t[fNRRows]; fListPhi = new Double_t[fNPhiSlices]; @@ -163,6 +164,8 @@ void AliTPCSpaceCharge3DCalc::InitAllocateMemory() { fMatrixChargeC[k] = new TMatrixD(fNRRows, fNZColumns); fMatrixChargeInverseA[k] = new TMatrixD(fNRRows, fNZColumns); fMatrixChargeInverseC[k] = new TMatrixD(fNRRows, fNZColumns); + fMatrixPotentialA[k] = new TMatrixD(fNRRows, fNZColumns); + fMatrixPotentialC[k] = new TMatrixD(fNRRows, fNZColumns); } fLookupIntDistA = @@ -292,20 +295,22 @@ void AliTPCSpaceCharge3DCalc::InitAllocateMemory() { fLookupIntCorrIrregularA->SetKernelType(fRBFKernelType); fLookupIntCorrIrregularC->SetKernelType(fRBFKernelType); - fFormulaBoundaryIFCA = NULL; - fFormulaBoundaryIFCC = NULL; - fFormulaBoundaryOFCA = NULL; - fFormulaBoundaryOFCC = NULL; - fFormulaBoundaryROCA = NULL; - fFormulaBoundaryROCC = NULL; - fFormulaBoundaryCE = NULL; + fFormulaBoundaryIFCA = NULL; //-> function define boundary values for IFC side A V(z) assuming symmetry in phi and r. + fFormulaBoundaryIFCC = NULL; //-> function define boundary values for IFC side C V(z) assuming symmetry in phi and r. + fFormulaBoundaryOFCA = NULL; //-> function define boundary values for OFC side A V(z) assuming symmetry in phi and r. + fFormulaBoundaryOFCC = NULL; ///<- function define boundary values for IFC side C V(z) assuming symmetry in phi and r. + fFormulaBoundaryROCA = NULL; ///<- function define boundary values for ROC side A V(r) assuming symmetry in phi and z. + fFormulaBoundaryROCC = NULL; ///<- function define boundary values for ROC side V V(t) assuming symmetry in phi and z. + fFormulaBoundaryCE = NULL; ///<- function define boundary values for CE V(z) assuming symmetry in phi and z. + + fFormulaPotentialV = NULL; ///<- potential V(r,rho,z) function + fFormulaChargeRho = NULL; ///<- charge density Rho(r,rho,z) function - fFormulaPotentialV = NULL; - fFormulaChargeRho = NULL; + // analytic formula for E + fFormulaEPhi = NULL; ///<- ePhi EPhi(r,rho,z) electric field (phi) function + fFormulaEr = NULL; ///<- er Er(r,rho,z) electric field (r) function + fFormulaEz = NULL; ///<- ez Ez(r,rho,z) electric field (z) function - fFormulaEPhi = NULL; - fFormulaEr = NULL; - fFormulaEz = NULL; } /// Destruction for AliTPCSpaceCharge3DCalc /// Deallocate memory for lookup table and charge distribution @@ -350,6 +355,9 @@ AliTPCSpaceCharge3DCalc::~AliTPCSpaceCharge3DCalc() { delete fMatrixChargeC[k]; delete fMatrixChargeInverseA[k]; delete fMatrixChargeInverseC[k]; + + delete fMatrixPotentialA[k]; + delete fMatrixPotentialC[k]; } delete[] fListR; delete[] fListPhi; @@ -381,6 +389,7 @@ AliTPCSpaceCharge3DCalc::~AliTPCSpaceCharge3DCalc() { delete[] fListPotentialBoundaryA; delete[] fListPotentialBoundaryC; } + /// Creating look-up tables of Correction/Distortion by integration following /// drift line, input from space charge 3d histogram (fSpaceCharge3D) and boundary values are filled with zeroes /// @@ -398,7 +407,7 @@ AliTPCSpaceCharge3DCalc::~AliTPCSpaceCharge3DCalc() { /// 2) Get the electric field \f$ \vec{E} = - \nabla \Phi(r,\phi,z) \f$ /// ~~~ /// ElectricField( matricesV, matricesEr, matricesEPhi, matricesEz, nRRow, nZColumn, phiSlice, -/// gridSizeR, gridSizePhi ,gridSizeZ,symmetry, AliTPCPoissonSolver::fgkIFCRadius); +/// gridSizeR, gridSizePhi ,gridSizeZ,symmetry, AliTPCPoissonSolver::fgkIFCRadius); /// ~~~ /// /// 3) Calculate local distortion and correction, using Langevin formula @@ -444,7 +453,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( // memory allocation for temporary matrices: // potential (boundary values), charge distribution - TMatrixD *matricesV[phiSlice], *matricesCharge[phiSlice]; + TMatrixD **matricesV, *matricesCharge[phiSlice]; TMatrixD *matricesEr[phiSlice], *matricesEPhi[phiSlice], *matricesEz[phiSlice]; TMatrixD *matricesDistDrDz[phiSlice], *matricesDistDPhiRDz[phiSlice], *matricesDistDz[phiSlice]; TMatrixD *matricesCorrDrDz[phiSlice], *matricesCorrDPhiRDz[phiSlice], *matricesCorrDz[phiSlice]; @@ -452,7 +461,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( TMatrixD *matricesGCorrDrDz[phiSlice], *matricesGCorrDPhiRDz[phiSlice], *matricesGCorrDz[phiSlice]; for (Int_t k = 0; k < phiSlice; k++) { - matricesV[k] = new TMatrixD(nRRow, nZColumn); + //matricesV[k] = new TMatrixD(nRRow, nZColumn); matricesCharge[k] = new TMatrixD(nRRow, nZColumn); matricesEr[k] = new TMatrixD(nRRow, nZColumn); matricesEPhi[k] = new TMatrixD(nRRow, nZColumn); @@ -527,6 +536,9 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( Double_t *potentialBoundary = NULL; TMatrixD *matrixV; TMatrixD *matrixCharge; + // for potential + TMatrixD **matricesVPotential; + Int_t pIndex = 0; // do if look up table haven't be initialized @@ -558,6 +570,9 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( matricesRIrregular = fMatrixRListIrregularA; matricesZIrregular = fMatrixZListIrregularA; matricesLookUpCharge = fMatrixChargeA; + + + matricesV = fMatrixPotentialA; chargeInterpolator = fInterpolatorChargeA; potentialInterpolator = fInterpolatorPotentialA; fLookupDistA->SetLookUpR(matricesDistDrDz); @@ -581,6 +596,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( matricesRIrregular = fMatrixRListIrregularC; matricesZIrregular = fMatrixZListIrregularC; matricesLookUpCharge = fMatrixChargeC; + matricesV = fMatrixPotentialC; chargeInterpolator = fInterpolatorChargeC; potentialInterpolator = fInterpolatorPotentialC; fLookupDistC->SetLookUpR(matricesDistDrDz); @@ -601,7 +617,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( // fill also charge //pIndex = 0; - //Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz",Step = 0: Fill Boundary and Charge Densities"); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step = 0: Fill Boundary and Charge Densities")); for (Int_t k = 0; k < phiSlice; k++) { phi0 = k * gridSizePhi; matrixV = matricesV[k]; @@ -612,7 +628,8 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( z0 = j * gridSizeZ; (*matrixCharge)(i, j) = chargeInterpolator->GetValue(rList[i], phiList[k], zList[j]); (*matrixV)(i, j) = 0.0; // fill zeros - if (fFormulaPotentialV == NULL) { + + if (fFormulaPotentialV == NULL) { // boundary IFC if (i == 0) { if (f1BoundaryIFC != NULL) { @@ -643,11 +660,11 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 0: Preparing Charge interpolator: %f\n", w.CpuTime())); AliTPCPoissonSolver::fgConvergenceError = stoppingConvergence; - fPoissonSolver->SetStrategy(AliTPCPoissonSolver::kMultiGrid); - (fPoissonSolver->fMgParameters).cycleType = AliTPCPoissonSolver::kFCycle; - (fPoissonSolver->fMgParameters).isFull3D = kFALSE; - (fPoissonSolver->fMgParameters).nMGCycle = maxIteration; - (fPoissonSolver->fMgParameters).maxLoop = 6; + //fPoissonSolver->SetStrategy(AliTPCPoissonSolver::kMultiGrid); + //(fPoissonSolver->fMgParameters).cycleType = AliTPCPoissonSolver::kFCycle; + //(fPoissonSolver->fMgParameters).isFull3D = kFALSE; + //(fPoissonSolver->fMgParameters).nMGCycle = maxIteration; + //(fPoissonSolver->fMgParameters).maxLoop = 6; w.Start(); @@ -660,12 +677,17 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 1: Poisson solver: %f\n", w.CpuTime())); + if (side == 0) myProfile.poissonSolverTime = w.CpuTime(); + if (side == 0) myProfile.iteration = fPoissonSolver->fIterations; + + w.Start(); ElectricField(matricesV, matricesEr, matricesEPhi, matricesEz, nRRow, nZColumn, phiSlice, gridSizeR, gridSizePhi, gridSizeZ, symmetry, AliTPCPoissonSolver::fgkIFCRadius); w.Stop(); + myProfile.electricFieldTime = w.CpuTime(); Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 2: Electric Field Calculation: %f\n", w.CpuTime())); w.Start(); LocalDistCorrDz(matricesEr, matricesEPhi, matricesEz, @@ -673,6 +695,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( matricesCorrDrDz, matricesCorrDPhiRDz, matricesCorrDz, nRRow, nZColumn, phiSlice, gridSizeZ, ezField); w.Stop(); + myProfile.localDistortionTime = w.CpuTime(); // copy to interpolator if (side == 0) { @@ -689,36 +712,47 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 3: Local distortion and correction: %f\n", w.CpuTime())); w.Start(); - - IntegrateDistCorrDriftLineDz( - lookupLocalDist, - matricesGDistDrDz, matricesGDistDPhiRDz, matricesGDistDz, - lookupLocalCorr, - matricesGCorrDrDz, matricesGCorrDPhiRDz, matricesGCorrDz, - matricesIrregularDrDz, matricesIrregularDPhiRDz, matricesIrregularDz, - matricesRIrregular, matricesPhiIrregular, matricesZIrregular, - nRRow, nZColumn, phiSlice, rList, phiList, zList - ); + if (fIntegrationStrategy == kNaive) + IntegrateDistCorrDriftLineDz( + lookupLocalDist, + matricesGDistDrDz, matricesGDistDPhiRDz, matricesGDistDz, + lookupLocalCorr, + matricesGCorrDrDz, matricesGCorrDPhiRDz, matricesGCorrDz, + matricesIrregularDrDz, matricesIrregularDPhiRDz, matricesIrregularDz, + matricesRIrregular, matricesPhiIrregular, matricesZIrregular, + nRRow, nZColumn, phiSlice, rList, phiList, zList + ); + else + IntegrateDistCorrDriftLineDzWithLookUp ( + lookupLocalDist, + matricesGDistDrDz, matricesGDistDPhiRDz, matricesGDistDz, + lookupLocalCorr, + matricesGCorrDrDz, matricesGCorrDPhiRDz, matricesGCorrDz, + nRRow, nZColumn, phiSlice, rList, phiList, zList + ); w.Stop(); Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 4: Global correction/distortion: %f\n", w.CpuTime())); - w.Start(); + myProfile.globalDistortionTime = w.CpuTime(); //// copy to 1D interpolator ///// lookupGlobalDist->CopyFromMatricesToInterpolator(); - lookupGlobalCorr->CopyFromMatricesToInterpolator(); + if (fCorrectionType==0) + lookupGlobalCorr->CopyFromMatricesToInterpolator(); //// w.Stop(); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 5: Filling up the look up: %f\n", w.CpuTime())); if (side == 0) { + + w.Start(); FillLookUpTable(lookupGlobalDist, fMatrixIntDistDrEzA, fMatrixIntDistDPhiREzA, fMatrixIntDistDzA, nRRow, nZColumn, phiSlice, rList, phiList, zList); - FillLookUpTable(lookupGlobalCorr, + if (fCorrectionType ==0) + FillLookUpTable(lookupGlobalCorr, fMatrixIntCorrDrEzA, fMatrixIntCorrDPhiREzA, fMatrixIntCorrDzA, nRRow, nZColumn, phiSlice, rList, phiList, zList); @@ -728,14 +762,19 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( else fLookupIntCorrIrregularA->CopyFromMatricesToInterpolator(); + w.Stop(); + myProfile.interpolationInitTime = w.CpuTime(); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 5: Filling up the look up: %f\n", w.CpuTime())); Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz"," A side done"); } if (side == 1) { + w.Start(); FillLookUpTable(lookupGlobalDist, fMatrixIntDistDrEzC, fMatrixIntDistDPhiREzC, fMatrixIntDistDzC, nRRow, nZColumn, phiSlice, rList, phiList, zList); - FillLookUpTable(lookupGlobalCorr, + if (fCorrectionType == 0) + FillLookUpTable(lookupGlobalCorr, fMatrixIntCorrDrEzC, fMatrixIntCorrDPhiREzC, fMatrixIntCorrDzC, nRRow, nZColumn, phiSlice, rList, phiList, zList); @@ -744,6 +783,9 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( fLookupIntCorrC->CopyFromMatricesToInterpolator(); else fLookupIntCorrIrregularC->CopyFromMatricesToInterpolator(); + w.Stop(); + myProfile.interpolationInitTime = w.CpuTime(); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 5: Filling up the look up: %f\n", w.CpuTime())); Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz"," C side done"); } @@ -756,7 +798,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( // memory de-allocation for temporary matrices for (Int_t k = 0; k < phiSlice; k++) { - delete matricesV[k]; + //delete matricesV[k]; delete matricesCharge[k]; delete matricesEr[k]; delete matricesEPhi[k]; @@ -822,7 +864,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( TMatrixD **matricesCorrDrDzA, TMatrixD **matricesCorrDPhiRDzA, TMatrixD **matricesCorrDzA, TMatrixD **matricesDistDrDzC, TMatrixD **matricesDistDPhiRDzC, TMatrixD **matricesDistDzC, TMatrixD **matricesCorrDrDzC, TMatrixD **matricesCorrDPhiRDzC, TMatrixD **matricesCorrDzC, - TFormula *intErDzTestFunction, TFormula *intEPhiRDzTestFunction, TFormula *intDzTestFunction) { + TFormula *intErDzTestFunction, TFormula *intEPhiRDzTestFunction, TFormula *intDzTestFunction, TFormula *ezFunction) { Int_t phiSlicesPerSector = phiSlice / kNumSector; const Float_t gridSizeR = (AliTPCPoissonSolver::fgkOFCRadius - AliTPCPoissonSolver::fgkIFCRadius) / (nRRow - 1); const Float_t gridSizeZ = AliTPCPoissonSolver::fgkTPCZ0 / (nZColumn - 1); @@ -1006,7 +1048,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( // fill also charge //pIndex = 0; - //Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","Step = 0: Fill Boundary and Charge Densities"); + //Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step = 0: Fill Boundary and Charge Densities")); for (Int_t k = 0; k < phiSlice; k++) { phi0 = k * gridSizePhi; matrixV = matricesV[k]; @@ -1071,7 +1113,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( w.Start(); - IntegrateDistCorrDriftLineDz(intErDzTestFunction,intEPhiRDzTestFunction, intDzTestFunction, ezField, + IntegrateDistCorrDriftLineDz(intErDzTestFunction,intEPhiRDzTestFunction, intDzTestFunction, ezFunction,ezField, matricesGDistDrDz, matricesGDistDPhiRDz, matricesGDistDz, matricesGCorrDrDz, matricesGCorrDPhiRDz, matricesGCorrDz, matricesIrregularDrDz, matricesIrregularDPhiRDz, matricesIrregularDz, @@ -1079,7 +1121,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( nRRow, nZColumn, phiSlice, rList, phiList, zList); w.Stop(); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","Step 4: Global correction/distortion: %f\n", w.CpuTime()); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 4: Global correction/distortion: %f\n", w.CpuTime())); w.Start(); //// copy to 1D interpolator ///// @@ -1089,7 +1131,7 @@ void AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz( w.Stop(); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","Step 5: Filling up the look up: %f\n", w.CpuTime()); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Step 5: Filling up the look up: %f\n", w.CpuTime())); if (side == 0) { FillLookUpTable(lookupGlobalDist, @@ -1249,10 +1291,10 @@ AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson(Int_t nRRow, Int_t nZColumn, I fInterpolationOrder ); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson","Step 1: Solving poisson solver"); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","Step 1: Solving poisson solver"); fPoissonSolver->PoissonSolver3D(matricesV, matricesCharge, nRRow, nZColumn, phiSlice, maxIteration, symmetry); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson","Step 2: Calculate electric field"); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","Step 2: Calculate electric field"); CalculateEField( matricesV, matricesEr, @@ -1265,14 +1307,14 @@ AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson(Int_t nRRow, Int_t nZColumn, I symmetry ); lookupEField->CopyFromMatricesToInterpolator(); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson","Step 3: Fill the look up table"); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","Step 3: Fill the look up table"); if (side == 0) { FillLookUpTable(lookupEField, fMatrixErOverEzA, fMatrixEPhiOverEzA, fMatrixDeltaEzA, nRRow, nZColumn, phiSlice, rList, phiList, zList); fLookupIntENoDriftA->CopyFromMatricesToInterpolator(); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson"," A side done"); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz"," A side done"); } if (side == 1) { FillLookUpTable(lookupEField, @@ -1280,7 +1322,7 @@ AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson(Int_t nRRow, Int_t nZColumn, I nRRow, nZColumn, phiSlice, rList, phiList, zList); fLookupIntENoDriftC->CopyFromMatricesToInterpolator(); - Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoisson"," C side done"); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz"," C side done"); } delete lookupEField; } @@ -1313,17 +1355,17 @@ void AliTPCSpaceCharge3DCalc::ForceInitSpaceCharge3DPoissonIntegralDz /// Electric field Calculation: /// /// -/// \param matricesV -/// \param matricesEr -/// \param matricesEPhi -/// \param matricesEz -/// \param nRRow -/// \param nZColumn -/// \param phiSlice -/// \param gridSizeR -/// \param gridSizePhi -/// \param gridSizeZ -/// \param symmetry +/// \param matricesV +/// \param matricesEr +/// \param matricesEPhi +/// \param matricesEz +/// \param nRRow +/// \param nZColumn +/// \param phiSlice +/// \param gridSizeR +/// \param gridSizePhi +/// \param gridSizeZ +/// \param symmetry /// \param innerRadius /// /// \pre Matrix matricesV is assumed had been calculated by Poisson solver @@ -1350,7 +1392,7 @@ void AliTPCSpaceCharge3DCalc::ForceInitSpaceCharge3DPoissonIntegralDz /// \f$ -\nabla_{r} V(r_{0},\phi_{j},z_{k}) \approx -( -0.5 V_{2,j,k} + 2 V_{1,j,k} - 1.5 * V_{0,j,k}) / h_{r} \f$ /// /// \f$ -\nabla_{r} V(r_{nRRow - 1},\phi_{j},z_{k}) \approx -( 1.5 V_{nRRow-1,j,k} - 2.0 V_{nRRow-2,j,k} + 0.5 V_{nRRow -3,j,k}) / h_{\phi} \f$ -/// +/// void AliTPCSpaceCharge3DCalc::ElectricField(TMatrixD **matricesV, TMatrixD **matricesEr, TMatrixD **matricesEPhi, TMatrixD **matricesEz, const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, @@ -1479,8 +1521,8 @@ void AliTPCSpaceCharge3DCalc::ElectricField(TMatrixD **matricesV, TMatrixD **mat /// ~~~ /// Double_t ezField = (AliTPCPoissonSolver::fgkCathodeV-AliTPCPoissonSolver::fgkGG)/AliTPCPoissonSolver::fgkTPCZ0; // = Electric Field (V/cm) Magnitude ~ -400 V/cm; /// -/// localIntErOverEz = (gridSizeZ/2.0)*((*eR)(i,j)+(*eR)(i,j+1))/(-1*ezField) ; -/// localIntEPhiOverEz = (gridSizeZ/2.0)*((*ePhi)(i,j)+(*ePhi)(i,j+1))/(-1*ezField) ; +/// localIntErOverEz = (gridSizeZ/2.0)*((*eR)(i,j)+(*eR)(i,j+1))/(ezField + (*eZ)(i,j)) ; +/// localIntEPhiOverEz = (gridSizeZ/2.0)*((*ePhi)(i,j)+(*ePhi)(i,j+1))/(ezField + (*eZ)(i,j)) ; /// localIntDeltaEz = (gridSizeZ/2.0)*((*eZ)(i,j)+(*eZ)(i,j+1)) ; /// ~~~ /// @@ -1503,7 +1545,7 @@ void AliTPCSpaceCharge3DCalc::ElectricField(TMatrixD **matricesV, TMatrixD **mat /// \f$ \hat{\delta}_{z}(r_{i},z_{j},\phi_{m}) = \int_{z_{j}}^{z_{j+1}} \frac{v^{\prime}(E)}{v_{0}} (E - E_{0}) dzDist\f$ /// /// ~~~ -/// (*distDz)(i,j) = localIntDeltaEz*AliTPCPoissonSolver::fgkdvdE; +/// (*distDz)(i,j) = localIntDeltaEz*-1*AliTPCPoissonSolver::fgkdvdE; /// ~~~ /// /// Where \f$c_{0}\f$ and \f$c_{1}\f$ are constants (see the ALICE-INT-2010-016 for further details). @@ -1633,14 +1675,13 @@ AliTPCSpaceCharge3DCalc::LocalDistCorrDz(TMatrixD **matricesEr, TMatrixD **matri for (Int_t j = 0; j < nZColumn - 1; j++) { for (Int_t i = 0; i < nRRow; i++) { - localIntErOverEz = (gridSizeZ / 2.0) * ((*eR)(i, j) + (*eR)(i, j + 1)) / (-1 * ezField); - localIntEPhiOverEz = (gridSizeZ / 2.0) * ((*ePhi)(i, j) + (*ePhi)(i, j + 1)) / (-1 * ezField); - localIntDeltaEz = (gridSizeZ / 2.0) * ((*eZ)(i, j) + (*eZ)(i, j + 1)); + localIntErOverEz = (gridSizeZ * 0.5) * ((*eR)(i, j) + (*eR)(i, j + 1)) / (ezField + (*eZ)(i,j)); + localIntEPhiOverEz = (gridSizeZ * 0.5) * ((*ePhi)(i, j) + (*ePhi)(i, j + 1)) / (ezField + (*eZ)(i,j)); + localIntDeltaEz = (gridSizeZ * 0.5) * ((*eZ)(i, j) + (*eZ)(i, j + 1)); (*distDrDz)(i, j) = fC0 * localIntErOverEz + fC1 * localIntEPhiOverEz; (*distDPhiRDz)(i, j) = fC0 * localIntEPhiOverEz - fC1 * localIntErOverEz; - (*distDz)(i, j) = localIntDeltaEz * AliTPCPoissonSolver::fgkdvdE * AliTPCPoissonSolver::fgkdvdE;// two times? - + (*distDz)(i, j) = localIntDeltaEz * -1 * AliTPCPoissonSolver::fgkdvdE; (*corrDrDz)(i, j + 1) = -1 * (*distDrDz)(i, j); (*corrDPhiRDz)(i, j + 1) = -1 * (*distDPhiRDz)(i, j); @@ -1687,7 +1728,7 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, const Double_t *rList, const Double_t *phiList, const Double_t *zList) { - Float_t drDist, dRPhi, dzDist, ddR, ddRPhi, ddZ; + Float_t drDist, dPhi, dzDist, ddR, ddRPhi, ddZ; Float_t radius0, phi0, z0, radius, phi, z, radiusCorrection; radiusCorrection = 0.0; radius = 0.0; @@ -1731,29 +1772,24 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( // follow the drift radius0 = rList[i]; phi = phi0; - radius = radius0; - - drDist = 0.0; - dRPhi = 0.0; - dzDist = 0.0; - ddRPhi = 0.0; /// - (*mDistDrDz)(i, j) = drDist; - (*mDistDPhiRDz)(i, j) = dRPhi; - (*mDistDz)(i, j) = dzDist; + (*mDistDrDz)(i, j) = 0.; + (*mDistDPhiRDz)(i, j) = 0.; + (*mDistDz)(i, j) = 0.; //////////////// use irregular grid look up table for correction - // set - (*mCorrIrregularDrDz)(i, j) = -drDist; - (*mCorrIrregularDPhiRDz)(i, j) = -dRPhi; - (*mCorrIrregularDz)(i, j) = -dzDist; + if (fCorrectionType == kIrregularInterpolator) { + (*mCorrIrregularDrDz)(i, j) = 0.0; + (*mCorrIrregularDPhiRDz)(i, j) = 0.0; + (*mCorrIrregularDz)(i, j) = -0.0; - // distorted point - (*mRIrregular)(i, j) = radius0 + drDist; - (*mPhiIrregular)(i, j) = phi0 + (dRPhi / radius0); - (*mZIrregular)(i, j) = z0 + dzDist; + // distorted point + (*mRIrregular)(i, j) = radius0; + (*mPhiIrregular)(i, j) = phi0; + (*mZIrregular)(i, j) = z0; + } /////////////// } @@ -1791,15 +1827,15 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( radius = radius0; drDist = 0.0; - dRPhi = 0.0; + dPhi = 0.0; dzDist = 0.0; - ddRPhi = 0.0; // follow the drift line from z=j --> nZColumn - 1 for (Int_t jj = j; jj < nZColumn; jj++) { // interpolation the local distortion for current position - phi += ddRPhi / radius; + // phi += ddRPhi / radius; + phi = phi0 + dPhi; radius = radius0 + drDist; z = zList[jj] + dzDist; @@ -1807,94 +1843,100 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( while (phi < 0.0) phi = TMath::TwoPi() + phi; while (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); + lookupLocalDist->GetValue(radius, phi, z, ddR, ddRPhi, ddZ); // add local distortion drDist += ddR; - dRPhi += ddRPhi; + dPhi += (ddRPhi/radius); dzDist += ddZ; } // set the global distortion after following the electron drift (*mDistDrDz)(i, j) = drDist; - (*mDistDPhiRDz)(i, j) = dRPhi; + (*mDistDPhiRDz)(i, j) = dPhi * radius0; (*mDistDz)(i, j) = dzDist; /////////////// use irregular grid look up table for correction // set - (*mCorrIrregularDrDz)(i, j) = -drDist; - (*mCorrIrregularDPhiRDz)(i, j) = -dRPhi; - (*mCorrIrregularDz)(i, j) = -dzDist; + if (fCorrectionType == kIrregularInterpolator) { + (*mCorrIrregularDrDz)(i, j) = - drDist; + (*mCorrIrregularDPhiRDz)(i, j) = -1 * dPhi * radius0; + (*mCorrIrregularDz)(i, j) = -dzDist; - // distorted point - (*mRIrregular)(i, j) = radius0 + drDist; - (*mPhiIrregular)(i, j) = phi0 + (dRPhi / radius0); - (*mZIrregular)(i, j) = z0 + dzDist; + // distorted point + (*mRIrregular)(i, j) = radius0 + drDist; + (*mPhiIrregular)(i, j) =phi0 + dPhi;; + (*mZIrregular)(i, j) = z0 + dzDist; + } /////////////// // put the radius to the original value - if (j == nZColumn - 2) radiusCorrection = radius0; + if (fCorrectionType == kRegularInterpolator) { + if (j == nZColumn - 2) radiusCorrection = radius0; - // get global correction from j+1 - drDist = (*mCorrDrDz)(i, j + 1); - dRPhi = (*mCorrDPhiRDz)(i, j + 1); - dzDist = (*mCorrDz)(i, j + 1); + // get global correction from j+1 + drDist = (*mCorrDrDz)(i, j + 1); + dPhi = (*mCorrDPhiRDz)(i, j + 1) / radius0; + dzDist = (*mCorrDz)(i, j + 1); - radiusCorrection = radius0 + drDist; - phi = phi0 + dRPhi / radiusCorrection; - z = zList[j + 1] + dzDist; + radiusCorrection = radius0 + drDist; + phi = phi0 + dPhi; + z = zList[j + 1] + dzDist; - while (phi < 0.0) phi = TMath::TwoPi() + phi; - while (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); + while (phi < 0.0) phi = TMath::TwoPi() + phi; + while (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); - lookupLocalCorr->GetValue(radiusCorrection, phi, z, ddR, ddRPhi, ddZ); + lookupLocalCorr->GetValue(radiusCorrection, phi, z, ddR, ddRPhi, ddZ); - drDist += ddR; - dzDist += ddZ; - dRPhi += ddRPhi; + drDist += ddR; + dzDist += ddZ; + dPhi += ddRPhi / radiusCorrection; - (*mCorrDrDz)(i, j) = drDist; - (*mCorrDPhiRDz)(i, j) = dRPhi; - (*mCorrDz)(i, j) = dzDist; + (*mCorrDrDz)(i, j) = drDist; + (*mCorrDPhiRDz)(i, j) = dPhi * radius0; + (*mCorrDz)(i, j) = dzDist; + } } } } } + /// follow the drift for exact function -/// -/// \param intDrDzF -/// \param intDPhiDzF -/// \param intDzDzF -/// \param ezField -/// \param matricesGDistDrDz -/// \param matricesGDistDPhiRDz -/// \param matricesGDistDz -/// \param matricesGCorrDrDz -/// \param matricesGCorrDPhiRDz -/// \param matricesGCorrDz -/// \param matricesGCorrIrregularDrDz -/// \param matricesGCorrIrregularDPhiRDz -/// \param matricesGCorrIrregularDz -/// \param matricesRIrregular -/// \param matricesPhiIrregular -/// \param matricesZIrregular -/// \param nRRow -/// \param nZColumn -/// \param phiSlice -/// \param rList -/// \param phiList -/// \param zList +/// +/// \param intDrDzF +/// \param intDPhiDzF +/// \param intDzDzF +/// \param ezField +/// \param matricesGDistDrDz +/// \param matricesGDistDPhiRDz +/// \param matricesGDistDz +/// \param matricesGCorrDrDz +/// \param matricesGCorrDPhiRDz +/// \param matricesGCorrDz +/// \param matricesGCorrIrregularDrDz +/// \param matricesGCorrIrregularDPhiRDz +/// \param matricesGCorrIrregularDz +/// \param matricesRIrregular +/// \param matricesPhiIrregular +/// \param matricesZIrregular +/// \param nRRow +/// \param nZColumn +/// \param phiSlice +/// \param rList +/// \param phiList +/// \param zList void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( - TFormula *intDrDzF, TFormula *intDPhiDzF, TFormula *intDzDzF, const Double_t ezField, + TFormula *intDrDzF, TFormula *intDPhiDzF, TFormula *intDzDzF, TFormula *ezF, const Double_t ezField, TMatrixD **matricesGDistDrDz, TMatrixD **matricesGDistDPhiRDz, TMatrixD **matricesGDistDz, TMatrixD **matricesGCorrDrDz, TMatrixD **matricesGCorrDPhiRDz, TMatrixD **matricesGCorrDz, TMatrixD **matricesGCorrIrregularDrDz, TMatrixD **matricesGCorrIrregularDPhiRDz, TMatrixD **matricesGCorrIrregularDz, TMatrixD **matricesRIrregular, TMatrixD **matricesPhiIrregular, TMatrixD **matricesZIrregular, const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, const Double_t *rList, const Double_t *phiList, const Double_t *zList) { - Float_t drDist, dRPhi, dzDist, ddR, ddRPhi, ddZ; + Float_t drDist, dPhi, dzDist, ddR, ddRPhi, ddZ; Float_t radius0, phi0, z0, radius, phi, z, radiusCorrection,z1; Float_t localIntErOverEz = 0.0; @@ -1939,33 +1981,26 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( mZIrregular = matricesZIrregular[m]; for (Int_t i = 0; i < nRRow; i++) { - // do from j to 0 // follow the drift radius0 = rList[i]; phi = phi0; radius = radius0; - - drDist = 0.0; - dRPhi = 0.0; - dzDist = 0.0; - ddRPhi = 0.0; - /// - (*mDistDrDz)(i, j) = drDist; - (*mDistDPhiRDz)(i, j) = dRPhi; - (*mDistDz)(i, j) = dzDist; + (*mDistDrDz)(i, j) = 0.0; + (*mDistDPhiRDz)(i, j) =0.0; + (*mDistDz)(i, j) = 0.0; //////////////// use irregular grid look up table for correction // set - (*mCorrIrregularDrDz)(i, j) = -drDist; - (*mCorrIrregularDPhiRDz)(i, j) = -dRPhi; - (*mCorrIrregularDz)(i, j) = -dzDist; + (*mCorrIrregularDrDz)(i, j) = 0.0; + (*mCorrIrregularDPhiRDz)(i, j) = 0.0; + (*mCorrIrregularDz)(i, j) = 0.0; // distorted point - (*mRIrregular)(i, j) = radius0 + drDist; - (*mPhiIrregular)(i, j) = phi0 + (dRPhi / radius0); - (*mZIrregular)(i, j) = z0 + dzDist; + (*mRIrregular)(i, j) = radius0; + (*mPhiIrregular)(i, j) = phi0 ; + (*mZIrregular)(i, j) = z0; /////////////// } @@ -2004,7 +2039,7 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( radius = radius0; drDist = 0.0; - dRPhi = 0.0; + dPhi = 0.0; dzDist = 0.0; ddRPhi = 0.0; @@ -2012,46 +2047,44 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( // follow the drift line from z=j --> nZColumn - 1 for (Int_t jj = j; jj < nZColumn; jj++) { // interpolation the local distortion for current position - phi += ddRPhi / radius; + phi = phi0 + dPhi; radius = radius0 + drDist; z = zList[jj] + dzDist; z1 = z + (zList[j+1] - zList[j]); - // regulate phi while (phi < 0.0) phi = TMath::TwoPi() + phi; while (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); //lookupLocalDist->GetValue(radius, phi, z, ddR, ddRPhi, ddZ); - localIntErOverEz = (intDrDzF->Eval(radius, phi, z1) - intDrDzF->Eval(radius, phi, z)) / (-1 * ezField); - localIntEPhiOverEz = (intDPhiDzF->Eval(radius, phi, z1) - intDPhiDzF->Eval(radius, phi, z)) / (-1 * ezField); + localIntErOverEz = (intDrDzF->Eval(radius, phi, z1) - intDrDzF->Eval(radius, phi, z)) / (ezField + ezF->Eval(radius,phi,z)); + localIntEPhiOverEz = (intDPhiDzF->Eval(radius, phi, z1) - intDPhiDzF->Eval(radius, phi, z)) / (ezField + ezF->Eval(radius,phi,z)); localIntDeltaEz = intDzDzF->Eval(radius, phi, z1) - intDzDzF->Eval(radius, phi, z); ddR = fC0 * localIntErOverEz + fC1 * localIntEPhiOverEz; ddRPhi = fC0 * localIntEPhiOverEz - fC1 * localIntErOverEz; - ddZ = localIntDeltaEz * AliTPCPoissonSolver::fgkdvdE * AliTPCPoissonSolver::fgkdvdE; // two times? - + ddZ = -1 * localIntDeltaEz * AliTPCPoissonSolver::fgkdvdE; - - // add local distortion drDist += ddR; - dRPhi += ddRPhi; + dPhi += (ddRPhi/radius); dzDist += ddZ; + // add local distortion + } // set the global distortion after following the electron drift (*mDistDrDz)(i, j) = drDist; - (*mDistDPhiRDz)(i, j) = dRPhi; + (*mDistDPhiRDz)(i, j) = dPhi *radius0; (*mDistDz)(i, j) = dzDist; /////////////// use irregular grid look up table for correction // set (*mCorrIrregularDrDz)(i, j) = -drDist; - (*mCorrIrregularDPhiRDz)(i, j) = -dRPhi; + (*mCorrIrregularDPhiRDz)(i, j) = -dPhi *radius0; (*mCorrIrregularDz)(i, j) = -dzDist; // distorted point (*mRIrregular)(i, j) = radius0 + drDist; - (*mPhiIrregular)(i, j) = phi0 + (dRPhi / radius0); + (*mPhiIrregular)(i, j) = phi0 + dPhi; (*mZIrregular)(i, j) = z0 + dzDist; /////////////// @@ -2060,11 +2093,11 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( // get global correction from j+1 drDist = (*mCorrDrDz)(i, j + 1); - dRPhi = (*mCorrDPhiRDz)(i, j + 1); dzDist = (*mCorrDz)(i, j + 1); radiusCorrection = radius0 + drDist; - phi = phi0 + dRPhi / radiusCorrection; + dPhi = (*mCorrDPhiRDz)(i, j + 1) /radius0; + phi = phi0 + dPhi; z = zList[j + 1] + dzDist; z1 = z - (zList[j+1] - zList[j]); @@ -2072,15 +2105,171 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( while (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); //lookupLocalCorr->GetValue(radiusCorrection, phi, z, ddR, ddRPhi, ddZ); - localIntErOverEz = (intDrDzF->Eval(radiusCorrection, phi, z1) - intDrDzF->Eval(radiusCorrection, phi, z)) / (-1 * ezField); - localIntEPhiOverEz = (intDPhiDzF->Eval(radiusCorrection, phi, z1) - intDPhiDzF->Eval(radiusCorrection, phi, z)) / (-1 * ezField); + localIntErOverEz = (intDrDzF->Eval(radiusCorrection, phi, z1) - intDrDzF->Eval(radiusCorrection, phi, z)) / (ezField + intDzDzF->Eval(radiusCorrection,phi,z)); + localIntEPhiOverEz = (intDPhiDzF->Eval(radiusCorrection, phi, z1) - intDPhiDzF->Eval(radiusCorrection, phi, z)) / (ezField + intDzDzF->Eval(radiusCorrection,phi,z)); localIntDeltaEz = intDzDzF->Eval(radiusCorrection, phi, z1) - intDzDzF->Eval(radiusCorrection, phi, z); ddR = fC0 * localIntErOverEz + fC1 * localIntEPhiOverEz; ddRPhi = fC0 * localIntEPhiOverEz - fC1 * localIntErOverEz; - ddZ = localIntDeltaEz *AliTPCPoissonSolver::fgkdvdE * AliTPCPoissonSolver::fgkdvdE; // two times? + ddZ = -1 * localIntDeltaEz * AliTPCPoissonSolver::fgkdvdE; // two times? + + + drDist += ddR; + dzDist += ddZ; + dPhi += ddRPhi/radiusCorrection; + + (*mCorrDrDz)(i, j) = drDist; + (*mCorrDPhiRDz)(i, j) = dPhi * radius0; + (*mCorrDz)(i, j) = dzDist; + + } + } + } +} + +/// See explanation at LocalDistCorrDz +/// +/// +/// \param matricesEr TMatrixD** electric field for \f$r\f$ component +/// \param matricesEPhi TMatrixD** electric field for \f$\phi\f$ component +/// \param matricesEz TMatrixD** electric field for \f$z\f$ component +/// \param matricesCorrDrDz TMatrixD** local correction \f$\hat{\delta}_{r}\f$ +/// \param matricesCorrDPhiRDz TMatrixD** local correction \f$r \hat{\delta}_{\phi}\f$ +/// \param matricesCorrDz TMatrixD** local correction \f$ \hat{\delta}_{z}\f$ +/// \param nRRow Int_t Number of nRRow in r-direction +/// \param nZColumn Int_t Number of nZColumn in z-direction +/// \param phiSlice Int_t Number of phi slices in \f$ phi \f$ direction +/// \param gridSizeZ const Float_t grid size in z direction +/// \param ezField const Double_t ezField calculate from the invoking operation +/// +/// \pre matricesEr, matricesEPhi, matrices Ez are provided +/// \post Local correction are computed according simplified Langevin equation +/// ~~~ +/// matricesCorrDz,matricesCorrDPhiRDz,matricesDistDz +/// ~~~ +void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDzWithLookUp ( AliTPCLookUpTable3DInterpolatorD *lookupLocalDist, TMatrixD** matricesGDistDrDz, TMatrixD** matricesGDistDPhiRDz, TMatrixD** matricesGDistDz, AliTPCLookUpTable3DInterpolatorD *lookupLocalCorr, TMatrixD** matricesGCorrDrDz, TMatrixD** matricesGCorrDPhiRDz, TMatrixD** matricesGCorrDz, const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, Double_t *rList, Double_t *phiList, Double_t *zList ) { + + Float_t drDist, dRPhi, dzDist, ddR, ddRPhi, ddZ; + Float_t radius0, phi0, z0, radius, phi, z, radiusCorrection; + radiusCorrection = 0.0; + radius = 0.0; + TMatrixD *mDistDrDz; + TMatrixD *mDistDPhiRDz; + TMatrixD *mDistDz; + TMatrixD *mCorrDrDz; + TMatrixD *mCorrDPhiRDz; + TMatrixD *mCorrDz; + Int_t j = nZColumn - 1; + + + // allocate look up for temporal + AliTPCLookUpTable3DInterpolatorD *lookupGlobalDistTemp = + new AliTPCLookUpTable3DInterpolatorD( + nRRow, matricesGDistDrDz, rList, phiSlice, matricesGDistDPhiRDz, phiList, nZColumn, matricesGDistDz, + zList, 2); + + + z0 = zList[j]; + + for (Int_t m = 0; m < phiSlice; m++) { + phi0 = phiList[m]; + + mDistDrDz = matricesGDistDrDz[m]; + mDistDPhiRDz = matricesGDistDPhiRDz[m]; + mDistDz = matricesGDistDz[m]; + + // + mCorrDrDz = matricesGCorrDrDz[m]; + mCorrDPhiRDz = matricesGCorrDPhiRDz[m]; + mCorrDz = matricesGCorrDz[m]; + + + for (Int_t i = 0; i < nRRow; i++) { + // do from j to 0 + // follow the drift + radius0 = rList[i]; + phi = phi0; + radius = radius0; + + drDist = 0.0; + dRPhi = 0.0; + dzDist = 0.0; + ddRPhi = 0.0; + + /// + (*mDistDrDz)(i, j) = drDist; + (*mDistDPhiRDz)(i, j) = dRPhi; + (*mDistDz)(i, j) = dzDist; + } + } + + // from j one column near end cap + for (j = nZColumn - 2; j >= 0; j--) { + + z0 = zList[j]; + for (Int_t m = 0; m < phiSlice; m++) { + phi0 = phiList[m]; + mDistDrDz = matricesGDistDrDz[m]; + mDistDPhiRDz = matricesGDistDPhiRDz[m]; + mDistDz = matricesGDistDz[m]; + + // + mCorrDrDz = matricesGCorrDrDz[m]; + mCorrDPhiRDz = matricesGCorrDPhiRDz[m]; + mCorrDz = matricesGCorrDz[m]; + + + + for (Int_t i = 0; i < nRRow; i++) { + + // do from j to 0 + // follow the drift + radius0 = rList[i]; + phi = phi0; + radius = radius0; + z = z0; + + lookupLocalDist->GetValue(radius,phi,z,ddR,ddRPhi,ddZ); + + phi += ddRPhi/radius; + radius = radius0 + ddR; + z = zList[j+1] + ddZ; + + if (j < nZColumn - 2) { + lookupGlobalDistTemp->GetValue(radius,phi,z,drDist,dRPhi,dzDist); + } + + (*mDistDrDz)(i,j) = drDist + ddR; + (*mDistDPhiRDz)(i,j) = dRPhi + ddRPhi; + (*mDistDz)(i,j) = dzDist + ddZ; + + if (j > 0) { + (*mDistDrDz)(i,j) = drDist + ddR; + (*mDistDPhiRDz)(i,j) = dRPhi + ddRPhi; + (*mDistDz)(i,j) = dzDist + ddZ; + } + // copy to 1D for being able to interpolate at next step + + + + // put the radius to the original value + if (j == nZColumn - 2) radiusCorrection = radius0; + + // get global correction from j+1 + drDist = (*mCorrDrDz)(i, j + 1); + dRPhi = (*mCorrDPhiRDz)(i, j + 1); + dzDist = (*mCorrDz)(i, j + 1); + + radiusCorrection = radius0 + drDist; + phi = phi0 + dRPhi / radiusCorrection; + z = zList[j + 1] + dzDist; + + while (phi < 0.0) phi = TMath::TwoPi() + phi; + while (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); + + lookupLocalCorr->GetValue(radiusCorrection, phi, z, ddR, ddRPhi, ddZ); drDist += ddR; dzDist += ddZ; @@ -2092,19 +2281,24 @@ void AliTPCSpaceCharge3DCalc::IntegrateDistCorrDriftLineDz( } } + + lookupGlobalDistTemp->CopyFromMatricesToInterpolator(j); + if (j > 0) lookupGlobalDistTemp->CopyFromMatricesToInterpolator(j-1); } + delete lookupGlobalDistTemp; } + /// -/// \param lookupGlobal -/// \param lookupRDz -/// \param lookupPhiRDz -/// \param lookupDz -/// \param nRRow -/// \param nZColumn -/// \param phiSlice -/// \param rList -/// \param phiList -/// \param zList +/// \param lookupGlobal +/// \param lookupRDz +/// \param lookupPhiRDz +/// \param lookupDz +/// \param nRRow +/// \param nZColumn +/// \param phiSlice +/// \param rList +/// \param phiList +/// \param zList void AliTPCSpaceCharge3DCalc::FillLookUpTable(AliTPCLookUpTable3DInterpolatorD *lookupGlobal, TMatrixD **lookupRDz, TMatrixD **lookupPhiRDz, TMatrixD **lookupDz, const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, const Double_t *rList, @@ -2212,9 +2406,6 @@ void AliTPCSpaceCharge3DCalc::GetCorrectionCylACIrregular(const Float_t x[], Sho r = x[0]; phi = x[1]; - if (phi < 0) phi += TMath::TwoPi(); // Table uses phi from 0 to 2*Pi - if (phi > TMath::TwoPi()) phi = phi - TMath::TwoPi(); // Table uses phi from 0 to 2*Pi - z = x[2]; // Create temporary copy of x[2] if ((roc % 36) < 18) { @@ -2227,24 +2418,19 @@ void AliTPCSpaceCharge3DCalc::GetCorrectionCylACIrregular(const Float_t x[], Sho if (sign == -1 && z > -AliTPCPoissonSolver::fgkZOffSet) z = -AliTPCPoissonSolver::fgkZOffSet; // Protect against discontinuity at CE - if ((sign == 1 && z < -1e-16) || (sign == -1 && z > -1e-16)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetCorrectionCylACIrregular","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + if ((sign == 1 && z < 0) || (sign == -1 && z > 0)) // just a consistency check + Error("AliTPCSpaceChargeCalc3D::GetCorrectionCylACIrregular","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); // get distortion from irregular table - Int_t iAnchor = TMath::FloorNint((r - AliTPCPoissonSolver::fgkIFCRadius) / gridSizeR); - Int_t kAnchor = TMath::FloorNint(phi / gridSizePhi); - Int_t zAnchor = TMath::FloorNint(z / gridSizeZ); - if (z > -1e-16) { - fLookupIntCorrIrregularA->GetValue(r, phi, z, dR, dRPhi, dZ, iAnchor, kAnchor, zAnchor, fNRRows / 8 + 1, - fNPhiSlices / 8 + 1, fNZColumns / 8 + 1, 0); + if (z > 0) { + fLookupIntCorrIrregularA->GetValue(r, phi, z, dR, dRPhi, dZ); } else { - fLookupIntCorrIrregularC->GetValue(r, phi, -z, dR, dRPhi, dZ, iAnchor, kAnchor, -zAnchor, fNRRows / 8 + 1, - fNPhiSlices / 8 + 1, fNZColumns / 8 + 1, 0); + fLookupIntCorrIrregularC->GetValue(r, phi, -z, dR, dRPhi, dZ); dZ = -1 * dZ; } @@ -2254,14 +2440,74 @@ void AliTPCSpaceCharge3DCalc::GetCorrectionCylACIrregular(const Float_t x[], Sho dZ; // z distortion - (scaled with drift velocity dependency on the Ez field and the overall scaling factor) } -/// Get correction regular grid by following electron + + + +/// Get Correction from irregular table /// /// \param x /// \param roc /// \param dx +void AliTPCSpaceCharge3DCalc::GetCorrectionCylACIrregular(const Float_t x[], Short_t roc, Float_t dx[], const Int_t side) { + if (!fInitLookUp) { + Info("AliTPCSpaceCharge3DCalc::GetCorrectionCylACIrregular","Lookup table was not initialized! Performing the initialization now ..."); + InitSpaceCharge3DPoissonIntegralDz(129, 129, 144, 100, 1e-8); + } + + Double_t dR, dRPhi, dZ; + Double_t r, phi, z; + Int_t sign; + const Float_t gridSizeR = (AliTPCPoissonSolver::fgkOFCRadius - AliTPCPoissonSolver::fgkIFCRadius) / (fNRRows - 1); + const Float_t gridSizeZ = AliTPCPoissonSolver::fgkTPCZ0 / (fNZColumns - 1); + const Float_t gridSizePhi = TMath::TwoPi() / fNPhiSlices; + + r = x[0]; + phi = x[1]; + 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 < AliTPCPoissonSolver::fgkZOffSet) z = AliTPCPoissonSolver::fgkZOffSet; // Protect against discontinuity at CE + if (sign == -1 && z > -AliTPCPoissonSolver::fgkZOffSet) z = -AliTPCPoissonSolver::fgkZOffSet; // Protect against discontinuity at CE + + + if ((sign == 1 && z < 0) || (sign == -1 && z > 0)) // just a consistency check + Error("AliTPCSpaceChargeCalc3D::GetCorrectionCylACIrregular","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + + + // get distortion from irregular table + + + + if (side == 0) + fLookupIntCorrIrregularA->GetValue(r, phi, z, dR, dRPhi, dZ); + else { + fLookupIntCorrIrregularC->GetValue(r, phi, -z, dR, dRPhi, dZ); + dZ = -1 * dZ; + } + + dx[0] = fCorrectionFactor * dR; + dx[1] = fCorrectionFactor * dRPhi; + dx[2] = fCorrectionFactor * + dZ; // z distortion - (scaled with drift velocity dependency on the Ez field and the overall scaling factor) + +} + + + + +/// Get correction regular grid by following electron +/// +/// \param x +/// \param roc +/// \param dx void AliTPCSpaceCharge3DCalc::GetCorrectionCylAC(const Float_t x[], Short_t roc, Float_t dx[]) { if (!fInitLookUp) { - Info("AliTPCSpaceCharge3DCalc::GetCorrectionCylAC","Lookup table was not initialized! Performing the initialization now ..."); + Info("AliTPCSpaceCharge3DCalc::GetDistortionCylAC","Lookup table was not initialized! Performing the initialization now ..."); InitSpaceCharge3DPoissonIntegralDz(129, 129, 144, 100, 1e-8); } @@ -2287,7 +2533,7 @@ void AliTPCSpaceCharge3DCalc::GetCorrectionCylAC(const Float_t x[], Short_t roc, if ((sign == 1 && z < -1e-16) || (sign == -1 && z > -1e-16)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetCorrectionCylAC","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceChargeCalc3D::GetCorrectionCylAC","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-16) fLookupIntCorrA->GetValue(r, phi, z, dR, dRPhi, dZ); @@ -2364,20 +2610,61 @@ void AliTPCSpaceCharge3DCalc::GetCorrection(const Float_t x[], Short_t roc, Floa pCyl[0] = TMath::Sqrt(x[0] * x[0] + x[1] * x[1]); pCyl[1] = TMath::ATan2(x[1], x[0]); - - // normalize phi - while (pCyl[1] > TMath::Pi()) pCyl[1] -= TMath::TwoPi(); - while (pCyl[1] < -TMath::Pi()) pCyl[1] += TMath::TwoPi(); - pCyl[2] = x[2]; // Create temporary copy of x[2] - if (fCorrectionType == kRegularInterpolator) + if (fCorrectionType == kRegularInterpolator) { + while (pCyl[1] > TMath::Pi()) pCyl[1] -= TMath::TwoPi(); + while (pCyl[1] < -TMath::Pi()) pCyl[1] += TMath::TwoPi(); + GetCorrectionCylAC(pCyl, roc, dCyl); - else + } else { GetCorrectionCylACIrregular(pCyl, roc, dCyl); + } + + // Calculate distorted position + if (pCyl[0] > 0.0) { + pCyl[0] = pCyl[0] + fCorrectionFactor * dCyl[0]; + pCyl[1] = pCyl[1] + fCorrectionFactor * dCyl[1] / pCyl[0]; + } + + dCyl[2] = fCorrectionFactor * dCyl[2]; + + // distortion in x,y and z + dx[0] = (pCyl[0] * TMath::Cos(pCyl[1]) - x[0]); + dx[1] = (pCyl[0] * TMath::Sin(pCyl[1]) - x[1]); + dx[2] = dCyl[2]; + +} +/// +/// \param x +/// \param roc +/// \param dx +void AliTPCSpaceCharge3DCalc::GetCorrection(const Float_t x[], Short_t roc, Float_t dx[],const Int_t side) { + if (!fInitLookUp) { + Info("AliTPCSpaceCharge3DCalc::GetCorrection","Lookup table was not initialized! Performing the initialization now ..."); + InitSpaceCharge3DPoissonIntegralDz(129, 129, 144, 100, 1e-8); + } + + Float_t pCyl[3]; // a point in cylindrical coordinate + Float_t dCyl[3]; // distortion + + pCyl[0] = TMath::Sqrt(x[0] * x[0] + x[1] * x[1]); + pCyl[1] = TMath::ATan2(x[1], x[0]); + pCyl[2] = x[2]; // Create temporary copy of x[2] + + + if (fCorrectionType == kRegularInterpolator) { + while (pCyl[1] > TMath::Pi()) pCyl[1] -= TMath::TwoPi(); + while (pCyl[1] < -TMath::Pi()) pCyl[1] += TMath::TwoPi(); + + GetCorrectionCylAC(pCyl, roc, dCyl); + } else { + GetCorrectionCylACIrregular(pCyl, roc, dCyl,side); + } + // Calculate distorted position if (pCyl[0] > 0.0) { pCyl[0] = pCyl[0] + fCorrectionFactor * dCyl[0]; @@ -2392,7 +2679,10 @@ void AliTPCSpaceCharge3DCalc::GetCorrection(const Float_t x[], Short_t roc, Floa dx[2] = dCyl[2]; } -// Use 3D space charge map as an optional input + + + +/// Use 3D space charge map as an optional input /// The layout of the input histogram is assumed to be: (phi,r,z) /// Density histogram is expected to bin in C/m^3 /// @@ -2403,7 +2693,51 @@ void AliTPCSpaceCharge3DCalc::GetCorrection(const Float_t x[], Short_t roc, Floa void AliTPCSpaceCharge3DCalc::SetInputSpaceCharge(TH3 *hisSpaceCharge3D, Double_t norm) { fHistogram3DSpaceCharge = hisSpaceCharge3D; fInitLookUp = kFALSE; + + + Info("AliTPCSpaceCharge3DCalc:SetInputSpaceCharge","Set Input Space Charge by 3D"); + Double_t rMin = hisSpaceCharge3D->GetYaxis()->GetBinCenter(0); + Double_t rMax = hisSpaceCharge3D->GetYaxis()->GetBinUpEdge(hisSpaceCharge3D->GetYaxis()->GetNbins()); + Double_t zMin = hisSpaceCharge3D->GetZaxis()->GetBinCenter(0); + Double_t zMax = hisSpaceCharge3D->GetZaxis()->GetBinCenter(hisSpaceCharge3D->GetZaxis()->GetNbins()); + + Double_t radius0, z0, phi0; + TMatrixD *charge; + for (Int_t iSide = 0;iSide < 2;iSide++) { + for (Int_t k = 0; k < fNPhiSlices; k++) { + if (iSide == 0) + charge = fMatrixChargeA[k]; + else + charge = fMatrixChargeC[k]; + + phi0 = fListPhi[k]; + + for (Int_t i = 0; i < fNRRows; i++) { + radius0 = fListR[i]; + + + for (Int_t j = 0; j < fNZColumns; j++) { + z0 = fListZ[j]; + + if (radius0 > rMin && radius0 < rMax && z0 > zMin && z0 < zMax) { + + (*charge)(i, j) = norm * InterpolatePhi(hisSpaceCharge3D, phi0, radius0, z0); + } + //} + } // end j + } // end i + } // end phi + } + + fInterpolatorChargeA->SetValue(fMatrixChargeA); + if (fInterpolationOrder > 2) + fInterpolatorChargeA->InitCubicSpline(); + fInterpolatorChargeC->SetValue(fMatrixChargeC); + if (fInterpolationOrder > 2) + fInterpolatorChargeC->InitCubicSpline(); + } + /// SetInputCharge /// /// \param hisSpaceCharge3D TH3* histogram for space charge @@ -2437,7 +2771,7 @@ void AliTPCSpaceCharge3DCalc::SetInputSpaceCharge(TH3 *hisSpaceCharge3D, Double_ z0 = fListZ[j]; if (radius0 > rMin && radius0 < rMax && z0 > zMin && z0 < zMax) { - (*charge)(i, j) = norm * InterpolatePhi(hisSpaceCharge3D, phi0, radius0, z0); + (*charge)(i, j) = norm * InterpolatePhi(hisSpaceCharge3D, phi0, radius0, z0); } } // end j } // end i @@ -2453,6 +2787,7 @@ void AliTPCSpaceCharge3DCalc::SetInputSpaceCharge(TH3 *hisSpaceCharge3D, Double_ fInitLookUp = kFALSE; } + /// InterpolationPhi is only used for reading from TH3F (since it is not cylindrical) /// /// \param r @@ -2820,7 +3155,7 @@ void AliTPCSpaceCharge3DCalc::InverseLocalDistortionToElectricField( localIntErOverEz = fC0 * (*distDrDz)(i, j) - fC1 * (*distDPhiRDz)(i, j); localIntErOverEz = localIntErOverEz / (fC0 * fC0 + fC1 * fC1); localIntEPhiOverEz = ((*distDrDz)(i, j) - (fC0 * localIntErOverEz)) / fC1; - localIntDeltaEz = (*distDz)(i, j) / (AliTPCPoissonSolver::fgkdvdE * AliTPCPoissonSolver::fgkdvdE); // two times? + localIntDeltaEz = -1 * (*distDz)(i, j) / AliTPCPoissonSolver::fgkdvdE; // two times? (*tDistDrDz)(i, j) = localIntErOverEz; (*tDistDPhiRDz)(i, j) = localIntEPhiOverEz; (*tDistDz)(i, j) = localIntDeltaEz; @@ -2844,11 +3179,9 @@ void AliTPCSpaceCharge3DCalc::InverseLocalDistortionToElectricField( (*mEPhi)(i, 0) = ((*distDPhiRDz)(i, 0) / gridSizeZ) * -1 * ezField; (*mEz)(i, 0) = ((*distDz)(i, 0) / gridSizeZ); (*mEr)(i, nZColumn - 1) = - ((-0.5 * (*distDrDz)(i, nZColumn - 3) + 1.5 * (*distDrDz)(i, nZColumn - 2)) / gridSizeZ) * -1 * - ezField; + ((-0.5 * (*distDrDz)(i, nZColumn - 3) + 1.5 * (*distDrDz)(i, nZColumn - 2)) / gridSizeZ) *ezField; (*mEPhi)(i, nZColumn - 1) = ((-0.5 * (*distDPhiRDz)(i, nZColumn - 3) + 1.5 * (*distDPhiRDz)(i, nZColumn - 2)) / gridSizeZ) * - -1 * ezField; (*mEz)(i, nZColumn - 1) = (-0.5 * (*distDz)(i, nZColumn - 3) + 1.5 * (*distDz)(i, nZColumn - 2)) / gridSizeZ; @@ -2856,9 +3189,9 @@ void AliTPCSpaceCharge3DCalc::InverseLocalDistortionToElectricField( for (Int_t i = 0; i < nRRow; i++) { for (Int_t j = 1; j < nZColumn - 1; j++) { - (*mEr)(i, j) = (((*distDrDz)(i, j) + (*distDrDz)(i, j - 1)) / (2 * gridSizeZ)) * -1 * + (*mEr)(i, j) = (((*distDrDz)(i, j) + (*distDrDz)(i, j - 1)) / (2 * gridSizeZ)) * ezField; // z direction - (*mEPhi)(i, j) = (((*distDPhiRDz)(i, j) + (*distDPhiRDz)(i, j - 1)) / (2 * gridSizeZ)) * -1 * + (*mEPhi)(i, j) = (((*distDPhiRDz)(i, j) + (*distDPhiRDz)(i, j - 1)) / (2 * gridSizeZ)) * ezField; // z direction (*mEz)(i, j) = ((*distDz)(i, j) + (*distDz)(i, j - 1)) / (2 * gridSizeZ); // z direction } @@ -3282,6 +3615,7 @@ void AliTPCSpaceCharge3DCalc::CalculateEField( const Float_t gridSizePhi = TMath::TwoPi() / phiSlice; TMatrixD *matricesEr[phiSlice], *matricesEz[phiSlice], *matricesEPhi[phiSlice]; + //Allocate memory for electric field r,z, phi direction for (Int_t k = 0; k < phiSlice; k++) { matricesEr[k] = new TMatrixD(nRRow, nZColumn); @@ -3297,7 +3631,8 @@ void AliTPCSpaceCharge3DCalc::CalculateEField( phiSlice, gridSizeR, gridSizePhi, gridSizeZ, symmetry, AliTPCPoissonSolver::fgkIFCRadius); w.Stop(); - Info("AliTPCSpaceCharge3DCalc::CalculateEField","%s",Form("Time for calculation E-field CPU = %f s\n", w.CpuTime())); + Info("AliTPCSpaceCharge3DCalc::InitSpaceCharge3DPoissonIntegralDz","%s",Form("Time for calculation E-field CPU = %f s\n", w.CpuTime())); + //Integrate E(r)/E(z) from point of origin to pad plane @@ -3386,7 +3721,7 @@ void AliTPCSpaceCharge3DCalc::GetCorrectionCylNoDrift(const Float_t x[], const S /// Calculates the correction due the Space Charge effect within the TPC drift volume if (!fInitLookUp) { - Info("AliTPCSpaceCharge3DCalc::GetCorrectionCylNoDrift","Lookup table was not initialized! Performing the initialization now ..."); + Info("AliTPCSpaceCharge3DCalc::","Lookup table was not initialized! Performing the initialization now ..."); // InitSpaceCharge3DDistortion(); return; } @@ -3413,7 +3748,7 @@ void AliTPCSpaceCharge3DCalc::GetCorrectionCylNoDrift(const Float_t x[], const S if ((sign == 1 && z < 0) || (sign == -1 && z > 0)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetCorrectionCylNoDrift","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (sign == -1 && z < 0.0) { printf("call C side\n"); @@ -3649,7 +3984,7 @@ Double_t AliTPCSpaceCharge3DCalc::GetChargeCylAC(const Float_t x[], Short_t roc) if ((sign == 1 && z < 0) || (sign == -1 && z > 0)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetChargeCylAC","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-16) return fInterpolatorChargeA->GetValue(r, phi, z); @@ -3682,7 +4017,7 @@ Double_t AliTPCSpaceCharge3DCalc::GetPotentialCylAC(const Float_t x[], Short_t r if ((sign == 1 && z < 0) || (sign == -1 && z > 0)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetPotentialCylAC","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-16) return fInterpolatorPotentialA->GetValue(r, phi, z); @@ -3716,7 +4051,7 @@ Double_t AliTPCSpaceCharge3DCalc::GetInverseChargeCylAC(const Float_t x[], Short if ((sign == 1 && z < 0) || (sign == -1 && z > 0)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetInverseChargeCylAC","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-6) return fInterpolatorInverseChargeA->GetValue(r, phi, z); @@ -3752,7 +4087,7 @@ void AliTPCSpaceCharge3DCalc::GetLocalDistortionCyl(const Float_t x[], Short_t r if ((sign == 1 && z < -1e-16) || (sign == -1 && z > -1e-16)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetLocalDistortionCyl","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-16) fLookupDistA->GetValue(r, phi, z, dR, dRPhi, dZ); @@ -3795,7 +4130,7 @@ void AliTPCSpaceCharge3DCalc::GetElectricFieldCyl(const Float_t x[], Short_t roc if ((sign == 1 && z < -1e-16) || (sign == -1 && z > -1e-16)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetElectricFieldCyl","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-16) fLookupElectricFieldA->GetValue(r, phi, z, eR, ePhi, eZ); @@ -3819,7 +4154,7 @@ void AliTPCSpaceCharge3DCalc::GetElectricFieldCyl(const Float_t x[], Short_t roc /// \param dx void AliTPCSpaceCharge3DCalc::GetInverseLocalDistortionCyl(const Float_t x[], Short_t roc, Float_t dx[]) { if (!fInitLookUp) { - Info("AliTPCSpaceCharge3DCalc::GetInverseLocalDistortionCyl","Lookup table was not initialized! Performing the initialization now ..."); + Info("AliTPCSpaceCharge3DCalc::","Lookup table was not initialized! Performing the initialization now ..."); InitSpaceCharge3DPoissonIntegralDz(129, 129, 144, 100, 1e-8); } @@ -3845,7 +4180,7 @@ void AliTPCSpaceCharge3DCalc::GetInverseLocalDistortionCyl(const Float_t x[], Sh if ((sign == 1 && z < -1e-16) || (sign == -1 && z > -1e-16)) // just a consistency check - Error("AliTPCSpaceCharge3DCalc::GetInverseLocalDistortionCyl","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); + Error("AliTPCSpaceCharge3DCalc::","ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!"); if (z > -1e-16) fLookupInverseDistA->GetValue(r, phi, z, dR, dRPhi, dZ); @@ -3910,3 +4245,37 @@ void AliTPCSpaceCharge3DCalc::SetPotentialBoundaryAndChargeFormula(TFormula *vTe fInterpolatorChargeC->SetValue(fMatrixChargeC); fInterpolatorChargeC->InitCubicSpline(); } + + +/// Set interpolation +void AliTPCSpaceCharge3DCalc::SetInterpolationOrder(Int_t order) { + fInterpolationOrder = order; + + + fInterpolatorChargeA->SetOrder(fInterpolationOrder); + fInterpolatorChargeC->SetOrder(fInterpolationOrder); + fInterpolatorPotentialA->SetOrder(fInterpolationOrder); + fInterpolatorPotentialC->SetOrder(fInterpolationOrder); + fInterpolatorInverseChargeA->SetOrder(fInterpolationOrder); + fInterpolatorInverseChargeC->SetOrder(fInterpolationOrder); + + fLookupDistA->SetOrder(fInterpolationOrder); + + fLookupDistC->SetOrder(fInterpolationOrder); + + fLookupInverseDistA->SetOrder(fInterpolationOrder); + + fLookupInverseDistC->SetOrder(fInterpolationOrder); + + fLookupElectricFieldA->SetOrder(fInterpolationOrder); + fLookupElectricFieldC->SetOrder(fInterpolationOrder); + fLookupIntDistA->SetOrder(fInterpolationOrder); + fLookupIntDistC->SetOrder(fInterpolationOrder); + fLookupIntCorrA->SetOrder(fInterpolationOrder); + fLookupIntCorrC->SetOrder(fInterpolationOrder); + fLookupIntENoDriftA->SetOrder(fInterpolationOrder); + fLookupIntENoDriftC->SetOrder(fInterpolationOrder); + fLookupIntCorrIrregularA->SetOrder(fInterpolationOrder); + + fLookupIntCorrIrregularC->SetOrder(fInterpolationOrder); +} diff --git a/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.h b/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.h index 327e7af4..89a34561 100644 --- a/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.h +++ b/TPCSpaceChargeBase/AliTPCSpaceCharge3DCalc.h @@ -42,7 +42,7 @@ class AliTPCSpaceCharge3DCalc { TMatrixD **matricesCorrDrDzA, TMatrixD **matricesCorrDPhiRDzA, TMatrixD **matricesCorrDzA, TMatrixD **matricesDistDrDzC, TMatrixD **matricesDistDPhiRDzC, TMatrixD **matricesDistDzC, TMatrixD **matricesCorrDrDzC, TMatrixD **matricesCorrDPhiRDzC, TMatrixD **matricesCorrDzC, - TFormula *intErDzTestFunction, TFormula *intEPhiRDzTestFunction, TFormula *intDzTestFunction); + TFormula *intErDzTestFunction, TFormula *intEPhiRDzTestFunction, TFormula *intDzTestFunction, TFormula *ezF); void InitSpaceCharge3DPoisson(Int_t nRRow, Int_t nZColumn, Int_t phiSlice, Int_t maxIteration, Double_t stopConvergence); @@ -53,9 +53,11 @@ class AliTPCSpaceCharge3DCalc { void GetCorrectionCyl(const Float_t x[], Short_t roc, Float_t dx[]); void GetCorrectionCylAC(const Float_t x[], Short_t roc, Float_t dx[]); void GetCorrectionCylACIrregular(const Float_t x[], Short_t roc, Float_t dx[]); + void GetCorrectionCylACIrregular(const Float_t x[], Short_t roc, Float_t dx[],const Int_t side); void GetDistortion(const Float_t x[], Short_t roc, Float_t dx[]); void GetCorrection(const Float_t x[], Short_t roc, Float_t dx[]); + void GetCorrection(const Float_t x[], Short_t roc, Float_t dx[],const Int_t side); Double_t GetChargeCylAC(const Float_t x[], Short_t roc); Double_t GetPotentialCylAC(const Float_t x[], Short_t roc); @@ -75,6 +77,10 @@ class AliTPCSpaceCharge3DCalc { kIrregularInterpolator = 1, ///< use irregular interpolator for correction look up table }; + enum IntegrationStrategy { + kNaive = 0, ///< use interpolation with regular interpolator for correction look up table + kOpt = 1, ///< use irregular interpolator for correction look up table + }; void SetInputSpaceCharge(TH3 *hisSpaceCharge3D, Double_t norm); void SetInputSpaceCharge(TH3 *hisSpaceCharge3D) { SetInputSpaceCharge(hisSpaceCharge3D, 1); } void SetInputSpaceCharge(TH3 *hisSpaceCharge3D, Double_t norm, Int_t side); @@ -109,7 +115,8 @@ class AliTPCSpaceCharge3DCalc { AliTPCPoissonSolver *GetPoissonSolver() { return fPoissonSolver; } - void SetInterpolationOrder(Int_t order) { fInterpolationOrder = order; } + void SetInterpolationOrder(Int_t order); + Int_t GetInterpolationOrder() { return fInterpolationOrder; } @@ -213,10 +220,22 @@ class AliTPCSpaceCharge3DCalc { Float_t GetSpaceChargeDensity(Float_t r, Float_t phi, Float_t z); Float_t GetPotential(Float_t r, Float_t phi, Float_t z); void GetElectricFieldCyl(const Float_t x[], Short_t roc, Double_t dx[]); + struct Profile { + Double_t poissonSolverTime; + Double_t electricFieldTime; + Double_t localDistortionTime; + Double_t globalDistortionTime; + Double_t interpolationInitTime; + Int_t iteration; + }; + Profile GetProfile() { return myProfile; } + void SetIntegrationStrategy(Int_t integrationStrategy) { + fIntegrationStrategy = integrationStrategy; + } private: static const Int_t kNMaxPhi = 360; - + Profile myProfile; Int_t fNRRows; ///< the maximum on row-slices so far ~ 2cm slicing Int_t fNPhiSlices; ///< the maximum of phi-slices so far = (8 per sector) Int_t fNZColumns; ///< the maximum on column-slices so ~ 2cm slicing @@ -237,7 +256,7 @@ class AliTPCSpaceCharge3DCalc { Int_t fInterpolationOrder; ///> Order of interpolation (1-> tri linear, 2->Lagrange interpolation order 2, 3> cubic spline) Int_t fIrregularGridSize; ///> Size of irregular grid cubes for interpolation (min 3) Int_t fRBFKernelType; ///> RBF kernel type - + Int_t fIntegrationStrategy; ///> Strategy for integration TMatrixD *fMatrixIntDistDrEzA[kNMaxPhi]; //[kNMaxPhi] Matrices for storing Global distortion \f$ R \f$ direction for Side A TMatrixD *fMatrixIntDistDPhiREzA[kNMaxPhi]; //[kNMaxPhi] Matrices for storing Global \f$ \phi R \f$ Distortion for Side A @@ -285,6 +304,9 @@ class AliTPCSpaceCharge3DCalc { TMatrixD *fMatrixChargeC[kNMaxPhi]; //[kNMaxPhi] Matrices for storing input charge densities side C TMatrixD *fMatrixChargeInverseA[kNMaxPhi]; //[kNMaxPhi] Matrices for storing charge densities from backward algorithm side A TMatrixD *fMatrixChargeInverseC[kNMaxPhi]; //[kNMaxPhi] Matrices for storing charge densities from backward algorithm side C + + TMatrixD *fMatrixPotentialA[kNMaxPhi]; //[kNMaxPhi] Matrices for storing potential side A + TMatrixD *fMatrixPotentialC[kNMaxPhi]; //[kNMaxPhi] Matrices for storing potential side C AliTPC3DCylindricalInterpolator *fInterpolatorChargeA; //-> interpolator for charge densities side A AliTPC3DCylindricalInterpolator *fInterpolatorChargeC; //-> interpolator for charge densities side C @@ -359,7 +381,7 @@ class AliTPCSpaceCharge3DCalc { const Double_t *rList, const Double_t *phiList, const Double_t *zList); void IntegrateDistCorrDriftLineDz( - TFormula *intErDzTestFunction, TFormula *intEPhiRDzTestFunction, TFormula *intDzTestFunction, + TFormula *intErDzTestFunction, TFormula *intEPhiRDzTestFunction, TFormula *intDzTestFunction, TFormula *ezF, const Double_t ezField, TMatrixD **matricesGDistDrDz, TMatrixD **matricesGDistDPhiRDz, TMatrixD **matricesGDistDz, TMatrixD **matricesGCorrDrDz, TMatrixD **matricesGCorrDPhiRDz, TMatrixD **matricesGCorrDz, @@ -369,6 +391,8 @@ class AliTPCSpaceCharge3DCalc { const Double_t *rList, const Double_t *phiList, const Double_t *zList); + void IntegrateDistCorrDriftLineDzWithLookUp ( AliTPCLookUpTable3DInterpolatorD *lookupLocalDist, TMatrixD** matricesGDistDrDz, TMatrixD** matricesGDistDPhiRDz, TMatrixD** matricesGDistDz, AliTPCLookUpTable3DInterpolatorD *lookupLocalCorr, TMatrixD** matricesGCorrDrDz, TMatrixD** matricesGCorrDPhiRDz, TMatrixD** matricesGCorrDz, const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, Double_t *rList, Double_t *phiList, Double_t *zList ); + void FillLookUpTable(AliTPCLookUpTable3DInterpolatorD *lookupGlobal, TMatrixD **lookupRDz, TMatrixD **lookupPhiRDz, TMatrixD **lookupDz, const Int_t nRRow, const Int_t nZColumn, const Int_t phiSlice, const Double_t *rList, const Double_t *phiList, const Double_t *zList);