+//----------------------------------------------------------------------------
+Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
+ Double_t xToGo) {
+ //----------------------------------------------------------------
+ // COPIED from AliTracker
+ //
+ // Propagates the track to the plane X=xk (cm).
+ //
+ // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
+ //----------------------------------------------------------------
+
+ const Double_t kEpsilon = 0.00001;
+ Double_t xpos = track->GetX();
+ Double_t dir = (xpos<xToGo) ? 1. : -1.;
+ Double_t maxStep = 5;
+ Double_t maxSnp = 0.8;
+ //
+ while ( (xToGo-xpos)*dir > kEpsilon){
+ Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
+ Double_t x = xpos+step;
+ Double_t xyz0[3],xyz1[3];
+ track->GetXYZ(xyz0); //starting global position
+
+ if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE; // no prolongation
+ xyz1[2]+=kEpsilon; // waiting for bug correction in geo
+
+ if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
+ if(!track->PropagateTo(x,GetFieldkG())) return kFALSE;
+
+ if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
+ track->GetXYZ(xyz0); // global position
+ Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
+ //
+ Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
+ sa=TMath::Sin(alphan-track->GetAlpha());
+ Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
+ Double_t sinNew = sf*ca - cf*sa;
+ if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
+ if(!track->Rotate(alphan)) return kFALSE;
+
+ xpos = track->GetX();
+ }
+ return kTRUE;
+}
+//---------------------------------------------------------------------------
+AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
+ const TObjArray *trkArray,
+ UShort_t *id,
+ const Float_t *diamondxy) const
+{
+//
+// Removes tracks in trksTree from fit of inVtx
+//
+
+ if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
+ printf("ERROR: primary vertex has no constraint: cannot remove tracks");
+ return 0x0;
+ }
+ if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
+ printf("WARNING: result of tracks' removal will be only approximately correct");
+
+ 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 sumWi(TMatrixD::kInverted,vV);
+ TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
+
+ Int_t nUsedTrks = inVtx->GetNIndices();
+ Double_t chi2 = inVtx->GetChi2();
+
+ AliExternalTrackParam *track = 0;
+ Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
+ for(Int_t i=0;i<ntrks;i++) {
+ track = (AliExternalTrackParam*)trkArray->At(i);
+ if(!inVtx->UsesTrack(id[i])) {
+ printf("track %d was not used in vertex fit",id[i]);
+ continue;
+ }
+ Double_t alpha = track->GetAlpha();
+ Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
+ track->PropagateTo(xl,GetFieldkG());
+ // vector of track global coordinates
+ TMatrixD ri(3,1);
+ // covariance matrix of ri
+ TMatrixD wWi(3,3);
+
+ // get space point from track
+ if(!TrackToPoint(track,ri,wWi)) continue;
+
+ TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
+
+ sumWi -= wWi;
+ sumWiri -= wWiri;
+
+ // track contribution to chi2
+ TMatrixD deltar = rv; deltar -= ri;
+ TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
+ Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
+ deltar(1,0)*wWideltar(1,0)+
+ deltar(2,0)*wWideltar(2,0);
+ // remove from total chi2
+ chi2 -= chi2i;
+
+ nUsedTrks--;
+ if(nUsedTrks<2) {
+ printf("Trying to remove too many tracks!");
+ return 0x0;
+ }
+ }
+
+ 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);
+
+ // store data in the vertex object
+ AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
+ outVtx->SetTitle(inVtx->GetTitle());
+ UShort_t *inindices = inVtx->GetIndices();
+ Int_t nIndices = nUsedTrks;
+ UShort_t *outindices = new UShort_t[nIndices];
+ Int_t j=0;
+ for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
+ Bool_t copyindex=kTRUE;
+ for(Int_t l=0; l<ntrks; l++) {
+ if(inindices[k]==id[l]) copyindex=kFALSE;
+ }
+ if(copyindex) {
+ outindices[j] = inindices[k]; j++;
+ }
+ }
+ outVtx->SetIndices(nIndices,outindices);
+ if (outindices) delete [] outindices;
+
+ /*
+ printf("Vertex before removing tracks:");
+ inVtx->PrintStatus();
+ inVtx->PrintIndices();
+ printf("Vertex after removing tracks:");
+ outVtx->PrintStatus();
+ outVtx->PrintIndices();
+ */
+
+ return outVtx;
+}
+//---------------------------------------------------------------------------
+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;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetCuts(Double_t *cuts)
+{
+//
+// Cut values
+//
+ SetDCAcut(cuts[0]);
+ SetDCAcutIter0(cuts[1]);
+ SetMaxd0z0(cuts[2]);
+ if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
+ SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
+ SetMinTracks((Int_t)(cuts[4]));
+ SetNSigmad0(cuts[5]);
+ SetMinDetFitter(cuts[6]);
+ SetMaxTgl(cuts[7]);
+ SetFiducialRZ(cuts[8],cuts[9]);
+ fAlgo=(Int_t)(cuts[10]);
+ fAlgoIter0=(Int_t)(cuts[11]);
+
+ return;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetITSMode(Double_t dcacut,
+ Double_t dcacutIter0,
+ Double_t maxd0z0,
+ Int_t minCls,
+ Int_t mintrks,
+ Double_t nsigma,
+ Double_t mindetfitter,
+ Double_t maxtgl,
+ Double_t fidR,
+ Double_t fidZ,
+ Int_t finderAlgo,
+ Int_t finderAlgoIter0)
+{
+//
+// Cut values for ITS mode
+//
+ fMode = 0;
+ if(minCls>0) {
+ SetITSrefitRequired();
+ } else {
+ SetITSrefitNotRequired();
+ }
+ SetDCAcut(dcacut);
+ SetDCAcutIter0(dcacutIter0);
+ SetMaxd0z0(maxd0z0);
+ SetMinClusters(TMath::Abs(minCls));
+ SetMinTracks(mintrks);
+ SetNSigmad0(nsigma);
+ SetMinDetFitter(mindetfitter);
+ SetMaxTgl(maxtgl);
+ SetFiducialRZ(fidR,fidZ);
+ fAlgo=finderAlgo;
+ fAlgoIter0=finderAlgoIter0;
+
+ return;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetTPCMode(Double_t dcacut,
+ Double_t dcacutIter0,
+ Double_t maxd0z0,
+ Int_t minCls,
+ Int_t mintrks,
+ Double_t nsigma,
+ Double_t mindetfitter,
+ Double_t maxtgl,
+ Double_t fidR,
+ Double_t fidZ,
+ Int_t finderAlgo,
+ Int_t finderAlgoIter0)
+{
+//
+// Cut values for TPC mode
+//
+ fMode = 1;
+ SetITSrefitNotRequired();
+ SetDCAcut(dcacut);
+ SetDCAcutIter0(dcacutIter0);
+ SetMaxd0z0(maxd0z0);
+ SetMinClusters(minCls);
+ SetMinTracks(mintrks);
+ SetNSigmad0(nsigma);
+ SetMinDetFitter(mindetfitter);
+ SetMaxTgl(maxtgl);
+ SetFiducialRZ(fidR,fidZ);
+ fAlgo=finderAlgo;
+ fAlgoIter0=finderAlgoIter0;
+
+ return;
+}