//---- AliRoot headers -----
#include "AliStrLine.h"
#include "AliVertexerTracks.h"
-#include "AliESD.h"
+#include "AliESDEvent.h"
#include "AliESDtrack.h"
ClassImp(AliVertexerTracks)
TObject(),
fVert(),
fCurrentVertex(0),
+fFieldkG(-999.),
fConstraint(kFALSE),
fOnlyFitter(kFALSE),
fMinTracks(1),
fMaxd0z0(0.5),
fITSin(kTRUE),
fITSrefit(kTRUE),
+fnSigmaForUi00(1.5),
fDebug(0)
{
//
SetNSigmad0();
SetMaxd0z0();
}
-//-----------------------------------------------------------------------------
-AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
+//----------------------------------------------------------------------------
+AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
TObject(),
fVert(),
fCurrentVertex(0),
+fFieldkG(-999.),
fConstraint(kFALSE),
fOnlyFitter(kFALSE),
fMinTracks(1),
fMaxd0z0(0.5),
fITSin(kTRUE),
fITSrefit(kTRUE),
+fnSigmaForUi00(1.5),
fDebug(0)
{
//
// Standard constructor
//
- SetVtxStart(xStart,yStart);
+ SetVtxStart();
SetVtxStartSigma();
SetMinTracks();
SetMinITSClusters();
SetNSigmad0();
- SetMaxd0z0();
+ SetMaxd0z0();
+ SetFieldkG(fieldkG);
}
//-----------------------------------------------------------------------------
AliVertexerTracks::~AliVertexerTracks()
// The objects pointed by the following pointer are not owned
// by this class and are not deleted
fCurrentVertex = 0;
- delete[] fTrksToSkip;
+ if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
}
//----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
{
//
// Primary vertex for current ESD event
//
fCurrentVertex = 0;
+ // accept 1-track case only if constraint is available
+ if(!fConstraint && fMinTracks==1) fMinTracks=2;
+
// read tracks from ESD
Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
- if (nTrksTot<=0) {
+ if(nTrksTot<=0) {
if(fDebug) printf("TooFewTracks\n");
TooFewTracks(esdEvent);
return fCurrentVertex;
} // end loop on the two iterations
- // take true pos from SPD vertex in ESD and write it in tracks' vertex
- Double_t tp[3];
- esdEvent->GetVertex()->GetTruePos(tp);
- fCurrentVertex->SetTruePos(tp);
-
if(fConstraint) {
if(fOnlyFitter) {
fCurrentVertex->SetTitle("VertexerTracksWithConstraintOnlyFitter");
fTrkArray.Delete();
- if(fTrksToSkip) delete [] fTrksToSkip;
+ if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
if(fDebug) fCurrentVertex->PrintStatus();
if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArray.GetEntries());
AliESDtrack *track1;
track1 = (AliESDtrack*)fTrkArray.At(0);
- Double_t field=GetField();
+ Double_t field=GetFieldkG();
Double_t alpha=track1->GetAlpha();
Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
Double_t pos[3],dir[3];
Double_t initPos[3];
initPos[2] = 0.;
for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
- Double_t field=GetField();
+ Double_t field=GetFieldkG();
Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
Double_t normdistx,normdisty;
Float_t d0z0[2],covd0z0[3];
Double_t sigma;
- Double_t field=GetField();
+ Double_t field=GetFieldkG();
AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
// Removes tracks in trksTree from fit of inVtx
//
+ if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
+ printf("ERROR: primary vertex has no constraint: cannot remove tracks\n");
+ return 0x0;
+ }
if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
printf("WARNING: result of tracks' removal will be only approximately correct\n");
}
Double_t alpha = track->GetAlpha();
Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
- track->AliExternalTrackParam::PropagateTo(xl,AliTracker::GetBz());
+ track->AliExternalTrackParam::PropagateTo(xl,GetFieldkG());
// vector of track global coordinates
TMatrixD ri(3,1);
// covariance matrix of ri
// store data in the vertex object
AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
outVtx->SetTitle(inVtx->GetTitle());
- Double_t tp[3];
- inVtx->GetTruePos(tp);
- outVtx->SetTruePos(tp);
UShort_t *inindices = inVtx->GetIndices();
UShort_t *outindices = new UShort_t[outVtx->GetNContributors()];
Int_t j=0;
}
//---------------------------------------------------------------------------
void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
+{
+ AliESDtrack *track1;
+ Double_t field=GetFieldkG();
+ const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
+ TClonesArray *linarray = new TClonesArray("AliStrLine",1000);
+ TClonesArray &lines = *linarray;
+ for(Int_t i=0; i<knacc; i++){
+ 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],dir[3],sigmasq[3];
+ track1->GetXYZAt(mindist,field,pos);
+ track1->GetPxPyPzAt(mindist,field,dir);
+ sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
+ sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
+ sigmasq[2]=track1->GetSigmaZ2();
+ TMatrixD ri(3,1);
+ TMatrixD wWi(3,3);
+ if(!TrackToPoint(track1,ri,wWi)) continue;
+ Double_t wmat[9];
+ Int_t iel=0;
+ for(Int_t ia=0;ia<3;ia++){
+ for(Int_t ib=0;ib<3;ib++){
+ wmat[iel]=wWi(ia,ib);
+ iel++;
+ }
+ }
+ new(lines[i]) AliStrLine(pos,sigmasq,wmat,dir);
+ }
+ fVert=TrackletVertexFinder(linarray,optUseWeights);
+ linarray->Delete();
+ delete linarray;
+}
+//---------------------------------------------------------------------------
+AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
{
// Calculate the point at minimum distance to prepared tracks
- Double_t initPos[3];
- initPos[2] = 0.;
- Double_t sigma=0;
- for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
- const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
- Double_t field=GetField();
+ const Int_t knacc = (Int_t)lines->GetEntriesFast();
+ Double_t initPos[3]={0.,0.,0.};
- AliESDtrack *track1;
Double_t (*vectP0)[3]=new Double_t [knacc][3];
Double_t (*vectP1)[3]=new Double_t [knacc][3];
Double_t sum[3][3];
Double_t dsum[3]={0,0,0};
- for(Int_t i=0;i<3;i++)
- for(Int_t j=0;j<3;j++)sum[i][j]=0;
- for(Int_t i=0; i<knacc; i++){
- 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],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);
+ TMatrixD sumWi(3,3);
+ for(Int_t i=0;i<3;i++){
+ for(Int_t j=0;j<3;j++){
+ sum[i][j]=0;
+ sumWi(i,j)=0.;
+ }
+ }
- Double_t p0[3],cd[3];
+ for(Int_t i=0; i<knacc; i++){
+ AliStrLine* line1 = (AliStrLine*)lines->At(i);
+ Double_t p0[3],cd[3],sigmasq[3];
+ Double_t wmat[9];
line1->GetP0(p0);
line1->GetCd(cd);
+ line1->GetSigma2P0(sigmasq);
+ line1->GetWMatrix(wmat);
+ TMatrixD wWi(3,3);
+ Int_t iel=0;
+ for(Int_t ia=0;ia<3;ia++){
+ for(Int_t ib=0;ib<3;ib++){
+ wWi(ia,ib)=wmat[iel];
+ iel++;
+ }
+ }
+
+ sumWi+=wWi;
+
Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
vectP0[i][0]=p0[0];
vectP0[i][1]=p0[1];
Double_t matr[3][3];
Double_t dknow[3];
if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
- if(optUseWeights==1){
- Double_t sigmasq[3];
- sigmasq[0]=track1->GetSigmaY2();
- sigmasq[1]=track1->GetSigmaY2();
- sigmasq[2]=track1->GetSigmaZ2();
- GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
- }
+ else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
+
for(Int_t iii=0;iii<3;iii++){
dsum[iii]+=dknow[iii];
for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
}
- delete line1;
}
-
+
+ TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
+ Double_t covmatrix[6];
+ covmatrix[0] = invsumWi(0,0);
+ covmatrix[1] = invsumWi(0,1);
+ covmatrix[2] = invsumWi(1,1);
+ covmatrix[3] = invsumWi(0,2);
+ covmatrix[4] = invsumWi(1,2);
+ covmatrix[5] = invsumWi(2,2);
+
Double_t vett[3][3];
Double_t det=GetDeterminant3X3(sum);
+ Double_t sigma=0;
if(det!=0){
for(Int_t zz=0;zz<3;zz++){
sigma=TMath::Sqrt(sigma);
}else{
- Warning("StrLinVertexFinderMinDist","Finder did not succed");
sigma=999;
}
+ AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
+ theVert.SetDispersion(sigma);
delete vectP0;
delete vectP1;
- fVert.SetXYZ(initPos);
- fVert.SetDispersion(sigma);
- fVert.SetNContributors(knacc);
+ return theVert;
}
//---------------------------------------------------------------------------
Bool_t AliVertexerTracks::TrackToPoint(AliESDtrack *t,
- TMatrixD &ri,TMatrixD &wWi) const
+ TMatrixD &ri,TMatrixD &wWi,
+ Bool_t uUi3by3) const
{
//
// Extract from the AliESDtrack the global coordinates ri and covariance matrix
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);
- Double_t cc[15];
- t->GetExternalCovariance(cc);
- uUi(0,0) = cc[0];
- uUi(0,1) = cc[1];
- uUi(1,0) = cc[1];
- uUi(1,1) = cc[2];
- if(uUi.Determinant() <= 0.) return kFALSE;
- TMatrixD uUiInv(TMatrixD::kInverted,uUi);
+
+
+ if(!uUi3by3) {
+ // 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
+ Double_t cc[15];
+ t->GetExternalCovariance(cc);
+ TMatrixD uUi(2,2);
+ uUi(0,0) = cc[0];
+ uUi(0,1) = cc[1];
+ uUi(1,0) = cc[1];
+ uUi(1,1) = cc[2];
+ //printf(" Ui :\n");
+ //printf(" %f %f\n",uUi(0,0),uUi(0,1));
+ //printf(" %f %f\n",uUi(1,0),uUi(1,1));
+
+ if(uUi.Determinant() <= 0.) return kFALSE;
+ TMatrixD uUiInv(TMatrixD::kInverted,uUi);
+
+ // weights matrix: wWi = qQiT * uUiInv * qQi
+ TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
+ TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
+ wWi = m;
+ } else {
+ if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
+ // matrix to go from global (x,y,z) to local (x,y,z);
+ TMatrixD qQi(3,3);
+ qQi(0,0) = cosRot;
+ qQi(0,1) = sinRot;
+ qQi(0,2) = 0.;
+ qQi(1,0) = -sinRot;
+ qQi(1,1) = cosRot;
+ qQi(1,2) = 0.;
+ qQi(2,0) = 0.;
+ qQi(2,1) = 0.;
+ qQi(2,2) = 1.;
+
+ // covariance of fVert along the track
+ Double_t p[3],pt,ptot;
+ t->GetPxPyPz(p);
+ pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
+ ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
+ Double_t cphi = p[0]/pt; //cos(phi)=px/pt
+ Double_t sphi = p[1]/pt; //sin(phi)=py/pt
+ Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
+ Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
+ Double_t covfVert[6];
+ fVert.GetCovMatrix(covfVert);
+ Double_t covfVertalongt =
+ covfVert[0]*cphi*cphi*clambda*clambda
+ +covfVert[1]*2.*cphi*sphi*clambda*clambda
+ +covfVert[3]*2.*cphi*clambda*slambda
+ +covfVert[2]*sphi*sphi*clambda*clambda
+ +covfVert[4]*2.*sphi*clambda*slambda
+ +covfVert[5]*slambda*slambda;
+ Double_t cc[15];
+ t->GetExternalCovariance(cc);
+ // covariance matrix of local (x,y,z) - inverted
+ TMatrixD uUi(3,3);
+ uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
+ if(fDebug) printf("=====> sqrtUi00 cm %f\n",TMath::Sqrt(uUi(0,0)));
+ uUi(0,1) = 0.;
+ uUi(0,2) = 0.;
+ uUi(1,0) = 0.;
+ uUi(1,1) = cc[0];
+ uUi(1,2) = cc[1];
+ uUi(2,0) = 0.;
+ uUi(2,1) = cc[1];
+ uUi(2,2) = cc[2];
+ //printf(" Ui :\n");
+ //printf(" %f %f\n",uUi(0,0),uUi(0,1));
+ //printf(" %f %f\n",uUi(1,0),uUi(1,1));
+
+ if(uUi.Determinant() <= 0.) return kFALSE;
+ TMatrixD uUiInv(TMatrixD::kInverted,uUi);
- // weights matrix: wWi = qQiT * uUiInv * qQi
- TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
- TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
- wWi = m;
+ // weights matrix: wWi = qQiT * uUiInv * qQi
+ TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
+ TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
+ wWi = m;
+ }
+
return kTRUE;
}
//---------------------------------------------------------------------------
-void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent)
+void AliVertexerTracks::TooFewTracks(const AliESDEvent* esdEvent)
{
//
// When the number of tracks is < fMinTracks
fCurrentVertex = new AliESDVertex(pos,err);
fCurrentVertex->SetNContributors(ncontr);
- Double_t tp[3];
- esdEvent->GetVertex()->GetTruePos(tp);
- fCurrentVertex->SetTruePos(tp);
if(fConstraint) {
fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
} else {
AliESDtrack *track2;
Double_t pos[3],dir[3];
Double_t alpha,mindist;
- Double_t field=GetField();
+ Double_t field=GetFieldkG();
for(Int_t i=0; i<nacc; i++){
track1 = (AliESDtrack*)fTrkArray.At(i);
{
//
// The optimal estimate of the vertex position is given by a "weighted
-// average of tracks positions"
+// average of tracks positions".
// Original method: V. Karimaki, CMS Note 97/0051
//
Double_t initPos[3];
- fVert.GetXYZ(initPos);
- Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
- if(arrEntries==1) useConstraint=kTRUE;
+ if(!fOnlyFitter) {
+ fVert.GetXYZ(initPos);
+ } else {
+ initPos[0]=fNominalPos[0];
+ initPos[1]=fNominalPos[1];
+ initPos[2]=fNominalPos[2];
+ }
+
+ Int_t nTrks = (Int_t)fTrkArray.GetEntries();
+ if(nTrks==1) useConstraint=kTRUE;
if(fDebug) {
- printf(" VertexFitter(): start\n");
- printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
+ printf("--- VertexFitter(): start\n");
+ printf(" Number of tracks in array: %d\n",nTrks);
printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
- printf("Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
+ printf(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
if(useConstraint) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2]));
}
+ // special treatment for few-tracks fits (e.g. secondary vertices)
+ Bool_t uUi3by3 = kFALSE; if(nTrks<5 && !useConstraint) uUi3by3 = kTRUE;
+
Int_t i,j,k,step=0;
TMatrixD rv(3,1);
TMatrixD vV(3,3);
chi2 = 0.;
nUsedTrks = 0;
+ if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
+
TMatrixD sumWiri(3,1);
TMatrixD sumWi(3,3);
for(i=0; i<3; i++) {
chi2 += chi2b;
}
-
// loop on tracks
- for(k=0; k<arrEntries; k++) {
+ for(k=0; k<nTrks; 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);
// to vtxSeed (from finder)
- t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz());
+ t->AliExternalTrackParam::PropagateTo(xlStart,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(t,ri,wWi)) continue;
-
+ if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
// track chi2
continue;
}
- // inverted of weights matrix
- TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
- vV = invsumWi;
-
- // position of primary vertex
- rv.Mult(vV,sumWiri);
-
+ if(step==0) {
+ // 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) {
if(fDebug) printf("TooFewTracks\n");
fCurrentVertex = new AliESDVertex(0.,0.,-1);
covmatrix[4] = vV(1,2);
covmatrix[5] = vV(2,2);
+ // for correct chi2/ndf, count constraint as additional "track"
+ if(fConstraint) nUsedTrks++;
// 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));
+ printf(" Vertex after fit:\n");
fCurrentVertex->PrintStatus();
+ printf("--- VertexFitter(): finish\n");
}
return;
// vertex fitter
if(optUseFitter){
+ //SetVtxStart(&fVert);
VertexFitter(fConstraint);
if(fDebug) printf(" Vertex fit completed\n");
}else{
Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
- Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
+ Double_t covmatrix[6];
+ fVert.GetCovMatrix(covmatrix);
Double_t chi2=99999.;
Int_t nUsedTrks=fVert.GetNContributors();
fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
indices[jj] = (UShort_t)eta->GetID();
if(optPropagate&&optUseFitter){
if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
- eta->RelateToVertex(fCurrentVertex,GetField(),100.);
+ eta->RelateToVertex(fCurrentVertex,GetFieldkG(),100.);
if(fDebug) printf("Track %d propagated to found vertex\n",jj);
}else{
AliWarning("Found vertex outside beam pipe!");
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff