X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliVertexerTracks.cxx;h=fc87a24f85b514875bc11caa695ea625a8c45db1;hb=8d3cf7a997c59e55b59713aad99d56a945c169ef;hp=97a7626595fcda9330624f9f5e52ca3e7198e3e0;hpb=60e55aee9905a1a5e66dce1cd084737911e78b04;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliVertexerTracks.cxx b/STEER/AliVertexerTracks.cxx index 97a7626595f..fc87a24f85b 100644 --- a/STEER/AliVertexerTracks.cxx +++ b/STEER/AliVertexerTracks.cxx @@ -120,11 +120,11 @@ 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 AliVEvent *vEvent) +AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent) { // // Primary vertex for current ESD or AOD event @@ -176,11 +176,13 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent) } if(skipThis) 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 - if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue; Double_t x,p[5],cov[15]; esdt->GetExternalParameters(x,p); esdt->GetExternalCovariance(cov); @@ -206,7 +208,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent) FindPrimaryVertex(&trkArrayOrig,idOrig); if(fMode==0) trkArrayOrig.Delete(); - delete [] idOrig; idOrig=NULL; + delete idOrig; idOrig=NULL; if(f) { f->Close(); delete f; f = NULL; @@ -214,6 +216,17 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent) 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; indGetTrack(indices[ind]); + esdt->SetVertexID(-1); + } + } + return fCurrentVertex; } //---------------------------------------------------------------------------- @@ -248,7 +261,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig, // 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 = fDCAcutIter0; // vertex finder @@ -288,7 +301,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig, // 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); AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel)); @@ -329,7 +342,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig, indices[jj] = fIdSel[jj]; fCurrentVertex->SetIndices(nIndices,indices); } - delete [] indices; indices=NULL; + if (indices) {delete indices; indices=NULL;} // // set vertex title @@ -344,9 +357,9 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig, 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; @@ -580,6 +593,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig, } 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())); @@ -701,7 +715,7 @@ Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track, AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, TObjArray *trkArray, UShort_t *id, - Float_t *diamondxy) + Float_t *diamondxy) const { // // Removes tracks in trksTree from fit of inVtx @@ -758,7 +772,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, 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)+ @@ -795,7 +809,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, 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; @@ -811,7 +825,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, } } outVtx->SetIndices(nIndices,outindices); - delete [] outindices; + if (outindices) delete outindices; /* printf("Vertex before removing tracks:"); @@ -825,6 +839,106 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, return outVtx; } //--------------------------------------------------------------------------- +AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx, + Float_t *diamondxyz, + Float_t *diamondcov) const +{ +// +// Removes diamond constraint from fit of inVtx +// + + if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) { + printf("ERROR: primary vertex has no constraint: cannot remove it\n"); + return 0x0; + } + if(inVtx->GetNContributors()<3) { + printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n"); + return 0x0; + } + + // diamond constraint + TMatrixD vVb(3,3); + vVb(0,0) = diamondcov[0]; + vVb(0,1) = diamondcov[1]; + vVb(0,2) = 0.; + vVb(1,0) = diamondcov[1]; + vVb(1,1) = diamondcov[2]; + vVb(1,2) = 0.; + vVb(2,0) = 0.; + vVb(2,1) = 0.; + vVb(2,2) = diamondcov[5]; + TMatrixD vVbInv(TMatrixD::kInverted,vVb); + TMatrixD rb(3,1); + rb(0,0) = diamondxyz[0]; + rb(1,0) = diamondxyz[1]; + rb(2,0) = diamondxyz[2]; + TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb); + + // input vertex + TMatrixD rv(3,1); + rv(0,0) = inVtx->GetXv(); + rv(1,0) = inVtx->GetYv(); + rv(2,0) = inVtx->GetZv(); + TMatrixD vV(3,3); + Double_t cov[6]; + inVtx->GetCovMatrix(cov); + vV(0,0) = cov[0]; + vV(0,1) = cov[1]; vV(1,0) = cov[1]; + vV(1,1) = cov[2]; + vV(0,2) = cov[3]; vV(2,0) = cov[3]; + vV(1,2) = cov[4]; vV(2,1) = cov[4]; + vV(2,2) = cov[5]; + TMatrixD vVInv(TMatrixD::kInverted,vV); + TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv); + + + TMatrixD sumWi = vVInv - vVbInv; + + + TMatrixD sumWiri = vVInvrv - vVbInvrb; + + TMatrixD rvnew(3,1); + TMatrixD vVnew(3,3); + + // new inverted of weights matrix + TMatrixD invsumWi(TMatrixD::kInverted,sumWi); + vVnew = invsumWi; + // new position of primary vertex + rvnew.Mult(vVnew,sumWiri); + + Double_t position[3]; + position[0] = rvnew(0,0); + position[1] = rvnew(1,0); + position[2] = rvnew(2,0); + cov[0] = vVnew(0,0); + cov[1] = vVnew(0,1); + cov[2] = vVnew(1,1); + cov[3] = vVnew(0,2); + cov[4] = vVnew(1,2); + cov[5] = vVnew(2,2); + + + Double_t chi2 = inVtx->GetChi2(); + + // diamond constribution to chi2 + TMatrixD deltar = rv; deltar -= rb; + TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar); + Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+ + deltar(1,0)*vVbInvdeltar(1,0)+ + deltar(2,0)*vVbInvdeltar(2,0); + // remove from total chi2 + chi2 -= chi2b; + + // store data in the vertex object + AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1); + outVtx->SetTitle("VertexerTracksNoConstraint"); + UShort_t *inindices = inVtx->GetIndices(); + Int_t nIndices = inVtx->GetNIndices(); + outVtx->SetIndices(nIndices,inindices); + + return outVtx; +} +//--------------------------------------------------------------------------- void AliVertexerTracks::SetCuts(Double_t *cuts) { // @@ -942,7 +1056,7 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights) { 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; iGetAlpha(); @@ -955,7 +1069,6 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights) sigmasq[2]=track1->GetSigmaZ2(); TMatrixD ri(3,1); TMatrixD wWi(3,3); - //if(!TrackToPoint(track1,ri,wWi)) {printf("WARNING\n");continue;} if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");} Double_t wmat[9]; Int_t iel=0; @@ -965,16 +1078,31 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights) 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; iGetEntriesFast(); + AliStrLine** lines2 = new AliStrLine* [knacc]; + for(Int_t i=0; iAt(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]; @@ -991,7 +1119,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t } for(Int_t i=0; iAt(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); @@ -1072,6 +1200,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t delete vectP1; return theVert; } + //--------------------------------------------------------------------------- Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t, TMatrixD &ri,TMatrixD &wWi, @@ -1209,8 +1338,8 @@ void AliVertexerTracks::TooFewTracks() } 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; } @@ -1540,7 +1669,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray, 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); @@ -1557,8 +1686,8 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray, } // clean up - delete [] indices; indices=NULL; - delete [] fIdSel; fIdSel=NULL; + if (indices) {delete indices; indices=NULL;} + delete fIdSel; fIdSel=NULL; fTrkArraySel.Delete(); return fCurrentVertex; @@ -1583,7 +1712,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray, VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate); - delete [] id; id=NULL; + delete id; id=NULL; return fCurrentVertex; }