+
+
+//_____________________________________________________________________
+Long64_t AliTPCcalibAlign::Merge(TCollection* const list) {
+ //
+ // merge function
+ //
+ if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
+ if (!list)
+ return 0;
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj = 0;
+ iter->Reset();
+ Int_t count=0;
+ //
+ TString str1(GetName());
+ while((obj = iter->Next()) != 0)
+ {
+ AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
+ if (entry == 0) continue;
+ if (str1.CompareTo(entry->GetName())!=0) continue;
+ Add(entry);
+ count++;
+ }
+ return count;
+}
+
+
+void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
+ //
+ // Add entry - used for merging of compoents
+ //
+ for (Int_t i=0; i<72;i++){
+ for (Int_t j=0; j<72;j++){
+ if (align->fPoints[GetIndex(i,j)]<1) continue;
+ fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
+ //
+ //
+ //
+ for (Int_t itype=0; itype<10; itype++){
+ TH1 * his0=0, *his1=0;
+ his0 = GetHisto((HistoType)itype,i,j);
+ his1 = align->GetHisto((HistoType)itype,i,j);
+ if (his1){
+ if (his0) his0->Add(his1);
+ else {
+ his0 = GetHisto((HistoType)itype,i,j,kTRUE);
+ his0->Add(his1);
+ }
+ }
+ }
+ }
+ }
+ TLinearFitter *f0=0;
+ TLinearFitter *f1=0;
+ for (Int_t i=0; i<72;i++){
+ for (Int_t j=0; j<72;j++){
+ if (align->fPoints[GetIndex(i,j)]<1) continue;
+ //
+ //
+ // fitter12
+ f0 = GetFitter12(i,j);
+ f1 = align->GetFitter12(i,j);
+ if (f1){
+ if (f0) f0->Add(f1);
+ else {
+ fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
+ }
+ }
+ //
+ // fitter9
+ f0 = GetFitter9(i,j);
+ f1 = align->GetFitter9(i,j);
+ if (f1){
+ if (f0) f0->Add(f1);
+ else {
+ fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
+ }
+ }
+ f0 = GetFitter6(i,j);
+ f1 = align->GetFitter6(i,j);
+ if (f1){
+ if (f0) f0->Add(f1);
+ else {
+ fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
+ }
+ }
+ }
+ }
+ //
+ // Add Kalman filter
+ //
+ for (Int_t i=0;i<36;i++){
+ TMatrixD *par0 = (TMatrixD*)fArraySectorIntParam.At(i);
+ if (!par0){
+ MakeSectorKalman();
+ par0 = (TMatrixD*)fArraySectorIntParam.At(i);
+ }
+ TMatrixD *par1 = (TMatrixD*)align->fArraySectorIntParam.At(i);
+ if (!par1) continue;
+ //
+ TMatrixD *cov0 = (TMatrixD*)fArraySectorIntCovar.At(i);
+ TMatrixD *cov1 = (TMatrixD*)align->fArraySectorIntCovar.At(i);
+ UpdateSectorKalman(*par0,*cov0,*par1,*cov1);
+ }
+ if (!fSectorParamA){
+ MakeKalman();
+ }
+ if (align->fSectorParamA){
+ UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA);
+ UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC);
+ }
+ if (!fClusterDelta[1]) MakeResidualHistos();
+
+ for (Int_t i=0; i<6; i++){
+ if (i==0 || i==3){
+ delete fClusterDelta[i]; // memory problem do not fit into memory
+ fClusterDelta[i]=0; //
+ delete align->fClusterDelta[i]; // memory problem do not fit into memory
+ align->fClusterDelta[i]=0; //
+ }
+ if (i==3) continue; // skip constrained histo z
+ if (i==0) continue; // skip non constrained histo y
+ if (align->fClusterDelta[i]) fClusterDelta[i]->Add(align->fClusterDelta[i]);
+ }
+}
+
+Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){
+ //
+ // GetTransformed value
+ //
+ //
+ // x2 = a00*x1 + a01*y1 + a02*z1 + a03
+ // y2 = a10*x1 + a11*y1 + a12*z1 + a13
+ // z2 = a20*x1 + a21*y1 + a22*z1 + a23
+ // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
+ // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
+
+
+ const TMatrixD * mat = GetTransformation(s1,s2,type);
+ if (!mat) {
+ if (value==0) return x1;
+ if (value==1) return y1;
+ if (value==2) return z1;
+ if (value==3) return dydx1;
+ if (value==4) return dzdx1;
+ //
+ if (value==5) return dydx1;
+ if (value==6) return dzdx1;
+ return 0;
+ }
+ Double_t valT=0;
+
+ if (value==0){
+ valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
+ }
+
+ if (value==1){
+ valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
+ }
+ if (value==2){
+ valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
+ }
+ if (value==3){
+ // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
+ valT = (*mat)(1,0) +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
+ valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
+ }
+
+ if (value==4){
+ // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
+ valT = (*mat)(2,0) +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
+ valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
+ }
+ //
+ if (value==5){
+ // onlys shift in angle
+ // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
+ valT = (*mat)(1,0) +(*mat)(1,1)*dydx1;
+ }
+
+ 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::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
+ //
+ // Update Kalman filter of Alignment
+ // IROC - OROC quadrants
+ //
+ if (!fClusterDelta[0]) MakeResidualHistos();
+ const Int_t kMinClusterF=40;
+ 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;
+ }
+ fyf.AddPoint(x,c->GetY(),0.1);
+ fzf.AddPoint(x,c->GetZ(),0.1);
+ }
+ nf = fyf.GetNpoints();
+ if (nf<kMinClusterF) return; // not enough points - skip
+ fyf.Eval();
+ fyf.GetParameters(pyf);
+ fyf.GetErrors(peyf);
+ fzf.Eval();
+ fzf.GetParameters(pzf);
+ fzf.GetErrors(pezf);
+ }
+ //
+ //
+ //
+ TVectorD vecX(2*nf+kdrow0Fit+kdrow1Fit+5); // x vector
+ TVectorD vecY(2*nf+kdrow0Fit+kdrow1Fit+5); // residuals vector
+ TVectorD vecZ(2*nf+kdrow0Fit+kdrow1Fit+5); // 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.1+TMath::Abs(vImpact[0]));
+ fzf.AddPoint(&xvertex,vPosL[2], 0.1+TMath::Abs(vImpact[1]));
+ 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 (nf>100&&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]= c->GetX();
+ resVector[3]= c->GetY()/c->GetX();
+ resVector[4]= c->GetZ()/c->GetX();
+ //
+ resVector[0]= c->GetY()-yfit;
+ //fClusterDelta[0]->Fill(resVector);
+ resVector[0]= c->GetZ()-zfit;
+ fClusterDelta[1]->Fill(resVector);
+ //
+ resVector[0]= c->GetY()-yfitC;
+ fClusterDelta[2]->Fill(resVector);
+ resVector[0]= c->GetZ()-zfitC;
+ //fClusterDelta[3]->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
+
+ matHk(0,6*quadrant1+4) = 1.; // delta y
+ matHk(1,6*quadrant1+0) = 1.; // delta ky
+ matHk(2,6*quadrant1+5) = 1.; // delta z
+ matHk(3,6*quadrant1+1) = 1.; // delta kz
+ // bug fix 24.02 - aware of sign in dx
+ matHk(0,6*quadrant1+3) = -(*p0)(1,0); // delta x to delta y - through ky
+ matHk(2,6*quadrant1+3) = -(*p0)(3,0); // delta x to delta z - thorugh kz
+ matHk(2,6*quadrant1+2) = ((*p0)(0,0)); // y to delta z - through psiz
+ //
+ matHk(0,6*quadrant0+4) = -1.; // delta y
+ matHk(1,6*quadrant0+0) = -1.; // delta ky
+ matHk(2,6*quadrant0+5) = -1.; // delta z
+ matHk(3,6*quadrant0+1) = -1.; // delta kz
+ // bug fix 24.02 be aware of sign in dx
+ matHk(0,6*quadrant0+3) = ((*p0)(1,0)); // delta x to delta y - through ky
+ matHk(2,6*quadrant0+3) = ((*p0)(3,0)); // delta x to delta z - thorugh kz
+ matHk(2,6*quadrant0+2) = -((*p0)(0,0)); // y to delta z - through psiz
+
+ //
+ //
+
+ vecZk =(*p1)-(*p0); // measurement
+ measR =(*c1)+(*c0); // error of 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
+ covXk2= (mat1-(matKk*matHk));
+ covXk = covXk2*covXk;
+ //
+ //
+ (*cOrig)=covXk;
+ (*vOrig)=vecXk;
+}
+
+void AliTPCcalibAlign::MakeSectorKalman(){
+ //
+ // Make a initial Kalman paramaters for IROC - Quadrants alignment
+ //
+ TMatrixD param(5*6,1);
+ TMatrixD covar(5*6,5*6);
+ //
+ // Set inital parameters
+ //
+ for (Int_t ip=0;ip<5*6;ip++) param(ip,0)=0; // mean alignment to 0
+ //
+ for (Int_t iq=0;iq<5;iq++){
+ // Initial uncertinty
+ covar(iq*6+0,iq*6+0) = 0.002*0.002; // 2 mrad
+ covar(iq*6+1,iq*6+1) = 0.002*0.002; // 2 mrad rotation
+ covar(iq*6+2,iq*6+2) = 0.002*0.002; // 2 mrad
+ //
+ covar(iq*6+3,iq*6+3) = 0.02*0.02; // 0.2 mm
+ covar(iq*6+4,iq*6+4) = 0.02*0.02; // 0.2 mm translation
+ covar(iq*6+5,iq*6+5) = 0.02*0.02; // 0.2 mm
+ }
+
+ for (Int_t isec=0;isec<36;isec++){
+ fArraySectorIntParam.AddAt(param.Clone(),isec);
+ fArraySectorIntCovar.AddAt(covar.Clone(),isec);
+ }
+}
+
+void AliTPCcalibAlign::UpdateSectorKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
+ //
+ // Update kalman vector para0 with vector par1
+ // Used for merging
+ //
+ static TMatrixD matHk(30,30); // vector to mesurement
+ static TMatrixD measR(30,30); // measurement error
+ static TMatrixD vecZk(30,1); // measurement
+ //
+ static TMatrixD vecYk(30,1); // Innovation or measurement residual
+ static TMatrixD matHkT(30,30); // helper matrix Hk transpose
+ static TMatrixD matSk(30,30); // Innovation (or residual) covariance
+ static TMatrixD matKk(30,30); // Optimal Kalman gain
+ static TMatrixD mat1(30,30); // update covariance matrix
+ static TMatrixD covXk2(30,30); // helper matrix
+ //
+ TMatrixD vecXk(par0); // X vector
+ TMatrixD covXk(cov0); // 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 = mat1; // unit matrix
+ //
+ vecZk = par1; // measurement
+ measR = cov1; // error of 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
+ //matKk.Print();
+ vecXk += matKk*vecYk; // updated vector
+ covXk2= (mat1-(matKk*matHk));
+ covXk = covXk2*covXk;
+ CheckCovariance(covXk);
+ CheckCovariance(cov1);
+ //
+ par0 = vecXk; // update measurement param
+ cov0 = covXk; // update measurement covar
+}
+
+Double_t AliTPCcalibAlign::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
+ //
+ // Get position correction for given sector
+ //
+
+ TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
+ if (!param) return 0;
+ Int_t quadrant=0;
+ if(lx>fXIO){
+ if (lx<fXquadrant) {
+ if (ly<0) quadrant=1;
+ if (ly>0) quadrant=2;
+ }
+ if (lx>fXquadrant) {
+ if (ly<0) quadrant=3;
+ if (ly>0) quadrant=4;
+ }
+ }
+ Double_t a10 = (*param)(quadrant*6+0,0);
+ Double_t a20 = (*param)(quadrant*6+1,0);
+ Double_t a21 = (*param)(quadrant*6+2,0);
+ Double_t dx = (*param)(quadrant*6+3,0);
+ Double_t dy = (*param)(quadrant*6+4,0);
+ Double_t dz = (*param)(quadrant*6+5,0);
+ Double_t deltaX = lx-fXIO;
+ if (coord==0) return dx;
+ if (coord==1) return dy+deltaX*a10;
+ if (coord==2) return dz+deltaX*a20+ly*a21;
+ return 0;
+}
+
+Double_t AliTPCcalibAlign::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
+ //
+ //
+ //
+ if (!Instance()) return 0;
+ return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz);
+}
+
+void AliTPCcalibAlign::MakeKalman(){
+ //
+ // Make a initial Kalman paramaters for sector Alignemnt
+ //
+ fSectorParamA = new TMatrixD(6*36+6,1);
+ fSectorParamC = new TMatrixD(6*36+6,1);
+ fSectorCovarA = new TMatrixD(6*36+6,6*36+6);
+ fSectorCovarC = new TMatrixD(6*36+6,6*36+6);
+ //
+ // set starting parameters at 0
+ //
+ for (Int_t isec=0;isec<37;isec++)
+ for (Int_t ipar=0;ipar<6;ipar++){
+ (*fSectorParamA)(isec*6+ipar,0) =0;
+ (*fSectorParamC)(isec*6+ipar,0) =0;
+ }
+ //
+ // set starting covariance
+ //
+ for (Int_t isec=0;isec<36;isec++)
+ for (Int_t ipar=0;ipar<6;ipar++){
+ if (ipar<3){
+ (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.002*0.002; // 2 mrad
+ (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.002*0.002;
+ }
+ if (ipar>=3){
+ (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.02*0.02; // 0.2 mm
+ (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.02*0.02;
+ }
+ }
+ (*fSectorCovarA)(36*6+0,36*6+0) =0.04; // common shift y up-up
+ (*fSectorCovarA)(36*6+1,36*6+1) =0.04; // common shift y down-down
+ (*fSectorCovarA)(36*6+2,36*6+2) =0.04; // common shift y up-down
+ (*fSectorCovarA)(36*6+3,36*6+3) =0.004; // common shift phi up-up
+ (*fSectorCovarA)(36*6+4,36*6+4) =0.004; // common shift phi down-down
+ (*fSectorCovarA)(36*6+5,36*6+5) =0.004; // common shift phi up-down
+ //
+ (*fSectorCovarC)(36*6+0,36*6+0) =0.04; // common shift y up-up
+ (*fSectorCovarC)(36*6+1,36*6+1) =0.04; // common shift y down-down
+ (*fSectorCovarC)(36*6+2,36*6+2) =0.04; // common shift y up-down
+ (*fSectorCovarC)(36*6+3,36*6+3) =0.004; // common shift phi up-up
+ (*fSectorCovarC)(36*6+4,36*6+4) =0.004; // common shift phi down-down
+ (*fSectorCovarC)(36*6+5,36*6+5) =0.004; // common shift phi up-down
+}
+
+void AliTPCcalibAlign::UpdateKalman(Int_t sector0, Int_t sector1, TMatrixD &p0, TMatrixD &c0, TMatrixD &p1, TMatrixD &c1){
+ //
+ // Update Kalman parameters
+ // Note numbering from 0..36 0..17 IROC 18..35 OROC
+ //
+ //
+ if (fSectorParamA==NULL) MakeKalman();
+ if (CheckCovariance(c0)>0) return;
+ if (CheckCovariance(c1)>0) return;
+ const Int_t nelem = 6*36+6;
+ //
+ //
+ static TMatrixD matHk(4,nelem); // vector to mesurement
+ static TMatrixD measR(4,4); // measurement error
+ static TMatrixD vecZk(4,1); // measurement
+ //
+ static TMatrixD vecYk(4,1); // Innovation or measurement residual
+ static TMatrixD matHkT(nelem,4); // helper matrix Hk transpose
+ static TMatrixD matSk(4,4); // Innovation (or residual) covariance
+ static TMatrixD matKk(nelem,4); // Optimal Kalman gain
+ static TMatrixD mat1(nelem,nelem); // update covariance matrix
+ static TMatrixD covXk2(nelem,nelem); // helper matrix
+ //
+ TMatrixD *vOrig = 0;
+ TMatrixD *cOrig = 0;
+ vOrig = (sector0%36>=18) ? fSectorParamA:fSectorParamC;
+ cOrig = (sector0%36>=18) ? fSectorCovarA:fSectorCovarC;
+ //
+ Int_t sec0= sector0%18;
+ Int_t sec1= sector1%18;
+ if (sector0>35) sec0+=18;
+ if (sector1>35) sec1+=18;
+ //
+ TMatrixD vecXk(*vOrig); // X vector
+ TMatrixD covXk(*cOrig); // X covariance
+ //
+ //Unit matrix
+ //
+ for (Int_t i=0;i<nelem;i++)
+ for (Int_t j=0;j<nelem;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<nelem;j++){
+ matHk(i,j)=0;
+ }
+ //
+ // Measurement
+ // 0 - y
+ // 1 - ky
+ // 2 - z
+ // 3 - kz
+
+ matHk(0,6*sec1+4) = 1.; // delta y
+ matHk(1,6*sec1+0) = 1.; // delta ky
+ matHk(2,6*sec1+5) = 1.; // delta z
+ matHk(3,6*sec1+1) = 1.; // delta kz
+ matHk(0,6*sec1+3) = p0(1,0); // delta x to delta y - through ky
+ matHk(2,6*sec1+3) = p0(3,0); // delta x to delta z - thorugh kz
+ matHk(2,6*sec1+2) = p0(0,0); // y to delta z - through psiz
+ //
+ matHk(0,6*sec0+4) = -1.; // delta y
+ matHk(1,6*sec0+0) = -1.; // delta ky
+ matHk(2,6*sec0+5) = -1.; // delta z
+ matHk(3,6*sec0+1) = -1.; // delta kz
+ matHk(0,6*sec0+3) = -p0(1,0); // delta x to delta y - through ky
+ matHk(2,6*sec0+3) = -p0(3,0); // delta x to delta z - thorugh kz
+ matHk(2,6*sec0+2) = -p0(0,0); // y to delta z - through psiz
+
+ Int_t dsec = (sector1%18)-(sector0%18);
+ if (dsec<-2) dsec+=18;
+ if (TMath::Abs(dsec)==1){
+ //
+ // Left right systematic fit part
+ //
+ Double_t dir = 0;
+ if (dsec>0) dir= 1.;
+ if (dsec<0) dir=-1.;
+ if (sector0>35&§or1>35){
+ matHk(0,36*6+0)=dir;
+ matHk(1,36*6+3+0)=dir;
+ }
+ if (sector0<36&§or1<36){
+ matHk(0,36*6+1)=dir;
+ matHk(1,36*6+3+1)=dir;
+ }
+ if (sector0<36&§or1>35){
+ matHk(0,36*6+2)=dir;
+ matHk(1,36*6+3+2)=dir;
+ }
+ if (sector0>35&§or1<36){
+ matHk(0,36*6+2)=-dir;
+ matHk(1,36*6+3+2)=-dir;
+ }
+ }
+ //
+ //
+ vecZk =(p1)-(p0); // measurement
+ measR =(c1)+(c0); // error of 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
+ covXk2= (mat1-(matKk*matHk));
+ covXk = covXk2*covXk;
+
+ if (CheckCovariance(covXk)>0) return;
+
+ //
+ //
+ (*cOrig)=covXk;
+ (*vOrig)=vecXk;
+}
+
+
+void AliTPCcalibAlign::UpdateKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
+ //
+ // Update kalman vector para0 with vector par1
+ // Used for merging
+ //
+ Int_t nelem = 6*36+6;
+ static TMatrixD matHk(nelem,nelem); // vector to mesurement
+ static TMatrixD measR(nelem,nelem); // measurement error
+ static TMatrixD vecZk(nelem,1); // measurement
+ //
+ static TMatrixD vecYk(nelem,1); // Innovation or measurement residual
+ static TMatrixD matHkT(nelem,nelem); // helper matrix Hk transpose
+ static TMatrixD matSk(nelem,nelem); // Innovation (or residual) covariance
+ static TMatrixD matKk(nelem,nelem); // Optimal Kalman gain
+ static TMatrixD mat1(nelem,nelem); // update covariance matrix
+ static TMatrixD covXk2(nelem,nelem); // helper matrix
+ //
+ TMatrixD vecXk(par0); // X vector
+ TMatrixD covXk(cov0); // X covariance
+
+ //
+ //Unit matrix
+ //
+ for (Int_t i=0;i<nelem;i++)
+ for (Int_t j=0;j<nelem;j++){
+ mat1(i,j)=0;
+ if (i==j) mat1(i,j)=1;
+ }
+ matHk = mat1; // unit matrix
+ //
+ vecZk = par1; // measurement
+ measR = cov1; // error of 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
+ //matKk.Print();
+ vecXk += matKk*vecYk; // updated vector
+ covXk2= (mat1-(matKk*matHk));
+ covXk = covXk2*covXk;
+ //
+ CheckCovariance(cov0);
+ CheckCovariance(cov1);
+ CheckCovariance(covXk);
+ //
+ par0 = vecXk; // update measurement param
+ cov0 = covXk; // update measurement covar
+}
+
+
+
+
+Int_t AliTPCcalibAlign::CheckCovariance(TMatrixD &covar){
+ //
+ // check the consistency of covariance matrix
+ //
+ Int_t ncols = covar.GetNcols();
+ Int_t nrows= covar.GetNrows();
+ const Float_t kEpsilon = 0.0001;
+ Int_t nerrors =0;
+ //
+ //
+ //
+ if (nrows!=ncols) {
+ printf("Error 0 - wrong matrix\n");
+ nerrors++;
+ }
+ //
+ // 1. Check that the non diagonal elements
+ //
+ for (Int_t i0=0;i0<nrows;i0++)
+ for (Int_t i1=i0+1;i1<ncols;i1++){
+ Double_t r0 = covar(i0,i1)/TMath::Sqrt(covar(i0,i0)*covar(i1,i1));
+ Double_t r1 = covar(i1,i0)/TMath::Sqrt(covar(i0,i0)*covar(i1,i1));
+ if (TMath::Abs(r0-r1)>kEpsilon){
+ printf("Error 1 - non symetric matrix %d\t%d\t%f",i0,i1,r1-r0);
+ nerrors++;
+ }
+ if (TMath::Abs(r0)>=1){
+ printf("Error 2 - Wrong correlation %d\t%d\t%f\n",i0,i1,r0);
+ nerrors++;
+ }
+ if (TMath::Abs(r1)>=1){
+ printf("Error 3 - Wrong correlation %d\t%d\t%f\n",i0,i1,r1);
+ nerrors++;
+ }
+ }
+ return nerrors;
+}
+
+
+
+void AliTPCcalibAlign::MakeReportDy(TFile *output){
+ //
+ // Draw histogram of dY
+ //
+ Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+
+ AliTPCcalibAlign *align = this;
+ TVectorD vecOOP(36);
+ TVectorD vecOOM(36);
+ TVectorD vecOIP(36);
+ TVectorD vecOIM(36);
+ TVectorD vecOIS(36);
+ TVectorD vecSec(36);
+ TCanvas * cOROCdy = new TCanvas("OROC dy","OROC dy",900,600);
+ cOROCdy->Divide(6,6);
+ TCanvas * cIROCdy = new TCanvas("IROC dy","IROC dy",900,600);
+ cIROCdy->Divide(6,6);
+ TCanvas * cDy = new TCanvas("Dy","Dy",600,700);
+ cDy->Divide(1,2);
+ for (Int_t isec=0;isec<36;isec++){
+ Bool_t isDraw=kFALSE;
+ vecSec(isec)=isec;
+ cOROCdy->cd(isec+1);
+ Int_t secPlus = (isec%18==17)? isec-17:isec+1;
+ Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
+ printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
+ TH1 * hisOOP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus+36);
+ TH1 * hisOOM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus+36);
+ TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec);
+
+ if (hisOIS) {
+ hisOIS = (TH1F*)(hisOIS->Clone());
+ hisOIS->SetDirectory(0);
+ hisOIS->Scale(1./(hisOIS->GetMaximum()+1));
+ hisOIS->SetLineColor(kmicolors[0]);
+ hisOIS->Draw();
+ isDraw = kTRUE;
+ vecOIS(isec)=10*hisOIS->GetMean();
+ }
+ if (hisOOP) {
+ hisOOP = (TH1F*)(hisOOP->Clone());
+ hisOOP->SetDirectory(0);
+ hisOOP->Scale(1./(hisOOP->GetMaximum()+1));
+ hisOOP->SetLineColor(kmicolors[1]);
+ if (isDraw) hisOOP->Draw("same");
+ if (!isDraw) {hisOOP->Draw(""); isDraw=kTRUE;}
+ vecOOP(isec)=10*hisOOP->GetMean();
+ }
+ if (hisOOM) {
+ hisOOM = (TH1F*)(hisOOM->Clone());
+ hisOOM->SetDirectory(0);
+ hisOOM->Scale(1/(hisOOM->GetMaximum()+1));
+ hisOOM->SetLineColor(kmicolors[3]);
+ if (isDraw) hisOOM->Draw("same");
+ if (!isDraw) {hisOOM->Draw(""); isDraw=kTRUE;}
+ vecOOM(isec)=10*hisOOM->GetMean();
+ }
+ }
+ //
+ //
+ for (Int_t isec=0;isec<36;isec++){
+ Bool_t isDraw=kFALSE;
+ cIROCdy->cd(isec+1);
+ Int_t secPlus = (isec%18==17)? isec-17:isec+1;
+ Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
+ printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
+ TH1 * hisOIP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus);
+ TH1 * hisOIM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus);
+ TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec);
+ if (hisOIS) {
+ hisOIS = (TH1F*)(hisOIS->Clone());
+ hisOIS->SetDirectory(0);
+ hisOIS->Scale(1./(hisOIS->GetMaximum()+1));
+ hisOIS->SetLineColor(kmicolors[0]);
+ hisOIS->Draw();
+ isDraw = kTRUE;
+ vecOIS(isec)=10*hisOIS->GetMean();
+ }
+ if (hisOIP) {
+ hisOIP = (TH1F*)(hisOIP->Clone());
+ hisOIP->SetDirectory(0);
+ hisOIP->Scale(1./(hisOIP->GetMaximum()+1));
+ hisOIP->SetLineColor(kmicolors[1]);
+ if (isDraw) hisOIP->Draw("same");
+ if (!isDraw) {hisOIP->Draw(""); isDraw=kTRUE;}
+ hisOIP->Draw("same");
+ vecOIP(isec)=10*hisOIP->GetMean();
+ }
+ if (hisOIM) {
+ hisOIM = (TH1F*)(hisOIM->Clone());
+ hisOIM->SetDirectory(0);
+ hisOIM->Scale(1/(hisOIM->GetMaximum()+1));
+ hisOIM->SetLineColor(kmicolors[3]);
+ if (isDraw) hisOIM->Draw("same");
+ if (!isDraw) {hisOIM->Draw(""); isDraw=kTRUE;}
+ vecOIM(isec)=10*hisOIM->GetMean();
+ }
+ }
+ TGraph* grOIM = new TGraph(36,vecSec.GetMatrixArray(),vecOIM.GetMatrixArray());
+ TGraph* grOIP = new TGraph(36,vecSec.GetMatrixArray(),vecOIP.GetMatrixArray());
+ TGraph* grOIS = new TGraph(36,vecSec.GetMatrixArray(),vecOIS.GetMatrixArray());
+ TGraph* grOOM = new TGraph(36,vecSec.GetMatrixArray(),vecOOM.GetMatrixArray());
+ TGraph* grOOP = new TGraph(36,vecSec.GetMatrixArray(),vecOOP.GetMatrixArray());
+ //
+ grOIS->SetMarkerStyle(kmimarkers[0]);
+ grOIP->SetMarkerStyle(kmimarkers[1]);
+ grOIM->SetMarkerStyle(kmimarkers[3]);
+ grOOP->SetMarkerStyle(kmimarkers[1]);
+ grOOM->SetMarkerStyle(kmimarkers[3]);
+ grOIS->SetMarkerColor(kmicolors[0]);
+ grOIP->SetMarkerColor(kmicolors[1]);
+ grOIM->SetMarkerColor(kmicolors[3]);
+ grOOP->SetMarkerColor(kmicolors[1]);
+ grOOM->SetMarkerColor(kmicolors[3]);
+ grOIS->SetLineColor(kmicolors[0]);
+ grOIP->SetLineColor(kmicolors[1]);
+ grOIM->SetLineColor(kmicolors[3]);
+ grOOP->SetLineColor(kmicolors[1]);
+ grOOM->SetLineColor(kmicolors[3]);
+ grOIS->SetMaximum(1.5);
+ grOIS->SetMinimum(-1.5);
+ grOIS->GetXaxis()->SetTitle("Sector number");
+ grOIS->GetYaxis()->SetTitle("#Delta_{y} (mm)");
+
+ cDy->cd(1);
+ grOIS->Draw("apl");
+ grOIM->Draw("pl");
+ grOIP->Draw("pl");
+ cDy->cd(2);
+ grOIS->Draw("apl");
+ grOOM->Draw("pl");
+ grOOP->Draw("pl");
+ cOROCdy->SaveAs("picAlign/OROCOROCdy.eps");
+ cOROCdy->SaveAs("picAlign/OROCOROCdy.gif");
+ cOROCdy->SaveAs("picAlign/OROCOROCdy.root");
+ //
+ cIROCdy->SaveAs("picAlign/OROCIROCdy.eps");
+ cIROCdy->SaveAs("picAlign/OROCIROCdy.gif");
+ cIROCdy->SaveAs("picAlign/OROCIROCdy.root");
+ //
+ cDy->SaveAs("picAlign/Sectordy.eps");
+ cDy->SaveAs("picAlign/Sectordy.gif");
+ cDy->SaveAs("picAlign/Sectordy.root");
+ if (output){
+ output->cd();
+ cOROCdy->Write("OROCOROCDy");
+ cIROCdy->Write("OROCIROCDy");
+ cDy->Write("SectorDy");
+ }
+}
+
+void AliTPCcalibAlign::MakeReportDyPhi(TFile */*output*/){
+ //
+ //
+ //
+ Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+
+ AliTPCcalibAlign *align = this;
+ TCanvas * cOROCdyPhi = new TCanvas("OROC dyphi","OROC dyphi",900,600);
+ cOROCdyPhi->Divide(6,6);
+ for (Int_t isec=0;isec<36;isec++){
+ cOROCdyPhi->cd(isec+1);
+ Int_t secPlus = (isec%18==17)? isec-17:isec+1;
+ Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
+ //printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
+ TH2F *htemp=0;
+ TProfile * profdyphiOOP=0,*profdyphiOOM=0,*profdyphiOOS=0;
+ htemp = (TH2F*) (align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secPlus+36));
+ if (htemp) profdyphiOOP= htemp->ProfileX();
+ htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secMinus+36));
+ if (htemp) profdyphiOOM= htemp->ProfileX();
+ htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,isec));
+ if (htemp) profdyphiOOS= htemp->ProfileX();
+
+ if (profdyphiOOS){
+ profdyphiOOS->SetLineColor(kmicolors[0]);
+ profdyphiOOS->SetMarkerStyle(kmimarkers[0]);
+ profdyphiOOS->SetMarkerSize(0.2);
+ profdyphiOOS->SetMaximum(0.5);
+ profdyphiOOS->SetMinimum(-0.5);
+ profdyphiOOS->SetXTitle("tan(#phi)");
+ profdyphiOOS->SetYTitle("#DeltaY (cm)");
+ }
+ if (profdyphiOOP){
+ profdyphiOOP->SetLineColor(kmicolors[1]);
+ profdyphiOOP->SetMarkerStyle(kmimarkers[1]);
+ profdyphiOOP->SetMarkerSize(0.2);
+ }
+ if (profdyphiOOM){
+ profdyphiOOM->SetLineColor(kmicolors[3]);
+ profdyphiOOM->SetMarkerStyle(kmimarkers[3]);
+ profdyphiOOM->SetMarkerSize(0.2);
+ }
+ if (profdyphiOOS){
+ profdyphiOOS->Draw();
+ }else{
+ if (profdyphiOOM) profdyphiOOM->Draw("");
+ if (profdyphiOOP) profdyphiOOP->Draw("");
+ }
+ if (profdyphiOOM) profdyphiOOM->Draw("same");
+ if (profdyphiOOP) profdyphiOOP->Draw("same");
+
+ }
+}
+