//---- Root headers --------
#include <TSystem.h>
+#include <TClonesArray.h>
#include <TDirectory.h>
#include <TFile.h>
-#include <TTree.h>
-#include <TMatrixD.h>
//---- AliRoot headers -----
#include "AliStrLine.h"
#include "AliExternalTrackParam.h"
-#include "AliESDEvent.h"
+#include "AliNeutralTrackParam.h"
+#include "AliVEvent.h"
+#include "AliVTrack.h"
#include "AliESDtrack.h"
#include "AliVertexerTracks.h"
fMinDetFitter(100.),
fMaxTgl(1000.),
fITSrefit(kTRUE),
+fITSpureSA(kFALSE),
fFiducialR(3.),
fFiducialZ(30.),
fnSigmaForUi00(1.5),
-fDebug(0),
-fAlgo(1)
+fAlgo(1),
+fAlgoIter0(4)
{
//
// Default constructor
fMinDetFitter(100.),
fMaxTgl(1000.),
fITSrefit(kTRUE),
+fITSpureSA(kFALSE),
fFiducialR(3.),
fFiducialZ(30.),
fnSigmaForUi00(1.5),
-fDebug(0),
-fAlgo(1)
+fAlgo(1),
+fAlgoIter0(4)
{
//
// Standard constructor
//
SetVtxStart();
SetVtxStartSigma();
- SetTPCMode();
}
//-----------------------------------------------------------------------------
AliVertexerTracks::~AliVertexerTracks()
if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
}
//----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
{
//
-// Primary vertex for current ESD event
+// Primary vertex for current ESD or AOD event
// (Two iterations:
// 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
// + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
//
fCurrentVertex = 0;
+ TString evtype = vEvent->IsA()->GetName();
+ Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
+
+ if(inputAOD && fMode==1) {
+ printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n");
+ TooFewTracks();
+ return fCurrentVertex;
+ }
+
// accept 1-track case only if constraint is available
if(!fConstraint && fMinTracks==1) fMinTracks=2;
- // read tracks from ESD
- Int_t nTrks = (Int_t)esdEvent->GetNumberOfTracks();
+ // read tracks from AlivEvent
+ Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
if(nTrks<fMinTracks) {
TooFewTracks();
return fCurrentVertex;
Int_t nTrksOrig=0;
AliExternalTrackParam *t=0;
- // loop on ESD tracks
+ // loop on tracks
for(Int_t i=0; i<nTrks; i++) {
- AliESDtrack *esdt = esdEvent->GetTrack(i);
+ AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
// check tracks to skip
Bool_t skipThis = kFALSE;
for(Int_t j=0; j<fNTrksToSkip; j++) {
- if(esdt->GetID()==fTrksToSkip[j]) {
- if(fDebug) printf("skipping track: %d\n",i);
+ if(track->GetID()==fTrksToSkip[j]) {
+ AliDebug(1,Form("skipping track: %d",i));
skipThis = kTRUE;
}
}
if(skipThis) continue;
- // check number of clusters in ITS or TPC
- if(esdt->GetNcls(fMode) < fMinClusters) continue;
-
- if(fMode==0) { // ITS mode
- if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
- Double_t x,p[5],cov[15];
- esdt->GetExternalParameters(x,p);
- esdt->GetExternalCovariance(cov);
- t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
- } else if(fMode==1) { // TPC mode
- t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
- if(!t) continue;
- Double_t radius = 2.8; //something less than the beam pipe radius
- if(!PropagateTrackTo(t,radius)) continue;
+ // skip pure ITS SA tracks (default)
+ if(!fITSpureSA && (track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
+ // or use only pure ITS SA tracks
+ if(fITSpureSA && !(track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
+
+ // kITSrefit
+ if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue;
+
+ if(!inputAOD) { // ESD
+ AliESDtrack* esdt = (AliESDtrack*)track;
+ if(esdt->GetNcls(fMode) < fMinClusters) continue;
+ if(fMode==0) { // ITS mode
+ Double_t x,p[5],cov[15];
+ esdt->GetExternalParameters(x,p);
+ esdt->GetExternalCovariance(cov);
+ t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
+ } else if(fMode==1) { // TPC mode
+ t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
+ if(!t) continue;
+ Double_t radius = 2.8; //something less than the beam pipe radius
+ if(!PropagateTrackTo(t,radius)) continue;
+ }
+ } else { // AOD (only ITS mode)
+ Int_t ncls0=0;
+ for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
+ if(ncls0 < fMinClusters) continue;
+ t = new AliExternalTrackParam(track);
}
trkArrayOrig.AddLast(t);
- idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
+ idOrig[nTrksOrig]=(UShort_t)track->GetID();
nTrksOrig++;
- } // end loop on ESD tracks
+ } // end loop on tracks
// call method that will reconstruct the vertex
FindPrimaryVertex(&trkArrayOrig,idOrig);
olddir->cd();
}
+ // set vertex ID for tracks used in the fit
+ // (only for ESD)
+ if(!inputAOD) {
+ Int_t nIndices = fCurrentVertex->GetNIndices();
+ UShort_t *indices = fCurrentVertex->GetIndices();
+ for(Int_t ind=0; ind<nIndices; ind++) {
+ AliESDtrack *esdt = (AliESDtrack*)vEvent->GetTrack(indices[ind]);
+ esdt->SetVertexID(-1);
+ }
+ }
+
return fCurrentVertex;
}
//----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig,
UShort_t *idOrig)
{
//
// read tracks from array
Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
- if(fDebug) printf("Initial number of tracks: %d\n",nTrksOrig);
+ AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
if(nTrksOrig<fMinTracks) {
- if(fDebug) printf("TooFewTracks\n");
+ AliDebug(1,"TooFewTracks");
TooFewTracks();
return fCurrentVertex;
}
if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
Double_t cutsave = fDCAcut;
fDCAcut = fDCAcutIter0;
- VertexFinder(1); // using weights, cutting dca < fDCAcutIter0
+ // vertex finder
+ switch (fAlgoIter0) {
+ 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;
+ }
fDCAcut = cutsave;
if(fVert.GetNContributors()>0) {
fVert.GetXYZ(fNominalPos);
fNominalPos[0] = fVert.GetXv();
fNominalPos[1] = fVert.GetYv();
fNominalPos[2] = fVert.GetZv();
- if(fDebug) printf("No mean vertex: VertexFinder gives (%f, %f, %f)\n",fNominalPos[0],fNominalPos[1],fNominalPos[2]);
+ AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
} else {
fNominalPos[0] = 0.;
fNominalPos[1] = 0.;
fNominalPos[2] = 0.;
- if(fDebug) printf("No mean vertex and VertexFinder failed\n");
+ AliDebug(1,"No mean vertex and VertexFinder failed");
}
}
if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
fIdSel = new UShort_t[nTrksOrig];
Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
- if(fDebug) printf("N tracks selected in iteration %d: %d\n",iter,nTrksSel);
+ AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
if(nTrksSel < fMinTracks) {
TooFewTracks();
return fCurrentVertex;
// vertex finder
if(!fOnlyFitter) {
if(nTrksSel==1) {
- if(fDebug) printf("Just one track\n");
+ AliDebug(1,"Just one track");
OneTrackVertFinder();
} else {
switch (fAlgo) {
case 3: HelixVertexFinder(); break;
case 4: VertexFinder(1); break;
case 5: VertexFinder(0); break;
- default: printf("Wrong algorithm\n"); break;
+ default: printf("Wrong algorithm"); break;
}
}
- if(fDebug) printf(" Vertex finding completed\n");
+ AliDebug(1," Vertex finding completed");
}
// vertex fitter
indices[jj] = fIdSel[jj];
fCurrentVertex->SetIndices(nIndices,indices);
}
- delete [] indices; indices=NULL;
+ if (indices) {delete [] indices; indices=NULL;}
//
// set vertex title
fCurrentVertex->SetTitle(title.Data());
//
- if(fDebug) {
- fCurrentVertex->PrintStatus();
- fCurrentVertex->PrintIndices();
- }
+ AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
// clean up
delete [] fIdSel; fIdSel=NULL;
return det;
}
//-------------------------------------------------------------------------
-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
+void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const 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(const Double_t *p0,const Double_t *p1,const Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
{
//
Double_t x12=p1[0]-p0[0];
}
//--------------------------------------------------------------------------
-Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
+Double_t AliVertexerTracks::GetStrLinMinDist(const Double_t *p0,const Double_t *p1,const Double_t *x0)
{
//
Double_t x12=p0[0]-p1[0];
void AliVertexerTracks::OneTrackVertFinder()
{
// find vertex for events with 1 track, using DCA to nominal beam axis
- if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries());
+ AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
AliExternalTrackParam *track1;
track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
Double_t alpha=track1->GetAlpha();
fVert.SetNContributors(ncombi);
}
//----------------------------------------------------------------------------
-Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
- UShort_t *idOrig,
+Int_t AliVertexerTracks::PrepareTracks(const TObjArray &trkArrayOrig,
+ const UShort_t *idOrig,
Int_t optImpParCut)
{
//
// Propagate tracks to initial vertex position and store them in a TObjArray
//
- if(fDebug) printf(" PrepareTracks()\n");
+ AliDebug(1," PrepareTracks()");
Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
Int_t nTrksSel = 0;
// loop on tracks
for(Int_t i=0; i<nTrksOrig; i++) {
- track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+ AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
+ if(trackOrig->Charge()!=0) { // normal tracks
+ track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+ } else { // neutral tracks (from a V0)
+ track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
+ }
+
// tgl cut
if(TMath::Abs(track->GetTgl())>fMaxTgl) {
- if(fDebug) printf(" rejecting track with tgl = %f\n",track->GetTgl());
+ AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
delete track; continue;
}
- Bool_t propagateOK = kFALSE;
+ Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE;
// propagate track to vertex
if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
fCurrentVertex->GetSigmaXYZ(sigmaCurr);
normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
+ AliDebug(1,Form("normdistx %f %f %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
+ AliDebug(1,Form("normdisty %f %f %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2])));
+ AliDebug(1,Form("sigmaCurr %f %f %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2])));
if(normdistx < 3. && normdisty < 3. &&
(sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
} else {
propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
+ if(fConstraint) cutond0z0=kFALSE;
}
}
if(!propagateOK) {
- if(fDebug) printf(" rejected\n");
+ AliDebug(1," rejected");
delete track; continue;
}
maxd0rphi = fNSigma*sigmad0;
if(optImpParCut==1) maxd0rphi *= 5.;
maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
- //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
- //maxd0z0 = 10.*fNSigma*sigmad0z0;
+ //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]);
- if(fDebug) printf("trk %d; id %d; |d0| = %f; d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0);
+ AliDebug(1,Form("trk %d; id %d; |d0| = %f; d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0));
//---- track selection based on impact parameters ----//
// always reject tracks outside fiducial volume
if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
- if(fDebug) printf(" rejected\n");
+ AliDebug(1," rejected");
delete track; continue;
}
// during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
- if(fDebug) printf(" rejected\n");
+ AliDebug(1," rejected");
delete track; continue;
}
// if fConstraint=kTRUE, during iteration 2,
// select tracks with d0oplusz0 < fMaxd0z0
if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
- ( fConstraint && optImpParCut==2)) {
+ ( fConstraint && optImpParCut==2 && cutond0z0)) {
if(nTrksOrig>=3 &&
TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
- if(fDebug) printf(" rejected\n");
+ AliDebug(1," rejected");
delete track; continue;
}
}
//
Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
sa=TMath::Sin(alphan-track->GetAlpha());
- Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
+ 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;
}
//---------------------------------------------------------------------------
AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
- TObjArray *trkArray,
+ const TObjArray *trkArray,
UShort_t *id,
- Float_t *diamondxy)
+ 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\n");
+ 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\n");
+ printf("WARNING: result of tracks' removal will be only approximately correct");
TMatrixD rv(3,1);
rv(0,0) = inVtx->GetXv();
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\n",id[i]);
+ printf("track %d was not used in vertex fit",id[i]);
continue;
}
Double_t alpha = track->GetAlpha();
sumWi -= wWi;
sumWiri -= wWiri;
- // track chi2
+ // track contribution to chi2
TMatrixD deltar = rv; deltar -= ri;
TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
nUsedTrks--;
if(nUsedTrks<2) {
- printf("Trying to remove too many tracks!\n");
+ printf("Trying to remove too many tracks!");
return 0x0;
}
}
cov[5] = vVnew(2,2);
// store data in the vertex object
- AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
+ 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;
}
}
outVtx->SetIndices(nIndices,outindices);
- delete [] outindices;
+ if (outindices) delete [] outindices;
- if(fDebug) {
- printf("Vertex before removing tracks:\n");
+ /*
+ printf("Vertex before removing tracks:");
inVtx->PrintStatus();
inVtx->PrintIndices();
- printf("Vertex after removing tracks:\n");
+ 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;
}
//---------------------------------------------------------------------------
SetDCAcut(cuts[0]);
SetDCAcutIter0(cuts[1]);
SetMaxd0z0(cuts[2]);
- SetMinClusters((Int_t)(cuts[3]));
+ 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;
}
Double_t mindetfitter,
Double_t maxtgl,
Double_t fidR,
- Double_t fidZ)
+ Double_t fidZ,
+ Int_t finderAlgo,
+ Int_t finderAlgoIter0)
{
//
// Cut values for ITS mode
//
fMode = 0;
- SetITSrefitRequired();
+ if(minCls>0) {
+ SetITSrefitRequired();
+ } else {
+ SetITSrefitNotRequired();
+ }
SetDCAcut(dcacut);
SetDCAcutIter0(dcacutIter0);
SetMaxd0z0(maxd0z0);
- SetMinClusters(minCls);
+ SetMinClusters(TMath::Abs(minCls));
SetMinTracks(mintrks);
SetNSigmad0(nsigma);
SetMinDetFitter(mindetfitter);
SetMaxTgl(maxtgl);
SetFiducialRZ(fidR,fidZ);
+ fAlgo=finderAlgo;
+ fAlgoIter0=finderAlgoIter0;
return;
}
Double_t mindetfitter,
Double_t maxtgl,
Double_t fidR,
- Double_t fidZ)
+ Double_t fidZ,
+ Int_t finderAlgo,
+ Int_t finderAlgoIter0)
{
//
// Cut values for TPC mode
SetMinDetFitter(mindetfitter);
SetMaxTgl(maxtgl);
SetFiducialRZ(fidR,fidZ);
+ fAlgo=finderAlgo;
+ fAlgoIter0=finderAlgoIter0;
return;
}
//---------------------------------------------------------------------------
-void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
+void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped)
{
//
// Mark the tracks not to be used in the vertex reconstruction.
{
AliExternalTrackParam *track1;
const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
- static TClonesArray linarray("AliStrLine",knacc);
+ AliStrLine **linarray = new AliStrLine* [knacc];
for(Int_t i=0; i<knacc; i++){
track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
Double_t alpha=track1->GetAlpha();
sigmasq[2]=track1->GetSigmaZ2();
TMatrixD ri(3,1);
TMatrixD wWi(3,3);
- if(!TrackToPoint(track1,ri,wWi)) continue;
+ if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
Double_t wmat[9];
Int_t iel=0;
for(Int_t ia=0;ia<3;ia++){
iel++;
}
}
- new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
+ linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
}
- fVert=TrackletVertexFinder(&linarray,optUseWeights);
- linarray.Clear("C");
+ fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
+ for(Int_t i=0; i<knacc; i++) delete linarray[i];
+ delete [] linarray;
}
//---------------------------------------------------------------------------
-AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
+AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
{
- // Calculate the point at minimum distance to prepared tracks
-
+ // Calculate the point at minimum distance to prepared tracks (TClonesArray)
const Int_t knacc = (Int_t)lines->GetEntriesFast();
+ AliStrLine** lines2 = new AliStrLine* [knacc];
+ for(Int_t i=0; i<knacc; i++){
+ lines2[i]= (AliStrLine*)lines->At(i);
+ }
+ AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
+ delete [] lines2;
+ return vert;
+}
+
+//---------------------------------------------------------------------------
+AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
+{
+ // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
+
Double_t initPos[3]={0.,0.,0.};
Double_t (*vectP0)[3]=new Double_t [knacc][3];
}
for(Int_t i=0; i<knacc; i++){
- AliStrLine* line1 = (AliStrLine*)lines->At(i);
+ AliStrLine *line1 = lines[i];
Double_t p0[3],cd[3],sigmasq[3];
Double_t wmat[9];
+ if(!line1) printf("ERROR %d %d\n",i,knacc);
line1->GetP0(p0);
line1->GetCd(cd);
line1->GetSigma2P0(sigmasq);
Double_t det=GetDeterminant3X3(sum);
Double_t sigma=0;
- if(det!=0){
+ if(TMath::Abs(det) > kAlmost0){
for(Int_t zz=0;zz<3;zz++){
for(Int_t ww=0;ww<3;ww++){
for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
}
AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
theVert.SetDispersion(sigma);
- delete vectP0;
- delete vectP1;
+ delete [] vectP0;
+ delete [] vectP1;
return theVert;
}
+
//---------------------------------------------------------------------------
Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
TMatrixD &ri,TMatrixD &wWi,
uUi(0,1) = t->GetSigmaZY();
uUi(1,0) = t->GetSigmaZY();
uUi(1,1) = t->GetSigmaZ2();
- //printf(" Ui :\n");
- //printf(" %f %f\n",uUi(0,0),uUi(0,1));
- //printf(" %f %f\n",uUi(1,0),uUi(1,1));
+ //printf(" Ui :");
+ //printf(" %f %f",uUi(0,0),uUi(0,1));
+ //printf(" %f %f",uUi(1,0),uUi(1,1));
if(uUi.Determinant() <= 0.) return kFALSE;
TMatrixD uUiInv(TMatrixD::kInverted,uUi);
// 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)));
+ AliDebug(1,Form("=====> sqrtUi00 cm %f",TMath::Sqrt(uUi(0,0))));
uUi(0,1) = 0.;
uUi(0,2) = 0.;
uUi(1,0) = 0.;
// When the number of tracks is < fMinTracks,
// deal with vertices not found and prepare to exit
//
- if(fDebug) printf("TooFewTracks\n");
+ AliDebug(1,"TooFewTracks");
Double_t pos[3],err[3];
pos[0] = fNominalPos[0];
pos[2] = fNominalPos[2];
err[2] = TMath::Sqrt(fNominalCov[5]);
Int_t ncontr = (err[0]>1. ? -1 : -3);
- fCurrentVertex = 0;
+ if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
fCurrentVertex = new AliESDVertex(pos,err);
fCurrentVertex->SetNContributors(ncontr);
Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
if(nTrksSel==1) useConstraint=kTRUE;
- if(fDebug) {
- printf("--- VertexFitter(): start\n");
- printf(" Number of tracks in array: %d\n",nTrksSel);
- 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(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]));
- }
+ AliDebug(1,"--- VertexFitter(): start");
+ AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
+ AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
+ AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
+ if(useConstraint) AliDebug(1,Form(" 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(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
// 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);
+ AliDebug(1,Form(" step = %d\n",step));
chi2 = 0.;
nTrksUsed = 0;
Double_t determinant = sumWi.Determinant();
if(determinant < fMinDetFitter) {
- printf("det(V) = %f (<%f)\n",determinant,fMinDetFitter);
+ AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
failed=1;
continue;
}
// for correct chi2/ndf, count constraint as additional "track"
if(fConstraint) nTrksUsed++;
// store data in the vertex object
+ if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
- if(fDebug) {
- printf(" Vertex after fit:\n");
- fCurrentVertex->PrintStatus();
- printf("--- VertexFitter(): finish\n");
- }
+ AliDebug(1," Vertex after fit:");
+ AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
+ AliDebug(1,"--- VertexFitter(): finish\n");
+
return;
}
//----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
+AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
UShort_t *id,
Bool_t optUseFitter,
- Bool_t optPropagate)
+ Bool_t optPropagate,
+ Bool_t optUseDiamondConstraint)
{
//
// Return vertex from tracks (AliExternalTrackParam) in array
//
-
- SetConstraintOff();
+ fCurrentVertex = 0;
+ // set optUseDiamondConstraint=TRUE only if you are reconstructing the
+ // primary vertex!
+ if(optUseDiamondConstraint) {
+ SetConstraintOn();
+ } else {
+ SetConstraintOff();
+ }
// get tracks and propagate them to initial vertex position
fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
- if(nTrksSel < TMath::Max(2,fMinTracks)) {
+ if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
+ (optUseDiamondConstraint && nTrksSel<1)) {
TooFewTracks();
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(nTrksSel==1) {
+ AliDebug(1,"Just one track");
+ OneTrackVertFinder();
+ } 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;
+ }
}
- if(fDebug) printf(" Vertex finding completed\n");
+ AliDebug(1," Vertex finding completed\n");
// vertex fitter
if(optUseFitter) {
Double_t d0z0[2],covd0z0[3];
AliExternalTrackParam *t = 0;
if(fCurrentVertex->GetNContributors()>0) {
- indices = new UShort_t[fCurrentVertex->GetNContributors()];
+ indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
indices[jj] = fIdSel[jj];
t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
if(optPropagate && optUseFitter) {
if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
- if(fDebug) printf("Track %d propagated to found vertex\n",jj);
+ AliDebug(1,Form("Track %d propagated to found vertex",jj));
} else {
AliWarning("Found vertex outside beam pipe!");
}
}
// clean up
- delete [] indices; indices=NULL;
+ if (indices) {delete [] indices; indices=NULL;}
delete [] fIdSel; fIdSel=NULL;
fTrkArraySel.Delete();
}
//----------------------------------------------------------------------------
AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
- Bool_t optUseFitter,
- Bool_t optPropagate)
+ Bool_t optUseFitter,
+ Bool_t optPropagate,
+ Bool_t optUseDiamondConstraint)
+
{
//
// Return vertex from array of ESD tracks
id[i] = (UShort_t)esdt->GetID();
}
- VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
+ VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
delete [] id; id=NULL;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff