X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCCorrection.cxx;h=1e0df2523bb390d6e271a6bd74bfc2626aa4a611;hb=04b2edad5f2807838a2d8f737640a16e749fc992;hp=2ba9416d1bd9fc3e6f86e37c856f874bc41df1d6;hpb=be67055b711d664cc10b27e88fd06cc2c5702871;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCCorrection.cxx b/TPC/AliTPCCorrection.cxx index 2ba9416d1bd..1e0df2523bb 100644 --- a/TPC/AliTPCCorrection.cxx +++ b/TPC/AliTPCCorrection.cxx @@ -50,20 +50,40 @@ #include #include #include +#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 "AliTPCExB.h" +#include "AliTPCCorrection.h" +#include "AliTPCRecoParam.h" #include "AliExternalTrackParam.h" #include "AliTrackPointArray.h" #include "TDatabasePDG.h" #include "AliTrackerBase.h" #include "AliTPCROC.h" +#include "THnSparse.h" +#include "AliTPCLaserTrack.h" #include "AliTPCCorrection.h" +#include "AliLog.h" ClassImp(AliTPCCorrection) // FIXME: the following values should come from the database -const Double_t AliTPCCorrection::fgkTPC_Z0 =249.7; // nominal gating grid position +const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm) const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm) const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point @@ -116,7 +136,7 @@ const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = { AliTPCCorrection::AliTPCCorrection() - : TNamed("correction_unity","unity"),fJLow(0),fKLow(0) + : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1) { // // default constructor @@ -124,7 +144,7 @@ AliTPCCorrection::AliTPCCorrection() } AliTPCCorrection::AliTPCCorrection(const char *name,const char *title) - : TNamed(name,title),fJLow(0),fKLow(0) +: TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1) { // // default constructor, that set the name and title @@ -220,13 +240,15 @@ void AliTPCCorrection::Print(Option_t* /*option*/) const { printf("TPC spacepoint correction: \"%s\"\n",GetTitle()); } -void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t /*t1*/,Float_t /*t2*/) { +void AliTPCCorrection:: 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 // - // SetOmegaTauT1T2(omegaTau, t1, t2); + fT1=t1; + fT2=t2; + //SetOmegaTauT1T2(omegaTau, t1, t2); } TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) { @@ -318,7 +340,6 @@ TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) { h->SetBinContent(iz,ir,r1-r0); } } - printf("SDF\n"); return h; } @@ -387,12 +408,11 @@ TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title, // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3) void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z, - const Double_t er[kNZ][kNR], Double_t &er_value ) -{ + const Double_t er[kNZ][kNR], Double_t &erValue ) { // // Interpolate table - 2D interpolation // - Double_t save_er[10] ; + Double_t saveEr[10] ; Search( kNZ, fgkZList, z, fJLow ) ; Search( kNR, fgkRList, r, fKLow ) ; @@ -402,16 +422,15 @@ void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ; for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) { - save_er[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ; + saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ; } - er_value = Interpolate( &fgkZList[fJLow], save_er, order, z ) ; + erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ; } Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[], - const Int_t order, const Double_t x ) -{ + const Int_t order, const Double_t x ) { // // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation. // @@ -430,8 +449,7 @@ Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t } -void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) -{ +void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) { // // Search an ordered table by starting at the most recently used point // @@ -478,8 +496,218 @@ void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Dou } +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; -AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){ + 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 coef1(rows) ; + std::vector 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){ // // Fit the track parameters - without and with distortion // 1. Space points in the TPC are simulated along the trajectory @@ -498,6 +726,7 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara // 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); @@ -513,14 +742,17 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara AliTrackPointArray pointArray0(npoints0); AliTrackPointArray pointArray1(npoints0); Double_t xyz[3]; - AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp); + AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp); // // simulate the track Int_t npoints=0; Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame for (Double_t radius=kRTPC0; radiusGaus(0,0.005); + xyz[1]+=gRandom->Gaus(0,0.005); + xyz[2]+=gRandom->Gaus(0,0.005); AliTrackPoint pIn0; // space point AliTrackPoint pIn1; Int_t sector= (xyz[2]>0)? 0:18; @@ -544,6 +776,7 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara npoints++; if (npoints>=npoints0) break; } + if (npoints0) ? ipoint: npoints-1-jpoint; + Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint; // AliTrackPoint pIn0; AliTrackPoint pIn1; @@ -574,9 +806,9 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point // - AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp); - AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp); - track.GetXYZ(xyz); + AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp); + AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp); + track.GetXYZ(xyz); // distorted track also propagated to the same reference radius // Double_t pointPos[2]={0,0}; Double_t pointCov[3]={0,0,0}; @@ -587,15 +819,19 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara pointCov[2]=prot0.GetCov()[5];//sigmaz^2 track0->Update(pointPos,pointCov); // - pointPos[0]=prot1.GetY();//local y - pointPos[1]=prot1.GetZ();//local z + Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X + Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X + Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X + + pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus + pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus pointCov[0]=prot1.GetCov()[3];//simay^2 pointCov[1]=prot1.GetCov()[4];//sigmayz pointCov[2]=prot1.GetCov()[5];//sigmaz^2 track1->Update(pointPos,pointCov); } - AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp); + AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp); track1->Rotate(track0->GetAlpha()); track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz()); @@ -678,7 +914,7 @@ TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){ -void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, TObjArray * corrArray, Int_t step, Bool_t debug ){ +void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){ // // Make a fit tree: // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX) @@ -695,7 +931,10 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t // corrArray - array with partial corrections // step - skipe entries - if 1 all entries processed - it is slow // debug 0 if debug on also space points dumped - it is slow - const Int_t kMinEntries=50; + const Double_t kMaxSnp = 0.85; + const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); + // const Double_t kB2C=-0.299792458e-3; + const Int_t kMinEntries=50; Double_t phi,theta, snp, mean,rms, entries; tinput->SetBranchAddress("theta",&theta); tinput->SetBranchAddress("phi", &phi); @@ -707,35 +946,35 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t // Int_t nentries=tinput->GetEntries(); Int_t ncorr=corrArray->GetEntries(); - Double_t corrections[100]; // + Double_t corrections[100]={0}; // Double_t tPar[5]; Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; Double_t refX=0; Int_t dir=0; - if (dtype==0) {refX=85; dir=-1;} - if (dtype==1) {refX=245; dir=1;} - if (dtype==2) {refX=0; dir=-1;} + if (dtype==0) {refX=85.; dir=-1;} + if (dtype==1) {refX=275.; dir=1;} + if (dtype==2) {refX=85.; dir=-1;} + if (dtype==3) {refX=360.; dir=-1;} // for (Int_t ientry=0; ientryGetEntry(ientry); + if (TMath::Abs(snp)>kMaxSnp) continue; tPar[0]=0; tPar[1]=theta*refX; tPar[2]=snp; tPar[3]=theta; - tPar[4]=0.00001; // should be calculated - non equal to 0 - cout<SetBranchAddress("dY.",&vecdY); + tree->SetBranchAddress("dZ.",&vecdZ); + tree->SetBranchAddress("eY.",&veceY); + tree->SetBranchAddress("eZ.",&veceZ); + tree->SetBranchAddress("LTr.",<r); + Int_t entries= tree->GetEntries(); + TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root"); + Double_t bz=AliTrackerBase::GetBz(); + // + + for (Int_t ientry=0; ientryGetEntry(ientry); + if (!ltr->GetVecGX()){ + ltr->UpdatePoints(); + } + TVectorD * delta= (itype==0)? vecdY:vecdZ; + TVectorD * err= (itype==0)? veceY:veceZ; + + for (Int_t irow=0; irow<159; irow++){ + Int_t nentries = 1000; + if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0; + if (veceY->GetMatrixArray()[irow]GetMatrixArray()[irow]GetVecPhi())[irow]; + Double_t theta =ltr->GetTgl(); + Double_t mean=delta->GetMatrixArray()[irow]; + Double_t gx=0,gy=0,gz=0; + Double_t snp = (*ltr->GetVecP2())[irow]; + Double_t rms = 0.1+err->GetMatrixArray()[irow]; + gx = (*ltr->GetVecGX())[irow]; + gy = (*ltr->GetVecGY())[irow]; + gz = (*ltr->GetVecGZ())[irow]; + Int_t bundle= ltr->GetBundle(); + Double_t dRrec=0; + // + // get delta R used in reconstruction + AliTPCcalibDB* calib=AliTPCcalibDB::Instance(); + AliTPCCorrection * correction = calib->GetTPCComposedCorrection(); + const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam(); + Double_t xyz0[3]={gx,gy,gz}; + Double_t oldR=TMath::Sqrt(gx*gx+gy*gy); + // + // old ExB correction + // + if(recoParam&&recoParam->GetUseExBCorrection()) { + Double_t xyz1[3]={gx,gy,gz}; + calib->GetExB()->Correct(xyz0,xyz1); + Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]); + dRrec=oldR-newR; + } + if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) { + Float_t xyz1[3]={gx,gy,gz}; + Int_t sector=(gz>0)?0:18; + correction->CorrectPoint(xyz1, sector); + Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]); + dRrec=oldR-newR; + } + + + (*pcstream)<<"fit"<< + "bz="<GetAxis(1)->GetNbins(); + Int_t first1=his0->GetAxis(1)->GetFirst(); + Int_t last1 =his0->GetAxis(1)->GetLast(); + // + Double_t bz=AliTrackerBase::GetBz(); + Int_t idim[4]={0,1,2,3}; + for (Int_t ibin1=first1; ibin1GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1)); + Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1); + THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1 + Int_t nbins3 = his1->GetAxis(3)->GetNbins(); + Int_t first3 = his1->GetAxis(3)->GetFirst(); + Int_t last3 = his1->GetAxis(3)->GetLast(); + // + + for (Int_t ibin3=first3-1; ibin3GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3)); + Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3); + if (ibin3GetAxis(3)->SetRangeUser(-1,1); + x3=0; + } + THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3 + Int_t nbins2 = his3->GetAxis(2)->GetNbins(); + Int_t first2 = his3->GetAxis(2)->GetFirst(); + Int_t last2 = his3->GetAxis(2)->GetLast(); + // + for (Int_t ibin2=first2; ibin2GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2)); + Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2); + TH1 * hisDelta = his3->Projection(0); + // + Double_t entries = hisDelta->GetEntries(); + Double_t mean=0, rms=0; + if (entries>kMinEntries){ + mean = hisDelta->GetMean(); + rms = hisDelta->GetRMS(); + } + (*pcstream)<