//---- Root headers --------
#include <TTree.h>
+#include <TMatrixD.h>
//---- AliRoot headers -----
#include "AliStrLine.h"
#include "AliVertexerTracks.h"
//----------------------------------------------------------------------------
AliVertexerTracks::AliVertexerTracks():
- TObject(),fVert()
+TObject(),
+fVert(),
+fCurrentVertex(0),
+fMinTracks(2),
+fMinITSClusters(5),
+fTrkArray(),
+fTrksToSkip(0),
+fNTrksToSkip(0),
+fDCAcut(0),
+fAlgo(1),
+fNSigma(3),
+fDebug(0)
{
//
// Default constructor
//
SetVtxStart();
+ SetVtxStartSigma();
SetMinTracks();
- fDCAcut=0;
- fAlgo=1;
+ SetMinITSClusters();
+ SetNSigmad0();
}
//-----------------------------------------------------------------------------
AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
- TObject(),fVert()
+TObject(),
+fVert(),
+fCurrentVertex(0),
+fMinTracks(2),
+fMinITSClusters(5),
+fTrkArray(),
+fTrksToSkip(0),
+fNTrksToSkip(0),
+fDCAcut(0),
+fAlgo(1),
+fNSigma(3),
+fDebug(0)
{
//
// Standard constructor
//
SetVtxStart(xStart,yStart);
+ SetVtxStartSigma();
SetMinTracks();
- fDCAcut=0;
- fAlgo=1;
+ SetMinITSClusters();
+ SetNSigmad0();
}
//-----------------------------------------------------------------------------
AliVertexerTracks::~AliVertexerTracks() {
// Default Destructor
- // The objects poited by the following pointers are not owned
+ // The objects pointed by the following pointer are not owned
// by this class and are not deleted
+ fCurrentVertex = 0;
+ delete[] fTrksToSkip;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
+//
+// Mark the tracks not ot be used in the vertex finding
+//
+ fNTrksToSkip = n;
+ fTrksToSkip = new Int_t[n];
+ for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
+ return;
}
//----------------------------------------------------------------------------
-Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree) {
+Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree, Int_t OptImpParCut) {
//
// Propagate tracks to initial vertex position and store them in a TObjArray
//
- Double_t alpha,xlStart;
+ Double_t maxd0rphi = 3.;
Int_t nTrks = 0;
+ Double_t sigmaCurr[3],sigmaVtx[3],posVtx[3];
+ Double_t covVtx[6],covVtxXY;
+ Float_t d0z0[2],covd0z0[3];
+ Double_t momentum[3],postrk[3];
+ Double_t pt,sigma,sigmavtxphi,phitrk;
+ Double_t field=GetField();
+
+ AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalSigma);
Int_t nEntries = (Int_t)trkTree.GetEntries();
if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
fTrkArray.Expand(nEntries);
+ if(fDebug) {
+ printf(" PrepareTracks()\n");
+ }
+
for(Int_t i=0; i<nEntries; i++) {
AliESDtrack *track = new AliESDtrack;
trkTree.SetBranchAddress("tracks",&track);
trkTree.GetEvent(i);
- // propagate track to vtxSeed
- alpha = track->GetAlpha();
- xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
- track->PropagateTo(xlStart,GetField()); // to vtxSeed
+
+
+ // propagate track to vertex
+ if(OptImpParCut==1) { // OptImpParCut==1
+ track->RelateToVertex(initVertex,field,100.);
+ initVertex->GetXYZ(posVtx);
+ initVertex->GetSigmaXYZ(sigmaVtx);
+ covVtxXY = 0.;
+ } else { // OptImpParCut==2
+ fCurrentVertex->GetSigmaXYZ(sigmaCurr);
+ if((sigmaCurr[0]+sigmaCurr[1])<(fNominalSigma[0]+fNominalSigma[1])) {
+ track->RelateToVertex(fCurrentVertex,field,100.);
+ fCurrentVertex->GetXYZ(posVtx);
+ fCurrentVertex->GetSigmaXYZ(sigmaVtx);
+ fCurrentVertex->GetCovMatrix(covVtx);
+ covVtxXY = covVtx[1];
+ } else {
+ track->RelateToVertex(initVertex,field,100.);
+ initVertex->GetXYZ(posVtx);
+ initVertex->GetSigmaXYZ(sigmaVtx);
+ covVtxXY = 0.;
+ }
+ }
+
+ // select tracks with d0rphi < maxd0rphi
+ track->GetImpactParameters(d0z0,covd0z0);
+
+ track->GetXYZ(postrk);
+ track->GetPxPyPz(momentum);
+ pt = TMath::Sqrt(momentum[0]*momentum[0]+momentum[1]*momentum[1]);
+
+ phitrk = TMath::ATan2(postrk[1]-posVtx[1],postrk[0]-posVtx[0]);
+ sigmavtxphi = TMath::Sqrt(sigmaVtx[0]*sigmaVtx[0]*
+ TMath::Cos(phitrk)*TMath::Cos(phitrk)+
+ sigmaVtx[1]*sigmaVtx[1]*
+ TMath::Sin(phitrk)*TMath::Sin(phitrk)+
+ covVtxXY*
+ TMath::Cos(phitrk)*TMath::Sin(phitrk));
+ sigma = TMath::Sqrt(Sigmad0rphi(pt)*Sigmad0rphi(pt)+
+ sigmavtxphi*sigmavtxphi);
+ maxd0rphi = fNSigma*sigma;
+ if(OptImpParCut==1) maxd0rphi *= 5.;
+
+ if(fDebug) printf("trk %d; lab %d; |d0| = %f; cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi);
+ if(TMath::Abs(d0z0[0]) > maxd0rphi) {
+ if(fDebug) printf(" rejected\n");
+ delete track; continue;
+ }
+
fTrkArray.AddLast(track);
nTrks++;
}
+
+ delete initVertex;
+
return nTrks;
}
//----------------------------------------------------------------------------
+Double_t AliVertexerTracks::Sigmad0rphi(Double_t pt) const {
+//
+// Impact parameter resolution in rphi [cm] vs pt [GeV/c]
+//
+ Double_t a_meas = 11.6;
+ Double_t a_scatt = 65.8;
+ Double_t b = 1.878;
+
+ Double_t sigma = a_meas*a_meas+a_scatt*a_scatt/TMath::Power(pt,b);
+ sigma = 0.0001*TMath::Sqrt(sigma);
+
+ return sigma;
+}
+//----------------------------------------------------------------------------
AliVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree) {
//
// Return vertex from tracks in trkTree
//
// get tracks and propagate them to initial vertex position
- Int_t nTrks = PrepareTracks(*trkTree);
+ Int_t nTrks = PrepareTracks(*trkTree,1);
if(nTrks < fMinTracks) {
printf("TooFewTracks\n");
Double_t vtx[3]={0,0,0};
if(fAlgo==5) VertexFinder(0);
return &fVert;
}
-
//----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindVertex(const AliESD *event) {
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
+{
//
-// This is a simple wrapping (by Jouri.Belikov@cern.ch) over the original
-// code by the authors of this class.
+// Primary vertex for current ESD event
+// (Two iterations:
+// 1st with 5*fNSigma*sigma(pt) cut w.r.t. to initial vertex;
+// 2nd with fNSigma*sigma(pt) cut w.r.t. to vertex found in 1st iteration)
//
- Int_t nt=event->GetNumberOfTracks(), nacc=0;
- while (nt--) {
- AliESDtrack *t=event->GetTrack(nt);
- if ((t->GetStatus()&AliESDtrack::kITSrefit)==0) continue;
- fTrkArray.AddLast(t);
- nacc++;
+ fCurrentVertex = 0;
+
+ Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
+ TTree *trkTree = new TTree("TreeT","tracks");
+ AliESDtrack *esdTrack = 0;
+ trkTree->Branch("tracks","AliESDtrack",&esdTrack);
+
+ Bool_t skipThis;
+ for(Int_t i=0; i<entr; i++) {
+ // check tracks to skip
+ skipThis = kFALSE;
+ for(Int_t j=0; j<fNTrksToSkip; j++) {
+ if(i==fTrksToSkip[j]) {
+ if(fDebug) printf("skipping track: %d\n",i);
+ skipThis = kTRUE;
+ }
+ }
+ if(skipThis) continue;
+ AliESDtrack *et = esdEvent->GetTrack(i);
+ esdTrack = new AliESDtrack(*et);
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue;
+ Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
+ if(nclus<fMinITSClusters) continue;
+ trkTree->Fill();
+ }
+ delete esdTrack;
+
+ // ITERATION 1
+ // propagate tracks to initVertex
+ // preselect them (reject for |d0|>5*fNSigma*sigma w.r.t. initVertex)
+ Int_t nTrks;
+ nTrks = PrepareTracks(*trkTree,1);
+ if(fDebug) printf(" tracks prepared - step 1: %d\n",nTrks);
+ if(nTrks < fMinTracks) {
+ printf("TooFewTracks\n");
+ TooFewTracks(esdEvent);
+ if(fDebug) fCurrentVertex->PrintStatus();
+ return fCurrentVertex;
}
- // get tracks and propagate them to initial vertex position
- if(nacc < fMinTracks) {
+ // vertex finder
+ switch (fAlgo) {
+ case 1: StrLinVertexFinderMinDist(1); break;
+ case 2: StrLinVertexFinderMinDist(0); break;
+ case 3: HelixVertexFinder(); break;
+ case 4: VertexFinder(1); break;
+ case 5: VertexFinder(0); break;
+ default: printf("Wrong algorithm\n"); break;
+ }
+ if(fDebug) printf(" vertex finding completed\n");
+
+ // vertex fitter
+ VertexFitter(kTRUE);
+ if(fDebug) printf(" vertex fit completed\n");
+
+ // ITERATION 2
+ // propagate tracks to best between initVertex and fCurrentVertex
+ // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
+ // between initVertex and fCurrentVertex)
+ nTrks = PrepareTracks(*trkTree,2);
+ delete trkTree;
+ if(fDebug) printf(" tracks prepared - step 2: %d\n",nTrks);
+ if(nTrks < fMinTracks) {
printf("TooFewTracks\n");
- Double_t vtx[3]={0,0,0};
- fVert.SetXYZ(vtx);
- fVert.SetDispersion(999);
- fVert.SetNContributors(-5);
- } else
- switch (fAlgo) {
- case 1: StrLinVertexFinderMinDist(1); break;
- case 2: StrLinVertexFinderMinDist(0); break;
- case 3: HelixVertexFinder(); break;
- case 4: VertexFinder(1); break;
- case 5: VertexFinder(0); break;
- default: printf("Wrong algorithm\n"); break;
- }
+ TooFewTracks(esdEvent);
+ if(fDebug) fCurrentVertex->PrintStatus();
+ return fCurrentVertex;
+ }
+
+ // vertex finder
+ switch (fAlgo) {
+ case 1: StrLinVertexFinderMinDist(1); break;
+ case 2: StrLinVertexFinderMinDist(0); break;
+ case 3: HelixVertexFinder(); break;
+ case 4: VertexFinder(1); break;
+ case 5: VertexFinder(0); break;
+ default: printf("Wrong algorithm\n"); break;
+ }
+ if(fDebug) printf(" vertex finding completed\n");
+
+ // fitter
+ VertexFitter(kTRUE);
+ if(fDebug) printf(" vertex fit completed\n");
+
fTrkArray.Clear();
- return &fVert;
+
+ // take true pos from ESD
+ Double_t tp[3];
+ esdEvent->GetVertex()->GetTruePos(tp);
+ fCurrentVertex->SetTruePos(tp);
+ if(fNominalSigma[0]>1.) {
+ fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
+ } else {
+ fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
+ }
+
+ if(fDebug) fCurrentVertex->PrintStatus();
+ if(fTrksToSkip) delete [] fTrksToSkip;
+
+ return fCurrentVertex;
}
+//----------------------------------------------------------------------------
+AliESDVertex* AliVertexerTracks::FindPrimaryVertexOld(const AliESD *esdEvent)
+{
+//
+// Primary vertex for current ESD event (one iteration with |d0rphi|<3cm)
+//
+ fCurrentVertex = 0;
+
+ Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
+ TTree *trkTree = new TTree("TreeT","tracks");
+ AliESDtrack *esdTrack = 0;
+ trkTree->Branch("tracks","AliESDtrack",&esdTrack);
+
+ for(Int_t i=0; i<entr; i++) {
+ AliESDtrack *et = esdEvent->GetTrack(i);
+ esdTrack = new AliESDtrack(*et);
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
+ if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue;
+ Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
+ if(nclus<fMinITSClusters) continue;
+
+ trkTree->Fill();
+ }
+ delete esdTrack;
+
+ // preselect tracks and propagate them to initial vertex position
+ Int_t nTrks = PrepareTracks(*trkTree,1);
+ delete trkTree;
+ if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+ if(nTrks < fMinTracks) {
+ printf("TooFewTracks\n");
+ fCurrentVertex = new AliESDVertex(0.,0.,-1);
+ return fCurrentVertex;
+ }
+ // VERTEX FINDER
+ switch (fAlgo) {
+ case 1: StrLinVertexFinderMinDist(1); break;
+ case 2: StrLinVertexFinderMinDist(0); break;
+ case 3: HelixVertexFinder(); break;
+ case 4: VertexFinder(1); break;
+ case 5: VertexFinder(0); break;
+ default: printf("Wrong algorithm\n"); break;
+ }
+
+ // VERTEX FITTER
+ VertexFitter();
+ if(fDebug) printf(" vertex fit completed\n");
+
+
+ Double_t tp[3];
+ esdEvent->GetVertex()->GetTruePos(tp);
+ fCurrentVertex->SetTruePos(tp);
+ fCurrentVertex->SetTitle("VertexerTracks");
+ fTrkArray.Clear();
+ return fCurrentVertex;
+}
//----------------------------------------------------------------------------
AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
//
Int_t ncombi = 0;
AliESDtrack *track1;
AliESDtrack *track2;
+ Double_t pos[3],dir[3];
Double_t alpha,mindist;
Double_t field=GetField();
track1 = (AliESDtrack*)fTrkArray.At(i);
alpha=track1->GetAlpha();
mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+ track1->GetXYZAt(mindist,field,pos);
+ track1->GetPxPyPzAt(mindist,field,dir);
+ AliStrLine *line1 = new AliStrLine(pos,dir);
- Double_t pos[3]; track1->GetXYZAt(mindist,field,pos);
- Double_t dir[3]; track1->GetPxPyPzAt(mindist,field,dir);
- AliStrLine *line1 = new AliStrLine(pos,dir);
+ // AliStrLine *line1 = new AliStrLine();
+ // track1->ApproximateHelixWithLine(mindist,field,line1);
for(Int_t j=i+1; j<nacc; j++){
track2 = (AliESDtrack*)fTrkArray.At(j);
alpha=track2->GetAlpha();
mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-
- Double_t pos[3]; track2->GetXYZAt(mindist,field,pos);
- Double_t dir[3]; track2->GetPxPyPzAt(mindist,field,dir);
- AliStrLine *line2 = new AliStrLine(pos,dir);
-
+ track2->GetXYZAt(mindist,field,pos);
+ track2->GetPxPyPzAt(mindist,field,dir);
+ AliStrLine *line2 = new AliStrLine(pos,dir);
+ // AliStrLine *line2 = new AliStrLine();
+ // track2->ApproximateHelixWithLine(mindist,field,line2);
Double_t distCA=line2->GetDCA(line1);
if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
Double_t pnt1[3],pnt2[3],crosspoint[3];
track1 = (AliESDtrack*)fTrkArray.At(i);
Double_t alpha=track1->GetAlpha();
Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-
- Double_t pos[3]; track1->GetXYZAt(mindist,field,pos);
- Double_t dir[3]; track1->GetPxPyPzAt(mindist,field,dir);
- AliStrLine *line1 = new AliStrLine(pos,dir);
+ Double_t pos[3],dir[3];
+ track1->GetXYZAt(mindist,field,pos);
+ track1->GetPxPyPzAt(mindist,field,dir);
+ AliStrLine *line1 = new AliStrLine(pos,dir);
+ // AliStrLine *line1 = new AliStrLine();
+ // track1->ApproximateHelixWithLine(mindist,field,line1);
Double_t p0[3],cd[3];
line1->GetP0(p0);
return det;
}
//____________________________________________________________________________
-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d){
+void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d){
//
Double_t x12=p0[0]-p1[0];
}
//____________________________________________________________________________
-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d){
+void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d){
//
Double_t x12=p1[0]-p0[0];
Double_t y12=p1[1]-p0[1];
d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
-}
+ }
//_____________________________________________________________________________
Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
//
Double_t z10=p0[2]-x0[2];
return ((x10*x10+y10*y10+z10*z10)*(x12*x12+y12*y12+z12*z12)-(x10*x12+y10*y12+z10*z12)*(x10*x12+y10*y12+z10*z12))/(x12*x12+y12*y12+z12*z12);
}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
+//
+// The optimal estimate of the vertex position is given by a "weighted
+// average of tracks positions"
+// Original method: CMS Note 97/0051
+//
+ Double_t initPos[3];
+ fVert.GetXYZ(initPos);
+ if(fDebug) {
+ printf(" VertexFitter(): start\n");
+ printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
+ printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
+ printf("Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
+ if(useNominalVtx) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],fNominalSigma[0],fNominalPos[1],fNominalSigma[1]);
+ }
+
+ Int_t i,j,k,step=0;
+ TMatrixD rv(3,1);
+ TMatrixD vV(3,3);
+ rv(0,0) = initPos[0];
+ rv(1,0) = initPos[1];
+ rv(2,0) = 0.;
+ Double_t xlStart,alpha;
+ Double_t rotAngle;
+ Double_t cosRot,sinRot;
+ Double_t cc[15];
+ Int_t nUsedTrks;
+ Double_t chi2,chi2i;
+ Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
+ AliESDtrack *t = 0;
+ Int_t failed = 0;
+
+ // 2 steps:
+ // 1st - estimate of vtx using all tracks
+ // 2nd - estimate of global chi2
+ for(step=0; step<2; step++) {
+ if(fDebug) printf(" step = %d\n",step);
+ chi2 = 0.;
+ nUsedTrks = 0;
+
+ TMatrixD sumWiri(3,1);
+ TMatrixD sumWi(3,3);
+ for(i=0; i<3; i++) {
+ sumWiri(i,0) = 0.;
+ for(j=0; j<3; j++) sumWi(j,i) = 0.;
+ }
+
+ if(useNominalVtx) {
+ for(i=0; i<3; i++) {
+ sumWiri(i,0) += fNominalPos[i]/fNominalSigma[i]/fNominalSigma[i];
+ sumWi(i,i) += 1./fNominalSigma[i]/fNominalSigma[i];
+ }
+ }
+
+
+ // loop on tracks
+ for(k=0; k<arrEntries; k++) {
+ // get track from track array
+ t = (AliESDtrack*)fTrkArray.At(k);
+ alpha = t->GetAlpha();
+ xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
+ t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed
+ rotAngle = alpha;
+ if(alpha<0.) rotAngle += 2.*TMath::Pi();
+ cosRot = TMath::Cos(rotAngle);
+ sinRot = TMath::Sin(rotAngle);
+
+ // vector of track global coordinates
+ TMatrixD ri(3,1);
+ ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
+ ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
+ ri(2,0) = t->GetZ();
+
+ // matrix to go from global (x,y,z) to local (y,z);
+ TMatrixD qQi(2,3);
+ qQi(0,0) = -sinRot;
+ qQi(0,1) = cosRot;
+ qQi(0,2) = 0.;
+ qQi(1,0) = 0.;
+ qQi(1,1) = 0.;
+ qQi(1,2) = 1.;
+
+ // covariance matrix of local (y,z) - inverted
+ TMatrixD uUi(2,2);
+ t->GetExternalCovariance(cc);
+ uUi(0,0) = cc[0];
+ uUi(0,1) = cc[1];
+ uUi(1,0) = cc[1];
+ uUi(1,1) = cc[2];
+
+ // weights matrix: wWi = qQiT * uUiInv * qQi
+ if(uUi.Determinant() <= 0.) continue;
+ TMatrixD uUiInv(TMatrixD::kInverted,uUi);
+ TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
+ TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
+
+ // track chi2
+ TMatrixD deltar = rv; deltar -= ri;
+ TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
+ chi2i = deltar(0,0)*wWideltar(0,0)+
+ deltar(1,0)*wWideltar(1,0)+
+ deltar(2,0)*wWideltar(2,0);
+
+
+ // add to total chi2
+ chi2 += chi2i;
+
+ TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
+
+ sumWiri += wWiri;
+ sumWi += wWi;
+
+ nUsedTrks++;
+ } // end loop on tracks
+
+ if(nUsedTrks < fMinTracks) {
+ failed=1;
+ continue;
+ }
+
+ Double_t determinant = sumWi.Determinant();
+ //cerr<<" determinant: "<<determinant<<endl;
+ if(determinant < 100.) {
+ printf("det(V) = 0\n");
+ failed=1;
+ continue;
+ }
+
+ // inverted of weights matrix
+ TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
+ vV = invsumWi;
+
+ // position of primary vertex
+ rv.Mult(vV,sumWiri);
+
+ } // end loop on the 2 steps
+
+ delete t;
+
+ if(failed) {
+ printf("TooFewTracks\n");
+ fCurrentVertex = new AliESDVertex(0.,0.,-1);
+ return;
+ }
+
+ Double_t position[3];
+ position[0] = rv(0,0);
+ position[1] = rv(1,0);
+ position[2] = rv(2,0);
+ Double_t covmatrix[6];
+ covmatrix[0] = vV(0,0);
+ covmatrix[1] = vV(0,1);
+ covmatrix[2] = vV(1,1);
+ covmatrix[3] = vV(0,2);
+ covmatrix[4] = vV(1,2);
+ covmatrix[5] = vV(2,2);
+
+ // store data in the vertex object
+ fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
+
+ if(fDebug) {
+ printf(" VertexFitter(): finish\n");
+ printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
+ fCurrentVertex->PrintStatus();
+ }
+
+ return;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) {
+//
+// When the number of tracks is < fMinTracks
+//
+
+ // deal with vertices not found
+ Double_t pos[3],err[3];
+ Int_t ncontr=0;
+ pos[0] = fNominalPos[0];
+ err[0] = fNominalSigma[0];
+ pos[1] = fNominalPos[1];
+ err[1] = fNominalSigma[1];
+ pos[2] = esdEvent->GetVertex()->GetZv();
+ err[2] = esdEvent->GetVertex()->GetZRes();
+ if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()<=0)
+ ncontr = -1; // (x,y,z) = (0,0,0)
+ if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()>0)
+ ncontr = -2; // (x,y,z) = (0,0,z_fromSPD)
+ if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()<=0)
+ ncontr = -3; // (x,y,z) = (x_mean,y_mean,0)
+ if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()>0)
+ ncontr = -4; // (x,y,z) = (x_mean,y_mean,z_fromSPD)
+ fCurrentVertex = 0;
+ fCurrentVertex = new AliESDVertex(pos,err);
+ fCurrentVertex->SetNContributors(ncontr);
+
+ Double_t tp[3];
+ esdEvent->GetVertex()->GetTruePos(tp);
+ fCurrentVertex->SetTruePos(tp);
+ fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
+ if(ncontr==-1||ncontr==-2)
+ fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
+
+ return;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff