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;
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);
+ t = new AliExternalTrackParam(); t->CopyFromVTrack(track);
}
// use TOF info about bunch crossing
//
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) {
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");
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;
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++;
return outVtx;
}
//---------------------------------------------------------------------------
-void AliVertexerTracks::SetCuts(Double_t *cuts)
+void AliVertexerTracks::SetCuts(Double_t *cuts, Int_t ncuts)
{
//
// 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]);
+ if (ncuts>0) SetDCAcut(cuts[0]);
+ if (ncuts>1) SetDCAcutIter0(cuts[1]);
+ if (ncuts>2) SetMaxd0z0(cuts[2]);
+ if (ncuts>3) if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
+ 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 (cuts[12]>1.) SetMVTukey2(cuts[12]);
- if (cuts[13]>1.) SetMVSig2Ini(cuts[13]);
- if (cuts[14]>0.1) SetMVMaxSigma2(cuts[14]);
- if (cuts[15]>1e-5) SetMVMinSig2Red(cuts[15]);
- if (cuts[16]>1e-5) SetMVMinDst(cuts[16]);
- if (cuts[17]>0.5) SetMVScanStep(cuts[17]);
- SetMVMaxWghNtr(cuts[18]);
- SetMVFinalWBinary(cuts[19]>0);
- if (cuts[20]>20.) SetBCSpacing(int(cuts[20]));
+ if (ncuts>12) if (cuts[12]>1.) SetMVTukey2(cuts[12]);
+ if (ncuts>13) if (cuts[13]>1.) SetMVSig2Ini(cuts[13]);
+ if (ncuts>14) if (cuts[14]>0.1) SetMVMaxSigma2(cuts[14]);
+ if (ncuts>15) if (cuts[15]>1e-5) SetMVMinSig2Red(cuts[15]);
+ if (ncuts>16) if (cuts[16]>1e-5) SetMVMinDst(cuts[16]);
+ 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 (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
+ 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 || fClusterize) 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++){
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);
+ if(!line1) {printf("ERROR %d %d\n",i,knacc); continue;}
line1->GetP0(p0);
line1->GetCd(cd);
line1->GetSigma2P0(sigmasq);
}
}
AliDebug(1," Vertex finding completed\n");
+ Double_t vdispersion=fVert.GetDispersion();
// vertex fitter
if(optUseFitter) {
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
if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
fCurrentVertex->SetIndices(nindices,indices);
// set vertex title
- TString title="VertexerTracksNoConstraintMV";
- if(fConstraint) title="VertexerTracksWithConstraintMV";
+ TString title="VertexerTracksMVNoConstraint";
+ if(fConstraint) title="VertexerTracksMVWithConstraint";
fCurrentVertex->SetTitle(title.Data());
fMVVertices->AddLast(fCurrentVertex);
AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
double step = fMVScanStep>1 ? fMVScanStep : 1.;
double zmx = 3*TMath::Sqrt(fNominalCov[5]);
double zmn = -zmx;
- int nz = (zmx-zmn)/step; if (nz<1) nz=1;
+ int nz = TMath::Nint((zmx-zmn)/step); if (nz<1) nz=1;
double dz = (zmx-zmn)/nz;
int izStart=0;
AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
//
}
+//______________________________________________________
+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->GetZv(),
+ 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;
+
+}