X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=STEER%2FAliVertexerTracks.cxx;h=c88f2144ba39372f7c7f699c5941086d23c9c95b;hb=b3f25f64f3d7f48e0e8bdd3e186e8e0f5eac89c5;hp=de569baa2fb38c98431aa2923e651755c86b337d;hpb=919e537f584de6126f21a3b31c67d9ab7479c09d;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliVertexerTracks.cxx b/STEER/AliVertexerTracks.cxx index de569baa2fb..c88f2144ba3 100644 --- a/STEER/AliVertexerTracks.cxx +++ b/STEER/AliVertexerTracks.cxx @@ -35,6 +35,9 @@ #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" @@ -67,7 +70,8 @@ fITSrefit(kTRUE), fFiducialR(3.), fFiducialZ(30.), fnSigmaForUi00(1.5), -fAlgo(1) +fAlgo(1), +fAlgoIter0(4) { // // Default constructor @@ -100,14 +104,14 @@ fITSrefit(kTRUE), fFiducialR(3.), fFiducialZ(30.), fnSigmaForUi00(1.5), -fAlgo(1) +fAlgo(1), +fAlgoIter0(4) { // // Standard constructor // SetVtxStart(); SetVtxStartSigma(); - SetTPCMode(); } //----------------------------------------------------------------------------- AliVertexerTracks::~AliVertexerTracks() @@ -116,14 +120,14 @@ 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 @@ -131,11 +135,20 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent) // 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(nTrksGetTrack(i); + AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i); // check tracks to skip Bool_t skipThis = kFALSE; for(Int_t j=0; jGetID()==fTrksToSkip[j]) { + 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; + // 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); if(fMode==0) trkArrayOrig.Delete(); - delete [] idOrig; idOrig=NULL; + delete idOrig; idOrig=NULL; if(f) { f->Close(); delete f; f = NULL; @@ -195,6 +216,17 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent) 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; } //---------------------------------------------------------------------------- @@ -229,10 +261,18 @@ 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; - 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); @@ -261,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)); @@ -302,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 @@ -317,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; @@ -547,7 +587,13 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig, // loop on tracks for(Int_t i=0; iCharge()!=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())); @@ -656,7 +702,7 @@ Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track, // 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; @@ -669,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 @@ -726,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)+ @@ -763,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; @@ -779,7 +825,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, } } outVtx->SetIndices(nIndices,outindices); - delete [] outindices; + if (outindices) delete outindices; /* printf("Vertex before removing tracks:"); @@ -793,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) { // @@ -808,6 +954,8 @@ void AliVertexerTracks::SetCuts(Double_t *cuts) SetMinDetFitter(cuts[6]); SetMaxTgl(cuts[7]); SetFiducialRZ(cuts[8],cuts[9]); + fAlgo=(Int_t)(cuts[10]); + fAlgoIter0=(Int_t)(cuts[11]); return; } @@ -821,7 +969,9 @@ void AliVertexerTracks::SetITSMode(Double_t dcacut, 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 @@ -841,6 +991,8 @@ void AliVertexerTracks::SetITSMode(Double_t dcacut, SetMinDetFitter(mindetfitter); SetMaxTgl(maxtgl); SetFiducialRZ(fidR,fidZ); + fAlgo=finderAlgo; + fAlgoIter0=finderAlgoIter0; return; } @@ -854,7 +1006,9 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut, 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 @@ -870,6 +1024,8 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut, SetMinDetFitter(mindetfitter); SetMaxTgl(maxtgl); SetFiducialRZ(fidR,fidZ); + fAlgo=finderAlgo; + fAlgoIter0=finderAlgoIter0; return; } @@ -900,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(); @@ -913,7 +1069,7 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights) 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++){ @@ -922,17 +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]; @@ -949,9 +1119,10 @@ 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); line1->GetP0(p0); line1->GetCd(cd); line1->GetSigma2P0(sigmasq); @@ -1029,6 +1200,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t delete vectP1; return theVert; } + //--------------------------------------------------------------------------- Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t, TMatrixD &ri,TMatrixD &wWi, @@ -1166,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; } @@ -1451,30 +1623,43 @@ void AliVertexerTracks::VertexFitter() AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray, UShort_t *id, Bool_t optUseFitter, - Bool_t optPropagate) + Bool_t optPropagate, + Bool_t optUseDiamondConstraint) { // // Return vertex from tracks (AliExternalTrackParam) in array // fCurrentVertex = 0; - SetConstraintOff(); + // 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 && nTrksSelGetNContributors()>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); @@ -1514,16 +1699,18 @@ 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; } //---------------------------------------------------------------------------- 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 @@ -1538,9 +1725,9 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray, id[i] = (UShort_t)esdt->GetID(); } - VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate); + VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint); - delete [] id; id=NULL; + delete id; id=NULL; return fCurrentVertex; }