#include <TTree.h>
#include <TMatrixD.h>
//---- AliRoot headers -----
+#include "AliLog.h"
#include "AliStrLine.h"
#include "AliExternalTrackParam.h"
+#include "AliNeutralTrackParam.h"
+#include "AliVEvent.h"
+#include "AliVTrack.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliVertexerTracks.h"
TObject(),
fVert(),
fCurrentVertex(0),
+fMode(0),
fFieldkG(-999.),
-fConstraint(kFALSE),
-fOnlyFitter(kFALSE),
-fMinTracks(1),
-fMinITSClusters(5),
fTrkArraySel(),
fIdSel(0),
fTrksToSkip(0),
fNTrksToSkip(0),
-fDCAcut(0),
-fAlgo(1),
-fNSigma(3),
+fConstraint(kFALSE),
+fOnlyFitter(kFALSE),
+fMinTracks(1),
+fMinClusters(5),
+fDCAcut(0.1),
+fDCAcutIter0(0.1),
+fNSigma(3.),
fMaxd0z0(0.5),
+fMinDetFitter(100.),
+fMaxTgl(1000.),
fITSrefit(kTRUE),
+fFiducialR(3.),
+fFiducialZ(30.),
fnSigmaForUi00(1.5),
-fDebug(0)
+fAlgo(1),
+fAlgoIter0(4)
{
//
// Default constructor
//
SetVtxStart();
SetVtxStartSigma();
- SetMinTracks();
- SetMinITSClusters();
- SetNSigmad0();
- SetMaxd0z0();
}
//----------------------------------------------------------------------------
AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
TObject(),
fVert(),
fCurrentVertex(0),
-fFieldkG(-999.),
-fConstraint(kFALSE),
-fOnlyFitter(kFALSE),
-fMinTracks(1),
-fMinITSClusters(5),
+fMode(0),
+fFieldkG(fieldkG),
fTrkArraySel(),
fIdSel(0),
fTrksToSkip(0),
fNTrksToSkip(0),
-fDCAcut(0),
-fAlgo(1),
-fNSigma(3),
+fConstraint(kFALSE),
+fOnlyFitter(kFALSE),
+fMinTracks(1),
+fMinClusters(5),
+fDCAcut(0.1),
+fDCAcutIter0(0.1),
+fNSigma(3.),
fMaxd0z0(0.5),
+fMinDetFitter(100.),
+fMaxTgl(1000.),
fITSrefit(kTRUE),
+fFiducialR(3.),
+fFiducialZ(30.),
fnSigmaForUi00(1.5),
-fDebug(0)
+fAlgo(1),
+fAlgoIter0(4)
{
//
// Standard constructor
//
SetVtxStart();
SetVtxStartSigma();
- SetMinTracks();
- SetMinITSClusters();
- SetNSigmad0();
- SetMaxd0z0();
- SetFieldkG(fieldkG);
}
//-----------------------------------------------------------------------------
AliVertexerTracks::~AliVertexerTracks()
// The objects pointed by the following pointer are not owned
// by this class and are not deleted
fCurrentVertex = 0;
- if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
- if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
+ if(fTrksToSkip) { delete fTrksToSkip; fTrksToSkip=NULL; }
+ if(fIdSel) { delete fIdSel; fIdSel=NULL; }
}
//----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(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();
- if(nTrks<1) {
+ // read tracks from AlivEvent
+ Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
+ if(nTrks<fMinTracks) {
TooFewTracks();
return fCurrentVertex;
}
//
TDirectory * olddir = gDirectory;
- TFile f("VertexerTracks.root","recreate");
+ TFile *f = 0;
+ if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
TObjArray trkArrayOrig(nTrks);
UShort_t *idOrig = new UShort_t[nTrks];
Int_t nTrksOrig=0;
+ AliExternalTrackParam *t=0;
+ // 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;
- if(fITSrefit) {
- if(!(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
- // check number of clusters in ITS
- if(esdt->GetNcls(0)<fMinITSClusters) 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);
}
- Double_t x,p[5],cov[15];
- esdt->GetExternalParameters(x,p);
- esdt->GetExternalCovariance(cov);
- AliExternalTrackParam *t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
trkArrayOrig.AddLast(t);
- idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
+ idOrig[nTrksOrig]=(UShort_t)track->GetID();
nTrksOrig++;
- }
+ } // end loop on tracks
+ // call method that will reconstruct the vertex
FindPrimaryVertex(&trkArrayOrig,idOrig);
- trkArrayOrig.Delete();
- delete [] idOrig; idOrig=NULL;
- f.Close();
- gSystem->Unlink("VertexerTracks.root");
- olddir->cd();
+ if(fMode==0) trkArrayOrig.Delete();
+ delete idOrig; idOrig=NULL;
+
+ if(f) {
+ f->Close(); delete f; f = NULL;
+ gSystem->Unlink("VertexerTracks.root");
+ 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;
}
// accept 1-track case only if constraint is available
if(!fConstraint && fMinTracks==1) fMinTracks=2;
- // read tracks from ESD
+ // read tracks from array
Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
- if(fDebug) printf("Initial number of tracks: %d\n",nTrksOrig);
- if(nTrksOrig<=0) {
- if(fDebug) printf("TooFewTracks\n");
+ AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
+ if(nTrksOrig<fMinTracks) {
+ AliDebug(1,"TooFewTracks");
TooFewTracks();
return fCurrentVertex;
}
// fill fTrkArraySel, for VertexFinder()
fIdSel = new UShort_t[nTrksOrig];
PrepareTracks(*trkArrayOrig,idOrig,0);
- if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
+ if(fIdSel) { delete fIdSel; fIdSel=NULL; }
Double_t cutsave = fDCAcut;
- fDCAcut = (fITSrefit ? 0.1 : 0.5);
- VertexFinder(1); // using weights, cutting dca < fDCAcut
+ fDCAcut = 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");
}
}
// between initVertex and fCurrentVertex)
for(Int_t iter=1; iter<=2; iter++) {
if(fOnlyFitter && iter==1) continue;
- if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
+ 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;
+ delete fIdSel; fIdSel=NULL;
fTrkArraySel.Delete();
- if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
+ if(fTrksToSkip) { delete fTrksToSkip; fTrksToSkip=NULL; }
//
return fCurrentVertex;
Double_t x10=p0[0]-x0[0];
Double_t y10=p0[1]-x0[1];
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);
+ // 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);
+ return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
+ (z10*x12-x10*z12)*(z10*x12-x10*z12)+
+ (x10*y12-y10*x12)*(x10*y12-y10*x12))
+ /(x12*x12+y12*y12+z12*z12);
}
//---------------------------------------------------------------------------
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 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];
- track1->GetXYZAt(mindist,field,pos);
- track1->GetPxPyPzAt(mindist,field,dir);
+ track1->GetXYZAt(mindist,GetFieldkG(),pos);
+ track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
AliStrLine *line1 = new AliStrLine(pos,dir);
Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
Double_t initPos[3];
initPos[2] = 0.;
for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
- Double_t field=GetFieldkG();
Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
for(Int_t j=i+1; j<nacc; j++){
track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
- distCA=track2->PropagateToDCA(track1,field);
+ distCA=track2->PropagateToDCA(track1,GetFieldkG());
if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
x=track1->GetX();
alpha=track1->GetAlpha();
//
// 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;
Double_t maxd0rphi;
- Double_t maxd0z0 = (fITSrefit ? fMaxd0z0 : 10.*fMaxd0z0);
Double_t sigmaCurr[3];
Double_t normdistx,normdisty;
Double_t d0z0[2],covd0z0[3];
Double_t sigmad0;
- Double_t field = GetFieldkG();
AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
AliExternalTrackParam *track=0;
+ // loop on tracks
for(Int_t i=0; i<nTrksOrig; i++) {
- track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
- // only TPC tracks in |eta|<0.9
- if(!fITSrefit && TMath::Abs(track->GetTgl())>1.) {
- if(fDebug) printf(" rejecting track with tgl = %f\n",track->GetTgl());
- delete track;
- continue;
+ 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) {
+ AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
+ delete track; continue;
+ }
+
+ Bool_t propagateOK = kFALSE;
// propagate track to vertex
if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
- track->PropagateToDCA(initVertex,field,100.,d0z0,covd0z0);
+ propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
} else { // optImpParCut==2
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]);
if(normdistx < 3. && normdisty < 3. &&
(sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
- track->PropagateToDCA(fCurrentVertex,field,100.,d0z0,covd0z0);
+ propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
} else {
- track->PropagateToDCA(initVertex,field,100.,d0z0,covd0z0);
+ propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
}
}
+ if(!propagateOK) {
+ AliDebug(1," rejected");
+ delete track; continue;
+ }
+
sigmad0 = TMath::Sqrt(covd0z0[0]);
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;
- 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]),maxd0z0);
+ 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) {
+ AliDebug(1," rejected");
+ delete track; continue;
+ }
+
+ // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
+ if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
+ AliDebug(1," rejected");
+ delete track; continue;
+ }
// if fConstraint=kFALSE, during iterations 1 and 2 ||
// if fConstraint=kTRUE, during iteration 2,
- // select tracks with d0oplusz0 < maxd0z0
+ // select tracks with d0oplusz0 < fMaxd0z0
if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
( fConstraint && optImpParCut==2)) {
if(nTrksOrig>=3 &&
- TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>maxd0z0) {
- if(fDebug) printf(" rejected\n");
+ TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
+ AliDebug(1," rejected");
delete track; continue;
}
}
-
- // select tracks with d0rphi < maxd0rphi
- if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
- if(fDebug) printf(" rejected\n");
- delete track; continue;
- }
-
+ // track passed all selections
fTrkArraySel.AddLast(track);
fIdSel[nTrksSel] = idOrig[i];
nTrksSel++;
- }
+ } // end loop on tracks
delete initVertex;
return nTrksSel;
}
+//----------------------------------------------------------------------------
+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,
TObjArray *trkArray,
//
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();
nUsedTrks--;
if(nUsedTrks<2) {
- printf("Trying to remove too many tracks!\n");
+ printf("Trying to remove too many tracks!");
return 0x0;
}
}
}
}
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;
}
//---------------------------------------------------------------------------
+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;
+}
+//---------------------------------------------------------------------------
void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
{
//
void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
{
AliExternalTrackParam *track1;
- Double_t field=GetFieldkG();
const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
- TClonesArray *linarray = new TClonesArray("AliStrLine",1000);
- TClonesArray &lines = *linarray;
+ AliStrLine **linarray = new AliStrLine* [knacc];
for(Int_t i=0; i<knacc; i++){
track1 = (AliExternalTrackParam*)fTrkArraySel.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);
+ track1->GetXYZAt(mindist,GetFieldkG(),pos);
+ track1->GetPxPyPzAt(mindist,GetFieldkG(),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;
+ 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(lines[i]) AliStrLine(pos,sigmasq,wmat,dir);
+ linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
}
- fVert=TrackletVertexFinder(linarray,optUseWeights);
- linarray->Delete();
- delete linarray;
+ 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)
{
- // 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);
- if (!line1) continue;
+ 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);
sigma+=GetStrLinMinDist(p0,p1,initPos);
}
- sigma=TMath::Sqrt(sigma);
+ if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
}else{
sigma=999;
}
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);
}
if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
- if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
- if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
+ if(fIdSel) {delete fIdSel; fIdSel=NULL;}
+ if(fTrksToSkip) {delete fTrksToSkip; fTrksToSkip=NULL;}
return;
}
AliExternalTrackParam *track2;
Double_t pos[3],dir[3];
Double_t alpha,mindist;
- Double_t field=GetFieldkG();
for(Int_t i=0; i<nacc; i++){
track1 = (AliExternalTrackParam*)fTrkArraySel.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);
+ track1->GetXYZAt(mindist,GetFieldkG(),pos);
+ track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
AliStrLine *line1 = new AliStrLine(pos,dir);
// AliStrLine *line1 = new AliStrLine();
- // track1->ApproximateHelixWithLine(mindist,field,line1);
+ // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
for(Int_t j=i+1; j<nacc; j++){
track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
alpha=track2->GetAlpha();
mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
- track2->GetXYZAt(mindist,field,pos);
- track2->GetPxPyPzAt(mindist,field,dir);
+ track2->GetXYZAt(mindist,GetFieldkG(),pos);
+ track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
AliStrLine *line2 = new AliStrLine(pos,dir);
// AliStrLine *line2 = new AliStrLine();
- // track2->ApproximateHelixWithLine(mindist,field,line2);
+ // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
Double_t distCA=line2->GetDCA(line1);
//printf("%d %d %f\n",i,j,distCA);
if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
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();
- //cerr<<" determinant: "<<determinant<<endl;
- if(determinant < 100.) {
- printf("det(V) = 0\n");
+ if(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;
}
//
// Return vertex from tracks (AliExternalTrackParam) in array
//
-
+ fCurrentVertex = 0;
SetConstraintOff();
// get tracks and propagate them to initial vertex position
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;
- delete [] fIdSel; fIdSel=NULL;
+ if (indices) {delete indices; indices=NULL;}
+ delete fIdSel; fIdSel=NULL;
fTrkArraySel.Delete();
return fCurrentVertex;
VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
- delete [] id; id=NULL;
+ delete id; id=NULL;
return fCurrentVertex;
}
//--------------------------------------------------------------------------
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff