+ if (value==6){
+ // only shift in angle
+ // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
+ valT = (*mat)(2,0) +(*mat)(2,1)*dydx1;
+ }
+ //
+ return valT;
+}
+
+
+void AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
+ //
+ // Update track parameters t1
+ //
+ TMatrixD vecXk(5,1); // X vector
+ TMatrixD covXk(5,5); // X covariance
+ TMatrixD matHk(1,5); // vector to mesurement
+ TMatrixD measR(1,1); // measurement error
+ //TMatrixD matQk(5,5); // prediction noise vector
+ TMatrixD vecZk(1,1); // measurement
+ //
+ TMatrixD vecYk(1,1); // Innovation or measurement residual
+ TMatrixD matHkT(5,1);
+ TMatrixD matSk(1,1); // Innovation (or residual) covariance
+ TMatrixD matKk(5,1); // Optimal Kalman gain
+ TMatrixD mat1(5,5); // update covariance matrix
+ TMatrixD covXk2(5,5); //
+ TMatrixD covOut(5,5);
+ //
+ Double_t *param1=(Double_t*) track1.GetParameter();
+ Double_t *covar1=(Double_t*) track1.GetCovariance();
+
+ //
+ // copy data to the matrix
+ for (Int_t ipar=0; ipar<5; ipar++){
+ vecXk(ipar,0) = param1[ipar];
+ for (Int_t jpar=0; jpar<5; jpar++){
+ covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
+ }
+ }
+ //
+ //
+ //
+ vecZk(0,0) = track2.GetParameter()[4]; // 1/pt measurement from track 2
+ measR(0,0) = track2.GetCovariance()[14]; // 1/pt measurement error
+ if (noField) {
+ measR(0,0)*=0.000000001;
+ vecZk(0,0)=0.;
+ }
+ //
+ matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;
+ matHk(0,3)= 0; matHk(0,4)= 1; // vector to measurement
+ //
+ //
+ //
+ vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
+ matHkT=matHk.T(); matHk.T();
+ matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
+ matSk.Invert();
+ matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
+ vecXk += matKk*vecYk; // updated vector
+ mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
+ covXk2 = (mat1-(matKk*matHk));
+ covOut = covXk2*covXk;
+ //
+ //
+ //
+ // copy from matrix to parameters
+ if (0) {
+ covOut.Print();
+ vecXk.Print();
+ covXk.Print();
+ track1.Print();
+ track2.Print();
+ }
+
+ for (Int_t ipar=0; ipar<5; ipar++){
+ param1[ipar]= vecXk(ipar,0) ;
+ for (Int_t jpar=0; jpar<5; jpar++){
+ covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
+ }
+ }
+
+}
+
+void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
+ //
+ // Global Align -combine the partial alignment of pair of sectors
+ // minPoints - minimal number of points - don't use sector alignment wit smaller number
+ // sysError - error added to the alignemnt error
+ //
+ AliTPCcalibAlign * align = this;
+ TMatrixD * arrayAlign[72];
+ TMatrixD * arrayAlignDiff[72];
+ //
+ for (Int_t i=0;i<72; i++) {
+ TMatrixD * mat = new TMatrixD(4,4);
+ mat->UnitMatrix();
+ arrayAlign[i]=mat;
+ arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
+ }
+
+ TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
+ for (Int_t iter=0; iter<niter;iter++){
+ printf("Iter=\t%d\n",iter);
+ for (Int_t is0=0;is0<72; is0++) {
+ //
+ //TMatrixD *mati0 = arrayAlign[is0];
+ TMatrixD matDiff(4,4);
+ Double_t sumw=0;
+ for (Int_t is1=0;is1<72; is1++) {
+ Bool_t invers=kFALSE;
+ Int_t npoints=0;
+ TMatrixD covar;
+ TVectorD errors;
+ const TMatrixD *mat = align->GetTransformation(is0,is1,0);
+ if (mat){
+ npoints = align->GetFitter6(is0,is1)->GetNpoints();
+ if (npoints>minPoints){
+ align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
+ align->GetFitter6(is0,is1)->GetErrors(errors);
+ }
+ }
+ else{
+ invers=kTRUE;
+ mat = align->GetTransformation(is1,is0,0);
+ if (mat) {
+ npoints = align->GetFitter6(is1,is0)->GetNpoints();
+ if (npoints>minPoints){
+ align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
+ align->GetFitter6(is1,is0)->GetErrors(errors);
+ }
+ }
+ }
+ if (!mat) continue;
+ if (npoints<minPoints) continue;
+ //
+ Double_t weight=1;
+ if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
+ if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
+ if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
+ if (is1%36!=is0%36) weight*=1/2.; //Not up-down
+ weight/=(errors[4]*errors[4]+sysError*sysError); // wieghting with error in Y
+ //
+ //
+ TMatrixD matT = *mat;
+ if (invers) matT.Invert();
+ TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
+ diffMat-=(*arrayAlign[is0]);
+ matDiff+=weight*diffMat;
+ sumw+=weight;
+
+ (*cstream)<<"LAlign"<<
+ "iter="<<iter<<
+ "s0="<<is0<<
+ "s1="<<is1<<
+ "npoints="<<npoints<<
+ "m60.="<<arrayAlign[is0]<<
+ "m61.="<<arrayAlign[is1]<<
+ "m01.="<<&matT<<
+ "diff.="<<&diffMat<<
+ "cov.="<<&covar<<
+ "err.="<<&errors<<
+ "\n";
+ }
+ if (sumw>0){
+ matDiff*=1/sumw;
+ matDiff(0,0)=0;
+ matDiff(1,1)=0;
+ matDiff(1,1)=0;
+ matDiff(1,1)=0;
+ (*arrayAlignDiff[is0]) = matDiff;
+ }
+ }
+ for (Int_t is0=0;is0<72; is0++) {
+ if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
+ if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
+ //
+ (*cstream)<<"GAlign"<<
+ "iter="<<iter<<
+ "s0="<<is0<<
+ "m6.="<<arrayAlign[is0]<<
+ "\n";
+ }
+ }
+
+ delete cstream;
+ for (Int_t isec=0;isec<72;isec++){
+ fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
+ delete arrayAlignDiff[isec];
+ }
+}
+
+
+ Int_t AliTPCcalibAlign::RefitLinear(const AliTPCseed * track, Int_t isec, Double_t *fitParam, Int_t refSector, TMatrixD &tparam, TMatrixD&tcovar, Double_t xRef, Bool_t both){
+ //
+ // Refit tracklet linearly using clusters at given sector isec
+ // Clusters are rotated to the reference frame of sector refSector
+ //
+ // fit parameters and errors retruning in the fitParam
+ //
+ // seed - acces to the original clusters
+ // isec - sector to be refited
+ // fitParam -
+ // 0 lx
+ // 1 ly
+ // 2 dy/dz
+ // 3 lz
+ // 4 dz/dx
+ // 5 sx
+ // 6 sy
+ // 7 sdydx
+ // 8 sz
+ // 9 sdzdx
+ // ref sector is the sector defining ref frame - rotation
+ // return value - number of used clusters
+
+ const Int_t kMinClusterF=15;
+ const Int_t kdrow1 =10; // rows to skip at the end
+ const Int_t kdrow0 =3; // rows to skip at beginning
+ const Float_t kedgeyIn=2.5;
+ const Float_t kedgeyOut=4.0;
+ const Float_t kMaxDist=5; // max distance -in sigma
+ const Float_t kMaxCorrY=0.05; // max correction
+ //
+ Double_t dalpha = 0;
+ if ((refSector%18)!=(isec%18)){
+ dalpha = -((refSector%18)-(isec%18))*TMath::TwoPi()/18.;
+ }
+ Double_t ca = TMath::Cos(dalpha);
+ Double_t sa = TMath::Sin(dalpha);
+ //
+ //
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ //
+ // full track fit parameters
+ //
+ static TLinearFitter fyf(2,"pol1"); // change to static - suggestion of calgrind - 30 % of time
+ static TLinearFitter fzf(2,"pol1"); // relative to time of given class
+ TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+ TMatrixD covY(4,4),covZ(4,4);
+ Double_t chi2FacY =1;
+ Double_t chi2FacZ =1;
+ Int_t nf=0;
+ //
+ //
+ //
+ Float_t erry=0.1; // initial cluster error estimate
+ Float_t errz=0.1; // initial cluster error estimate
+ for (Int_t iter=0; iter<2; iter++){
+ fyf.ClearPoints();
+ fzf.ClearPoints();
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ //
+ if (c->GetDetector()%36!=(isec%36)) continue;
+ if (!both && c->GetDetector()!=isec) continue;
+
+ if (c->GetRow()<kdrow0) continue;
+ //cluster position in reference frame
+ Double_t lxR = ca*c->GetX()-sa*c->GetY();
+ Double_t lyR = +sa*c->GetX()+ca*c->GetY();
+ Double_t lzR = c->GetZ();
+
+ Double_t dx = lxR -xRef; // distance to reference X
+ Double_t x[2]={dx, dx*dx};
+
+ Double_t yfitR = pyf[0]+pyf[1]*dx; // fit value Y in ref frame
+ Double_t zfitR = pzf[0]+pzf[1]*dx; // fit value Z in ref frame
+ //
+ Double_t yfit = -sa*lxR + ca*yfitR; // fit value Y in local frame
+ //
+ if (iter==0 &&c->GetType()<0) continue;
+ if (iter>0){
+ if (TMath::Abs(lyR-yfitR)>kMaxDist*erry) continue;
+ if (TMath::Abs(lzR-zfitR)>kMaxDist*errz) continue;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (isec<36 && dedge<kedgeyIn) continue;
+ if (isec>35 && dedge<kedgeyOut) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ Double_t corrclY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),c->GetY(), c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs((corrtrY+corrclY)*0.5)>kMaxCorrY) continue;
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ }
+ fyf.AddPoint(x,lyR,erry);
+ fzf.AddPoint(x,lzR,errz);
+ }
+ nf = fyf.GetNpoints();
+ if (nf<kMinClusterF) return 0; // not enough points - skip
+ fyf.Eval();
+ fyf.GetParameters(pyf);
+ fyf.GetErrors(peyf);
+ fzf.Eval();
+ fzf.GetParameters(pzf);
+ fzf.GetErrors(pezf);
+ chi2FacY = TMath::Sqrt(fyf.GetChisquare()/(fyf.GetNpoints()-2.));
+ chi2FacZ = TMath::Sqrt(fzf.GetChisquare()/(fzf.GetNpoints()-2.));
+ peyf[0]*=chi2FacY;
+ peyf[1]*=chi2FacY;
+ pezf[0]*=chi2FacZ;
+ pezf[1]*=chi2FacZ;
+ erry*=chi2FacY;
+ errz*=chi2FacZ;
+ fyf.GetCovarianceMatrix(covY);
+ fzf.GetCovarianceMatrix(covZ);
+ for (Int_t i0=0;i0<2;i0++)
+ for (Int_t i1=0;i1<2;i1++){
+ covY(i0,i1)*=chi2FacY*chi2FacY;
+ covZ(i0,i1)*=chi2FacZ*chi2FacZ;
+ }
+ }
+ fitParam[0] = xRef;
+ //
+ fitParam[1] = pyf[0];
+ fitParam[2] = pyf[1];
+ fitParam[3] = pzf[0];
+ fitParam[4] = pzf[1];
+ //
+ fitParam[5] = 0;
+ fitParam[6] = peyf[0];
+ fitParam[7] = peyf[1];
+ fitParam[8] = pezf[0];
+ fitParam[9] = pezf[1];
+ //
+ //
+ tparam(0,0) = pyf[0];
+ tparam(1,0) = pyf[1];
+ tparam(2,0) = pzf[0];
+ tparam(3,0) = pzf[1];
+ //
+ tcovar(0,0) = covY(0,0);
+ tcovar(1,1) = covY(1,1);
+ tcovar(1,0) = covY(1,0);
+ tcovar(0,1) = covY(0,1);
+ tcovar(2,2) = covZ(0,0);
+ tcovar(3,3) = covZ(1,1);
+ tcovar(3,2) = covZ(1,0);
+ tcovar(2,3) = covZ(0,1);
+ return nf;
+}
+
+void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){
+ //
+ // Update the cluster residula histograms for setup with field
+ // Kalman track fitting is used
+ // Only high momenta primary tracks used
+ //
+ // 1. Apply selection
+ // 2. Refit the track - in-out
+ // - update the cluster delta in upper part
+ // 3. Refit the track - out-in
+ // - update the cluster delta histo lower part
+ //
+ const Double_t kPtCut=1.0; // pt
+ const Double_t kSnpCut=0.2; // snp cut
+ const Double_t kNclCut=120; //
+ const Double_t kVertexCut=1;
+ const Double_t kMaxDist=0.5; // max distance between tracks and cluster
+ const Double_t kEdgeCut = 2.5;
+ if (!fCurrentTrack) return;
+ if (!fCurrentFriendTrack) return;
+ Float_t vertexXY=0,vertexZ=0;
+ fCurrentTrack->GetImpactParameters(vertexXY,vertexZ);
+ if (TMath::Abs(vertexXY)>kVertexCut) return;
+ if (TMath::Abs(vertexZ)>kVertexCut) return;
+ if (TMath::Abs(seed->Pt())<kPtCut) return;
+ if (seed->GetNumberOfClusters()<kNclCut) return;
+ if (TMath::Abs(seed->GetSnp())>kSnpCut) return;
+ if (!fClusterDelta[0]) MakeResidualHistos();
+
+ Int_t detector=-1;
+ //
+ //
+ AliExternalTrackParam trackIn = *(fCurrentTrack->GetInnerParam());
+ AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
+ static Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ //
+ Int_t ncl=0;
+ for (Int_t irow=0; irow<160; irow++){
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ if (detector!=cl->GetDetector()%36) return; // cluster from different sectors
+ // skip such tracks
+ ncl++;
+ }
+ if (ncl<kNclCut) return;
+
+ Int_t nclIn=0,nclOut=0;
+ Double_t xyz[3];
+ //
+ // Refit out - store residual maps
+ //
+ for (Int_t irow=0; irow<160; irow++){
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ Int_t sector = cl->GetDetector();
+ Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
+ if (TMath::Abs(dalpha)>0.01){
+ if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
+ }
+ Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
+ Double_t cov[3]={0.01,0.,0.01};
+ AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
+ Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackOut.GetY());
+ cov[0]+=1./(irow+1.); // bigger error at boundary
+ cov[0]+=1./(160.-irow); // bigger error at boundary
+ cov[2]+=1./(irow+1.); // bigger error at boundary
+ cov[2]+=1./(160.-irow); // bigger error at boundary
+ cov[0]+=0.5/dedge; // bigger error close to the boundary
+ cov[2]+=0.5/dedge; // bigger error close to the boundary
+ cov[0]*=cov[0];
+ cov[2]*=cov[2];
+ if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
+
+ if (TMath::Abs(dedge)<kEdgeCut) continue;
+
+ if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
+ nclOut++;
+ trackOut.Update(&r[1],cov);
+ }
+ if (nclOut<kNclCut/2) continue;
+ if (cl->GetDetector()%36!=detector) continue;
+ //
+ // fill residual histogram
+ //
+ Double_t resVector[5];
+ trackOut.GetXYZ(xyz);
+ resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
+ if (resVector[1]<0) resVector[1]+=18;
+ resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
+ resVector[3]= cl->GetZ()/resVector[2];
+ //
+ resVector[0]= cl->GetY()-trackOut.GetY();
+ fClusterDelta[0]->Fill(resVector);
+ resVector[0]= cl->GetZ()-trackOut.GetZ();
+ fClusterDelta[1]->Fill(resVector);
+ }
+ //
+ // Refit in - store residual maps
+ //
+ for (Int_t irow=159; irow>=0; irow--){
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ Int_t sector = cl->GetDetector();
+ Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
+ if (TMath::Abs(dalpha)>0.01){
+ if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
+ }
+ Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
+ Double_t cov[3]={0.01,0.,0.01};
+ AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
+ Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackIn.GetY());
+ cov[0]+=1./(irow+1.); // bigger error at boundary
+ cov[0]+=1./(160.-irow); // bigger error at boundary
+ cov[2]+=1./(irow+1.); // bigger error at boundary
+ cov[2]+=1./(160.-irow); // bigger error at boundary
+ cov[0]+=0.5/dedge; // bigger error close to the boundary +-
+ cov[2]+=0.5/dedge; // bigger error close to the boundary +-
+ cov[0]*=cov[0];
+ cov[2]*=cov[2];
+ if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
+ if (TMath::Abs(dedge)<kEdgeCut) continue;
+
+
+ if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
+ nclIn++;
+ trackIn.Update(&r[1],cov);
+ }
+ if (nclIn<kNclCut/2) continue;
+ if (cl->GetDetector()%36!=detector) continue;
+ //
+ // fill residual histogram
+ //
+ Double_t resVector[5];
+ trackIn.GetXYZ(xyz);
+ resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
+ if (resVector[1]<0) resVector[1]+=18;
+ resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
+ resVector[3]= cl->GetZ()/resVector[2];
+ //
+ resVector[0]= cl->GetY()-trackIn.GetY();
+ fClusterDelta[0]->Fill(resVector);
+ resVector[0]= cl->GetZ()-trackIn.GetZ();
+ fClusterDelta[1]->Fill(resVector);
+ }
+
+
+}
+
+
+void AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
+ //
+ // Update Kalman filter of Alignment - only setup without filed
+ // IROC - OROC quadrants
+ //
+ if (TMath::Abs(AliTracker::GetBz())>0.5) return;
+ if (!fClusterDelta[0]) MakeResidualHistos();
+ const Int_t kMinClusterF=40;
+ const Int_t kMinClusterFit=10;
+ const Int_t kMinClusterQ=10;
+ //
+ const Int_t kdrow1Fit =5; // rows to skip from fit at the end
+ const Int_t kdrow0Fit =10; // rows to skip from fit at beginning
+ const Float_t kedgey=3.0;
+ const Float_t kMaxDist=1;
+ const Float_t kMaxCorrY=0.05;
+ const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
+ isec = isec%36; // use the hardware numbering
+ //
+ //
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ //
+ // full track fit parameters
+ //
+ static TLinearFitter fyf(2,"pol1"); // make it static - too much time for comiling of formula
+ static TLinearFitter fzf(2,"pol1"); // calgrind recomendation
+ TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+ TVectorD pyfc(2),pzfc(2);
+ Int_t nf=0;
+ //
+ // make full fit as reference
+ //
+ for (Int_t iter=0; iter<2; iter++){
+ fyf.ClearPoints();
+ fzf.ClearPoints();
+ for (Int_t irow=kdrow0Fit;irow<159-kdrow1Fit;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if ((c->GetDetector()%36)!=isec) continue;
+ if (c->GetRow()<kdrow0Fit) continue;
+ Double_t dx = c->GetX()-fXmiddle;
+ Double_t x[2]={dx, dx*dx};
+ if (iter==0 &&c->GetType()<0) continue;
+ if (iter==1){
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t zfit = pzf[0]+pzf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ }
+ if (TMath::Abs(x[0])<10){
+ fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
+ }
+ fzf.AddPoint(x,c->GetZ(),0.1);
+ }
+ nf = fyf.GetNpoints();
+ if (fyf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
+ if (fzf.GetNpoints()<kMinClusterF) return; // not enough points - skip
+ fyf.Eval();
+ fyf.GetParameters(pyf);
+ fyf.GetErrors(peyf);
+ fzf.Eval();
+ fzf.GetParameters(pzf);
+ fzf.GetErrors(pezf);
+ }
+ //
+ //
+ //
+ TVectorD vecX(160); // x vector
+ TVectorD vecY(160); // residuals vector
+ TVectorD vecZ(160); // residuals vector
+ TVectorD vPosG(3); //vertex position
+ TVectorD vPosL(3); // vertex position in the TPC local system
+ TVectorF vImpact(2); //track impact parameter
+ // Double_t tofSignal= fCurrentTrack->GetTOFsignal(); // tof signal
+ TVectorF tpcPosG(3); // global position of track at the middle of fXmiddle
+ Double_t lphi = TMath::ATan2(pyf[0],fXmiddle); // expected local phi angle - if vertex at 0
+ Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi; // expected global phi if vertex at 0
+ Double_t ky = pyf[0]/fXmiddle;
+ Double_t kyE =0, kzE=0; // ky and kz expected
+ Double_t alpha =2.*TMath::Pi()*(isec%18+0.5)/18.;
+ Double_t scos=TMath::Cos(alpha);
+ Double_t ssin=TMath::Sin(alpha);
+ const AliESDVertex* vertex = fCurrentEvent->GetPrimaryVertexTracks();
+ vertex->GetXYZ(vPosG.GetMatrixArray());
+ fCurrentTrack->GetImpactParameters(vImpact[0],vImpact[1]); // track impact parameters
+ //
+ tpcPosG[0]= scos*fXmiddle-ssin*pyf[0];
+ tpcPosG[1]=+ssin*fXmiddle+scos*pyf[0];
+ tpcPosG[2]=pzf[0];
+ vPosL[0]= scos*vPosG[0]+ssin*vPosG[1];
+ vPosL[1]=-ssin*vPosG[0]+scos*vPosG[1];
+ vPosL[2]= vPosG[2];
+ kyE = (pyf[0]-vPosL[1])/(fXmiddle-vPosL[0]);
+ kzE = (pzf[0]-vPosL[2])/(fXmiddle-vPosL[0]);
+ //
+ // get constrained parameters
+ //
+ Double_t xvertex=vPosL[0]-fXmiddle;
+ fyf.AddPoint(&xvertex,vPosL[1], 0.00001);
+ fzf.AddPoint(&xvertex,vPosL[2], 2.);
+ fyf.Eval();
+ fyf.GetParameters(pyfc);
+ fzf.Eval();
+ fzf.GetParameters(pzfc);
+ //
+ //
+ // Make Fitters and params for 5 fitters
+ // 1-4 OROC quadrants
+ // 0 IROC
+ //
+ static TLinearFitter *fittersY[5]={0,0,0,0,0}; // calgrind recomendation - fater to clear points
+ static TLinearFitter *fittersZ[5]={0,0,0,0,0}; // than create the fitter
+ if (fittersY[0]==0){
+ for (Int_t i=0;i<5;i++) {
+ fittersY[i] = new TLinearFitter(2,"pol1");
+ fittersZ[i] = new TLinearFitter(2,"pol1");
+ }
+ }
+ //
+ Int_t npoints[5];
+ TVectorD paramsY[5];
+ TVectorD errorsY[5];
+ TMatrixD covY[5];
+ Double_t chi2CY[5];
+ TVectorD paramsZ[5];
+ TVectorD errorsZ[5];
+ TMatrixD covZ[5];
+ Double_t chi2CZ[5];
+ for (Int_t i=0;i<5;i++) {
+ npoints[i]=0;
+ paramsY[i].ResizeTo(2);
+ errorsY[i].ResizeTo(2);
+ covY[i].ResizeTo(2,2);
+ paramsZ[i].ResizeTo(2);
+ errorsZ[i].ResizeTo(2);
+ covZ[i].ResizeTo(2,2);
+ fittersY[i]->ClearPoints();
+ fittersZ[i]->ClearPoints();
+ }
+ //
+ // Update fitters
+ //
+ Int_t countRes=0;
+ for (Int_t irow=0;irow<159;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if ((c->GetDetector()%36)!=isec) continue;
+ Double_t dx = c->GetX()-fXmiddle;
+ Double_t x[2]={dx, dx*dx};
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t zfit = pzf[0]+pzf[1]*dx;
+ Double_t yfitC = pyfc[0]+pyfc[1]*dx;
+ Double_t zfitC = pzfc[0]+pzfc[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ //
+ vecX[countRes]=c->GetX();
+ vecY[countRes]=c->GetY()-yfit;
+ vecZ[countRes]=c->GetZ()-zfit;
+ countRes++;
+ //
+ // Fill THnSparse cluster residuals
+ // use only primary candidates with ITS signal
+ if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){
+ Double_t resVector[5];
+ resVector[1]= 9.*gphi/TMath::Pi();
+ resVector[2]= TMath::Sqrt(c->GetX()*c->GetX()+c->GetY()*c->GetY());
+ resVector[3]= c->GetZ()/resVector[2];
+ //
+ //
+ resVector[0]= c->GetY()-yfitC;
+ fClusterDelta[0]->Fill(resVector);
+ resVector[0]= c->GetZ()-zfitC;
+ fClusterDelta[1]->Fill(resVector);
+ }
+ if (c->GetRow()<kdrow0Fit) continue;
+ if (c->GetRow()>159-kdrow1Fit) continue;
+ //
+
+ if (c->GetDetector()>35){
+ if (c->GetX()<fXquadrant){
+ if (yfit<-kPRFWidth) fittersY[1]->AddPoint(x,c->GetY(),0.1);
+ if (yfit<-kPRFWidth) fittersZ[1]->AddPoint(x,c->GetZ(),0.1);
+ if (yfit>kPRFWidth) fittersY[2]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fittersZ[2]->AddPoint(x,c->GetZ(),0.1);
+ }
+ if (c->GetX()>fXquadrant){
+ if (yfit<-kPRFWidth) fittersY[3]->AddPoint(x,c->GetY(),0.1);
+ if (yfit<-kPRFWidth) fittersZ[3]->AddPoint(x,c->GetZ(),0.1);
+ if (yfit>kPRFWidth) fittersY[4]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fittersZ[4]->AddPoint(x,c->GetZ(),0.1);
+ }
+ }
+ if (c->GetDetector()<36){
+ fittersY[0]->AddPoint(x,c->GetY(),0.1);
+ fittersZ[0]->AddPoint(x,c->GetZ(),0.1);
+ }
+ }
+ //
+ //
+ //
+ for (Int_t i=0;i<5;i++) {
+ npoints[i] = fittersY[i]->GetNpoints();
+ if (npoints[i]>=kMinClusterQ){
+ // Y fit
+ fittersY[i]->Eval();
+ Double_t chi2FacY = TMath::Sqrt(fittersY[i]->GetChisquare()/(fittersY[i]->GetNpoints()-2));
+ chi2CY[i]=chi2FacY;
+ fittersY[i]->GetParameters(paramsY[i]);
+ fittersY[i]->GetErrors(errorsY[i]);
+ fittersY[i]->GetCovarianceMatrix(covY[i]);
+ // renormalize errors
+ errorsY[i][0]*=chi2FacY;
+ errorsY[i][1]*=chi2FacY;
+ covY[i](0,0)*=chi2FacY*chi2FacY;
+ covY[i](0,1)*=chi2FacY*chi2FacY;
+ covY[i](1,0)*=chi2FacY*chi2FacY;
+ covY[i](1,1)*=chi2FacY*chi2FacY;
+ // Z fit
+ fittersZ[i]->Eval();
+ Double_t chi2FacZ = TMath::Sqrt(fittersZ[i]->GetChisquare()/(fittersZ[i]->GetNpoints()-2));
+ chi2CZ[i]=chi2FacZ;
+ fittersZ[i]->GetParameters(paramsZ[i]);
+ fittersZ[i]->GetErrors(errorsZ[i]);
+ fittersZ[i]->GetCovarianceMatrix(covZ[i]);
+ // renormalize errors
+ errorsZ[i][0]*=chi2FacZ;
+ errorsZ[i][1]*=chi2FacZ;
+ covZ[i](0,0)*=chi2FacZ*chi2FacZ;
+ covZ[i](0,1)*=chi2FacZ*chi2FacZ;
+ covZ[i](1,0)*=chi2FacZ*chi2FacZ;
+ covZ[i](1,1)*=chi2FacZ*chi2FacZ;
+ }
+ }
+ //
+ // void UpdateSectorKalman
+ //
+ for (Int_t i0=0;i0<5;i0++){
+ for (Int_t i1=i0+1;i1<5;i1++){
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ TMatrixD v0(4,1),v1(4,1); // measurement
+ TMatrixD cov0(4,4),cov1(4,4); // covariance
+ //
+ v0(0,0)= paramsY[i0][0]; v1(0,0)= paramsY[i1][0];
+ v0(1,0)= paramsY[i0][1]; v1(1,0)= paramsY[i1][1];
+ v0(2,0)= paramsZ[i0][0]; v1(2,0)= paramsZ[i1][0];
+ v0(3,0)= paramsZ[i0][1]; v1(3,0)= paramsZ[i1][1];
+ //covariance i0
+ cov0(0,0) = covY[i0](0,0);
+ cov0(1,1) = covY[i0](1,1);
+ cov0(1,0) = covY[i0](1,0);
+ cov0(0,1) = covY[i0](0,1);
+ cov0(2,2) = covZ[i0](0,0);
+ cov0(3,3) = covZ[i0](1,1);
+ cov0(3,2) = covZ[i0](1,0);
+ cov0(2,3) = covZ[i0](0,1);
+ //covariance i1
+ cov1(0,0) = covY[i1](0,0);
+ cov1(1,1) = covY[i1](1,1);
+ cov1(1,0) = covY[i1](1,0);
+ cov1(0,1) = covY[i1](0,1);
+ cov1(2,2) = covZ[i1](0,0);
+ cov1(3,3) = covZ[i1](1,1);
+ cov1(3,2) = covZ[i1](1,0);
+ cov1(2,3) = covZ[i1](0,1);
+ //
+ // And now update
+ //
+ if (TMath::Abs(pyf[1])<0.8){ //angular cut
+ UpdateSectorKalman(isec, i0,i1, &v0,&cov0,&v1,&cov1);
+ }
+ }
+ }
+
+ //
+ // Dump debug information
+ //
+ if (fStreamLevel>0){
+ // get vertex position
+ //
+ TTreeSRedirector *cstream = GetDebugStreamer();
+
+
+ if (cstream){
+ for (Int_t i0=0;i0<5;i0++){
+ for (Int_t i1=i0;i1<5;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ (*cstream)<<"sectorAlign"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "triggerClass="<<&fTriggerClass<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "isec="<<isec<< // current sector
+ //
+ "xref="<<fXmiddle<< // reference X
+ "vPosG.="<<&vPosG<< // vertex position in global system
+ "vPosL.="<<&vPosL<< // vertex position in local system
+ "vImpact.="<<&vImpact<< // track impact parameter
+ //"tofSignal="<<tofSignal<< // tof signal
+ "tpcPosG.="<<&tpcPosG<< // global position of track at the middle of fXmiddle
+ "lphi="<<lphi<< // expected local phi angle - if vertex at 0
+ "gphi="<<gphi<< // expected global phi if vertex at 0
+ "ky="<<ky<<
+ "kyE="<<kyE<< // expect ky - assiming pirmary track
+ "kzE="<<kzE<< // expected kz - assuming primary tracks
+ "salpha="<<alpha<< // sector alpha
+ "scos="<<scos<< // tracking cosinus
+ "ssin="<<ssin<< // tracking sinus
+ //
+ // full fit
+ //
+ "nf="<<nf<< // total number of points
+ "pyf.="<<&pyf<< // full OROC fit y
+ "pzf.="<<&pzf<< // full OROC fit z
+ "vX.="<<&vecX<< // x cluster
+ "vY.="<<&vecY<< // y residual cluster
+ "vZ.="<<&vecZ<< // z residual cluster
+ // quadrant and IROC fit
+ "i0="<<i0<< // quadrant number
+ "i1="<<i1<<
+ "n0="<<npoints[i0]<< // number of points
+ "n1="<<npoints[i1]<<
+ //
+ "py0.="<<¶msY[i0]<< // parameters
+ "py1.="<<¶msY[i1]<<
+ "ey0.="<<&errorsY[i0]<< // errors
+ "ey1.="<<&errorsY[i1]<<
+ "chiy0="<<chi2CY[i0]<< // chi2s
+ "chiy1="<<chi2CY[i1]<<
+ //
+ "pz0.="<<¶msZ[i0]<< // parameters
+ "pz1.="<<¶msZ[i1]<<
+ "ez0.="<<&errorsZ[i0]<< // errors
+ "ez1.="<<&errorsZ[i1]<<
+ "chiz0="<<chi2CZ[i0]<< // chi2s
+ "chiz1="<<chi2CZ[i1]<<
+ "\n";
+ }
+ }
+ }
+ }
+}
+
+void AliTPCcalibAlign::UpdateSectorKalman(Int_t sector, Int_t quadrant0, Int_t quadrant1, TMatrixD *const p0, TMatrixD *const c0, TMatrixD *const p1, TMatrixD *const c1 ){
+ //
+ //
+ // tracks are refitted at sector middle
+ //
+ if (fArraySectorIntParam.At(0)==NULL) MakeSectorKalman();
+ //
+ //
+ static TMatrixD matHk(4,30); // vector to mesurement
+ static TMatrixD measR(4,4); // measurement error
+ // static TMatrixD matQk(2,2); // prediction noise vector
+ static TMatrixD vecZk(4,1); // measurement
+ //
+ static TMatrixD vecYk(4,1); // Innovation or measurement residual
+ static TMatrixD matHkT(30,4); // helper matrix Hk transpose
+ static TMatrixD matSk(4,4); // Innovation (or residual) covariance
+ static TMatrixD matKk(30,4); // Optimal Kalman gain
+ static TMatrixD mat1(30,30); // update covariance matrix
+ static TMatrixD covXk2(30,30); // helper matrix
+ //
+ TMatrixD *vOrig = (TMatrixD*)(fArraySectorIntParam.At(sector));
+ TMatrixD *cOrig = (TMatrixD*)(fArraySectorIntCovar.At(sector));
+ //
+ TMatrixD vecXk(*vOrig); // X vector
+ TMatrixD covXk(*cOrig); // X covariance
+ //
+ //Unit matrix
+ //
+ for (Int_t i=0;i<30;i++)
+ for (Int_t j=0;j<30;j++){
+ mat1(i,j)=0;
+ if (i==j) mat1(i,j)=1;
+ }
+ //
+ //
+ // matHk - vector to measurement
+ //
+ for (Int_t i=0;i<4;i++)
+ for (Int_t j=0;j<30;j++){
+ matHk(i,j)=0;
+ }
+ //
+ // Measurement
+ // 0 - y
+ // 1 - ky
+ // 2 - z
+ // 3 - kz