+AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
+ Float_t *diamondxyz,
+ Float_t *diamondcov) const
+{
+//
+// Removes diamond constraint from fit of inVtx
+//
+
+ if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
+ printf("ERROR: primary vertex has no constraint: cannot remove it\n");
+ return 0x0;
+ }
+ if(inVtx->GetNContributors()<3) {
+ printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
+ return 0x0;
+ }
+
+ // diamond constraint
+ TMatrixD vVb(3,3);
+ vVb(0,0) = diamondcov[0];
+ vVb(0,1) = diamondcov[1];
+ vVb(0,2) = 0.;
+ vVb(1,0) = diamondcov[1];
+ vVb(1,1) = diamondcov[2];
+ vVb(1,2) = 0.;
+ vVb(2,0) = 0.;
+ vVb(2,1) = 0.;
+ vVb(2,2) = diamondcov[5];
+ TMatrixD vVbInv(TMatrixD::kInverted,vVb);
+ TMatrixD rb(3,1);
+ rb(0,0) = diamondxyz[0];
+ rb(1,0) = diamondxyz[1];
+ rb(2,0) = diamondxyz[2];
+ TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
+
+ // input vertex
+ TMatrixD rv(3,1);
+ rv(0,0) = inVtx->GetXv();
+ rv(1,0) = inVtx->GetYv();
+ rv(2,0) = inVtx->GetZv();
+ TMatrixD vV(3,3);
+ Double_t cov[6];
+ inVtx->GetCovMatrix(cov);
+ vV(0,0) = cov[0];
+ vV(0,1) = cov[1]; vV(1,0) = cov[1];
+ vV(1,1) = cov[2];
+ vV(0,2) = cov[3]; vV(2,0) = cov[3];
+ vV(1,2) = cov[4]; vV(2,1) = cov[4];
+ vV(2,2) = cov[5];
+ TMatrixD vVInv(TMatrixD::kInverted,vV);
+ TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
+
+
+ TMatrixD sumWi = vVInv - vVbInv;
+
+
+ TMatrixD sumWiri = vVInvrv - vVbInvrb;
+
+ TMatrixD rvnew(3,1);
+ TMatrixD vVnew(3,3);
+
+ // new inverted of weights matrix
+ TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
+ vVnew = invsumWi;
+ // new position of primary vertex
+ rvnew.Mult(vVnew,sumWiri);
+
+ Double_t position[3];
+ position[0] = rvnew(0,0);
+ position[1] = rvnew(1,0);
+ position[2] = rvnew(2,0);
+ cov[0] = vVnew(0,0);
+ cov[1] = vVnew(0,1);
+ cov[2] = vVnew(1,1);
+ cov[3] = vVnew(0,2);
+ cov[4] = vVnew(1,2);
+ cov[5] = vVnew(2,2);
+
+
+ Double_t chi2 = inVtx->GetChi2();
+
+ // diamond constribution to chi2
+ TMatrixD deltar = rv; deltar -= rb;
+ TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
+ Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
+ deltar(1,0)*vVbInvdeltar(1,0)+
+ deltar(2,0)*vVbInvdeltar(2,0);
+ // remove from total chi2
+ chi2 -= chi2b;
+
+ // store data in the vertex object
+ AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
+ outVtx->SetTitle("VertexerTracksNoConstraint");
+ UShort_t *inindices = inVtx->GetIndices();
+ Int_t nIndices = inVtx->GetNIndices();
+ outVtx->SetIndices(nIndices,inindices);
+
+ return outVtx;
+}
+//---------------------------------------------------------------------------