fConstraint(kFALSE),
fOnlyFitter(kFALSE),
fMinTracks(1),
-fMinClusters(5),
+fMinClusters(3),
fDCAcut(0.1),
fDCAcutIter0(0.1),
fNSigma(3.),
fMVMaxWghNtr(10.),
fMVFinalWBinary(kTRUE),
fBCSpacing(50),
-fMVVertices(0)
+fMVVertices(0),
+fClusterize(kFALSE),
+fDeltaZCutForCluster(0.1),
+fnSigmaZCutForCluster(999999.)
{
//
// Default constructor
fConstraint(kFALSE),
fOnlyFitter(kFALSE),
fMinTracks(1),
-fMinClusters(5),
+fMinClusters(3),
fDCAcut(0.1),
fDCAcutIter0(0.1),
fNSigma(3.),
fMVMaxWghNtr(10.),
fMVFinalWBinary(kTRUE),
fBCSpacing(50),
-fMVVertices(0)
+fMVVertices(0),
+fClusterize(kFALSE),
+fDeltaZCutForCluster(0.1),
+fnSigmaZCutForCluster(999999.)
{
//
// Standard constructor
if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
TObjArray trkArrayOrig(nTrks);
UShort_t *idOrig = new UShort_t[nTrks];
+ Double_t *zTr = new Double_t[nTrks];
+ Double_t *err2zTr = new Double_t[nTrks];
Int_t nTrksOrig=0;
AliExternalTrackParam *t=0;
if(!PropagateTrackTo(t,radius)) continue;
}
} else { // AOD (only ITS mode)
+ if(track->GetID()<0) continue; // exclude global constrained and TPC only tracks (filter bits 128 and 512)
Int_t ncls0=0;
for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
if(ncls0 < fMinClusters) continue;
//
trkArrayOrig.AddLast(t);
idOrig[nTrksOrig]=(UShort_t)track->GetID();
+ zTr[nTrksOrig]=t->GetZ();
+ err2zTr[nTrksOrig]=t->GetSigmaZ2();
+
nTrksOrig++;
} // end loop on tracks
// call method that will reconstruct the vertex
- FindPrimaryVertex(&trkArrayOrig,idOrig);
-
+ if(fClusterize) FindAllVertices(nTrksOrig,&trkArrayOrig,zTr,err2zTr,idOrig);
+ else FindPrimaryVertex(&trkArrayOrig,idOrig);
if(!inputAOD) AnalyzePileUp((AliESDEvent*)vEvent);
if(fMode==0) trkArrayOrig.Delete();
delete [] idOrig; idOrig=NULL;
+ delete [] zTr; zTr=NULL;
+ delete [] err2zTr; err2zTr=NULL;
if(f) {
f->Close(); delete f; f = NULL;
esdt->SetVertexID(-1);
}
}
-
+
return fCurrentVertex;
}
//----------------------------------------------------------------------------
case 3: HelixVertexFinder(); break;
case 4: VertexFinder(1); break;
case 5: VertexFinder(0); break;
- default: printf("Wrong algorithm\n"); break;
+ default: {AliFatal(Form("Wrong seeder algorithm %d",fAlgoIter0));} break;
}
fDCAcut = cutsave;
if(fVert.GetNContributors()>0) {
fVert.GetXYZ(fNominalPos);
- fNominalPos[0] = fVert.GetXv();
- fNominalPos[1] = fVert.GetYv();
- fNominalPos[2] = fVert.GetZv();
+ fNominalPos[0] = fVert.GetX();
+ fNominalPos[1] = fVert.GetY();
+ fNominalPos[2] = fVert.GetZ();
AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
} else {
fNominalPos[0] = 0.;
case kVertexFinder1 : VertexFinder(1); break;
case kVertexFinder0 : VertexFinder(0); break;
case kMultiVertexer : FindVerticesMV(); multiMode = kTRUE; break;
- default: printf("Wrong algorithm"); break;
+ default: {AliFatal(Form("Wrong vertexer algorithm %d",fAlgo));} break;
}
}
AliDebug(1," Vertex finding completed");
}
fCurrentVertex->SetTitle(title.Data());
//
- AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
+ AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ(),fCurrentVertex->GetNContributors()));
}
// clean up
delete [] fIdSel; fIdSel=NULL;
m[1][1]=2-2/kk*y12*y12;
m[1][2]=-2/kk*y12*z12;
m[2][0]=-2/kk*x12*z12;
- m[2][1]=-2*y12*z12;
+ m[2][1]=-2/kk*y12*z12;
m[2][2]=2-2/kk*z12*z12;
d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
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]);
- 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])));
+ normdistx = TMath::Abs(fCurrentVertex->GetX()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
+ normdisty = TMath::Abs(fCurrentVertex->GetY()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
+ AliDebug(1,Form("normdistx %f %f %f",fCurrentVertex->GetX(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
+ AliDebug(1,Form("normdisty %f %f %f",fCurrentVertex->GetY(),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]))) {
printf("WARNING: result of tracks' removal will be only approximately correct");
TMatrixD rv(3,1);
- rv(0,0) = inVtx->GetXv();
- rv(1,0) = inVtx->GetYv();
- rv(2,0) = inVtx->GetZv();
+ rv(0,0) = inVtx->GetX();
+ rv(1,0) = inVtx->GetY();
+ rv(2,0) = inVtx->GetZ();
TMatrixD vV(3,3);
Double_t cov[6];
inVtx->GetCovMatrix(cov);
for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
Bool_t copyindex=kTRUE;
for(Int_t l=0; l<ntrks; l++) {
- if(inindices[k]==id[l]) copyindex=kFALSE;
+ if(inindices[k]==id[l]) {copyindex=kFALSE; break;}
}
if(copyindex) {
outindices[j] = inindices[k]; j++;
// input vertex
TMatrixD rv(3,1);
- rv(0,0) = inVtx->GetXv();
- rv(1,0) = inVtx->GetYv();
- rv(2,0) = inVtx->GetZv();
+ rv(0,0) = inVtx->GetX();
+ rv(1,0) = inVtx->GetY();
+ rv(2,0) = inVtx->GetZ();
TMatrixD vV(3,3);
Double_t cov[6];
inVtx->GetCovMatrix(cov);
if (ncuts>1) SetDCAcutIter0(cuts[1]);
if (ncuts>2) SetMaxd0z0(cuts[2]);
if (ncuts>3) if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
- if (ncuts>4) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
- if (ncuts>5) SetMinTracks((Int_t)(cuts[4]));
- if (ncuts>6) SetNSigmad0(cuts[5]);
- if (ncuts>7) SetMinDetFitter(cuts[6]);
- if (ncuts>8) SetMaxTgl(cuts[7]);
+ if (ncuts>3) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
+ if (ncuts>4) SetMinTracks((Int_t)(cuts[4]));
+ if (ncuts>5) SetNSigmad0(cuts[5]);
+ if (ncuts>6) SetMinDetFitter(cuts[6]);
+ if (ncuts>7) SetMaxTgl(cuts[7]);
if (ncuts>9) SetFiducialRZ(cuts[8],cuts[9]);
if (ncuts>10) fAlgo=(Int_t)(cuts[10]);
if (ncuts>11) fAlgoIter0=(Int_t)(cuts[11]);
if (ncuts>17) if (cuts[17]>0.5) SetMVScanStep(cuts[17]);
if (ncuts>18) SetMVMaxWghNtr(cuts[18]);
if (ncuts>19) SetMVFinalWBinary(cuts[19]>0);
- if (ncuts>20) if (cuts[20]>20.) SetBCSpacing(int(cuts[20]));
+ if (ncuts>20) SetBCSpacing(int(cuts[20]));
+ //
+ if (ncuts>21) if (cuts[21]>0.5) SetUseTrackClusterization(kTRUE);
+ if (ncuts>22) SetDeltaZCutForCluster(cuts[22]);
+ if (ncuts>23) SetnSigmaZCutForCluster(cuts[23]);
//
- if (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
+ if ( (fAlgo==kMultiVertexer || fClusterize) && fBCSpacing>0) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
else SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
//
return;
sigmasq[2]=track1->GetSigmaZ2();
TMatrixD ri(3,1);
TMatrixD wWi(3,3);
- if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
+ if(!TrackToPoint(track1,ri,wWi)) {
+ optUseWeights=kFALSE;
+ AliDebug(1,"WARNING: TrackToPoint failed");
+ }
Double_t wmat[9];
Int_t iel=0;
for(Int_t ia=0;ia<3;ia++){
fCurrentVertex->SetBC(currBC);
fVert = *fCurrentVertex; // RS
AliDebug(1," Vertex after fit:");
- AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
+ AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ(),fCurrentVertex->GetNContributors()));
AliDebug(1,"--- VertexFitter(): finish\n");
}
}
AliDebug(1," Vertex finding completed\n");
+ Double_t vdispersion=fVert.GetDispersion();
// vertex fitter
if(optUseFitter) {
VertexFitter();
} else {
- Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
+ Double_t position[3]={fVert.GetX(),fVert.GetY(),fVert.GetZ()};
Double_t covmatrix[6];
fVert.GetCovMatrix(covmatrix);
Double_t chi2=99999.;
Int_t nTrksUsed=fVert.GetNContributors();
fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
}
- fCurrentVertex->SetDispersion(fVert.GetDispersion());
+ fCurrentVertex->SetDispersion(vdispersion);
// set indices of used tracks and propagate track to found vertex
indices[jj] = fIdSel[jj];
t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
if(optPropagate && optUseFitter) {
- if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
+ if(TMath::Sqrt(fCurrentVertex->GetX()*fCurrentVertex->GetX()+fCurrentVertex->GetY()*fCurrentVertex->GetY())<3.) {
t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
AliDebug(1,Form("Track %d propagated to found vertex",jj));
} else {
// create indices
int ntrk = fTrkArraySel.GetEntries();
int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
+ if (nindices<1) {
+ delete fCurrentVertex;
+ fCurrentVertex = 0;
+ return kFALSE;
+ }
UShort_t *indices = 0;
if (nindices>0) indices = new UShort_t[nindices];
int nadded = 0;
t->SetBit(kBitAccounted);
indices[nadded++] = fIdSel[itr];
}
- if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
- fCurrentVertex->SetIndices(nindices,indices);
+ if (nadded!=nindices) {
+ printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
+ }
+ fCurrentVertex->SetIndices(nadded,indices);
// set vertex title
TString title="VertexerTracksMVNoConstraint";
if(fConstraint) title="VertexerTracksMVWithConstraint";
//
}
+//______________________________________________________
+void AliVertexerTracks::FindAllVertices(Int_t nTrksOrig,
+ const TObjArray *trkArrayOrig,
+ Double_t* zTr,
+ Double_t* err2zTr,
+ UShort_t* idOrig){
+
+ // clusterize tracks using z coordinates of intersection with beam axis
+ // and compute all vertices
+
+ UShort_t* posOrig=new UShort_t[nTrksOrig];
+ for(Int_t iTr=0; iTr<nTrksOrig; iTr++) posOrig[iTr]=iTr;
+
+
+ // sort points along Z
+ AliDebug(1,Form("Sort points along Z, used tracks %d",nTrksOrig));
+ for(Int_t iTr1=0; iTr1<nTrksOrig; iTr1++){
+ for(Int_t iTr2=iTr1+1; iTr2<nTrksOrig; iTr2++){
+ if(zTr[iTr1]>zTr[iTr2]){
+ Double_t tmpz=zTr[iTr2];
+ Double_t tmperr=err2zTr[iTr2];
+ UShort_t tmppos=posOrig[iTr2];
+ UShort_t tmpid=idOrig[iTr2];
+ zTr[iTr2]=zTr[iTr1];
+ err2zTr[iTr2]=err2zTr[iTr1];
+ posOrig[iTr2]=posOrig[iTr1];
+ idOrig[iTr2]=idOrig[iTr1];
+ zTr[iTr1]=tmpz;
+ err2zTr[iTr1]=tmperr;
+ idOrig[iTr1]=tmpid;
+ posOrig[iTr1]=tmppos;
+ }
+ }
+ }
+
+ // clusterize
+ Int_t nClusters=0;
+ Int_t* firstTr=new Int_t[nTrksOrig];
+ Int_t* lastTr=new Int_t[nTrksOrig];
+
+ firstTr[0]=0;
+ for(Int_t iTr=0; iTr<nTrksOrig-1; iTr++){
+ Double_t distz=zTr[iTr+1]-zTr[iTr];
+ Double_t errdistz=TMath::Sqrt(err2zTr[iTr+1]+err2zTr[iTr]);
+ if(errdistz<=0.000001) errdistz=0.000001;
+ if(distz>fDeltaZCutForCluster || (distz/errdistz)>fnSigmaZCutForCluster){
+ lastTr[nClusters]=iTr;
+ firstTr[nClusters+1]=iTr+1;
+ nClusters++;
+ }
+ }
+ lastTr[nClusters]=nTrksOrig-1;
+
+ // Compute vertices
+ AliDebug(1,Form("Number of found clusters %d",nClusters+1));
+ Int_t nFoundVertices=0;
+
+ if (!fMVVertices) fMVVertices = new TObjArray(nClusters+1);
+
+ fMVVertices->Clear();
+ TObjArray cluTrackArr(nTrksOrig);
+ UShort_t *idTrkClu=new UShort_t[nTrksOrig];
+ // Int_t maxContr=0;
+ // Int_t maxPos=-1;
+
+ for(Int_t iClu=0; iClu<=nClusters; iClu++){
+ Int_t nCluTracks=lastTr[iClu]-firstTr[iClu]+1;
+ cluTrackArr.Clear();
+ AliDebug(1,Form(" Vertex #%d tracks %d first tr %d last track %d",iClu,nCluTracks,firstTr[iClu],lastTr[iClu]));
+ Int_t nSelTr=0;
+ for(Int_t iTr=firstTr[iClu]; iTr<=lastTr[iClu]; iTr++){
+ AliExternalTrackParam* t=(AliExternalTrackParam*)trkArrayOrig->At(posOrig[iTr]);
+ if(TMath::Abs(t->GetZ()-zTr[iTr])>0.00001){
+ AliError(Form("Clu %d Track %d zTrack=%f zVec=%f\n",iClu,iTr,t->GetZ(),zTr[iTr]));
+ }
+ cluTrackArr.AddAt(t,nSelTr);
+ idTrkClu[nSelTr]=idOrig[iTr];
+ AliDebug(1,Form(" Add track %d: id %d, z=%f",iTr,idOrig[iTr],zTr[iTr]));
+ nSelTr++;
+ }
+ AliESDVertex* vert=FindPrimaryVertex(&cluTrackArr,idTrkClu);
+ AliDebug(1,Form("Found vertex in z=%f with %d contributors",vert->GetZ(),
+ vert->GetNContributors()));
+
+ fCurrentVertex=0;
+ if(vert->GetNContributors()>0){
+ nFoundVertices++;
+ fMVVertices->AddLast(vert);
+ }
+ // if(vert->GetNContributors()>maxContr){
+ // maxContr=vert->GetNContributors();
+ // maxPos=nFoundVertices-1;
+ // }
+ }
+
+ AliDebug(1,Form("Number of found vertices %d (%d)",nFoundVertices,fMVVertices->GetEntriesFast()));
+ // if(maxPos>=0 && maxContr>0){
+ // AliESDVertex* vtxMax=(AliESDVertex*)fMVVertices->At(maxPos);
+ // if(fCurrentVertex) delete fCurrentVertex;
+ // fCurrentVertex=new AliESDVertex(*vtxMax);
+ // }
+
+ delete [] firstTr;
+ delete [] lastTr;
+ delete [] idTrkClu;
+ delete [] posOrig;
+
+ return;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff