+ fCurrentVertex = 0;
+ if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
+}
+//----------------------------------------------------------------------------
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
+{
+//
+// Primary vertex for current ESD event
+// (Two iterations:
+// 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
+// + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
+// 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
+// All ESD tracks with inside the beam pipe are then propagated to found vertex
+//
+ fCurrentVertex = 0;
+
+ // accept 1-track case only if constraint is available
+ if(!fConstraint && fMinTracks==1) fMinTracks=2;
+
+ // read tracks from ESD
+ Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
+ if(nTrksTot<=0) {
+ if(fDebug) printf("TooFewTracks\n");
+ TooFewTracks(esdEvent);
+ return fCurrentVertex;
+ }
+
+ TTree *trkTree = new TTree("TreeT","tracks");
+ AliESDtrack *esdTrack = 0;
+ trkTree->Branch("tracks","AliESDtrack",&esdTrack);
+
+ Bool_t skipThis;
+ for(Int_t i=0; i<nTrksTot; i++) {
+ AliESDtrack *et = esdEvent->GetTrack(i);
+ esdTrack = new AliESDtrack(*et);
+ // check tracks to skip
+ skipThis = kFALSE;
+ for(Int_t j=0; j<fNTrksToSkip; j++) {
+ if(et->GetID()==fTrksToSkip[j]) {
+ if(fDebug) printf("skipping track: %d\n",i);
+ skipThis = kTRUE;
+ }
+ }
+ if(skipThis) {delete esdTrack;continue;}
+ if(fITSin) {
+ if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) {delete esdTrack;continue;}
+ if(fITSrefit && !(esdTrack->GetStatus()&AliESDtrack::kITSrefit)) {delete esdTrack;continue;}
+ Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
+ if(nclus<fMinITSClusters) {delete esdTrack;continue;}
+ }
+ trkTree->Fill();
+ delete esdTrack;
+ }
+
+ // If fConstraint=kFALSE
+ // run VertexFinder(1) to get rough estimate of initVertex (x,y)
+ if(!fConstraint) {
+ // fill fTrkArray, for VertexFinder()
+ if(!fTrkArray.IsEmpty()) fTrkArray.Delete();
+ PrepareTracks(*trkTree,0);
+ Double_t cutsave = fDCAcut; fDCAcut = 0.1; // 1 mm
+ VertexFinder(1); // using weights, cutting dca < fDCAcut
+ fDCAcut = cutsave;
+ fTrkArray.Delete();
+ 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]);
+ } else {
+ fNominalPos[0] = 0.;
+ fNominalPos[1] = 0.;
+ fNominalPos[2] = 0.;
+ if(fDebug) printf("No mean vertex and VertexFinder failed\n");
+ }
+ }
+
+ // TWO ITERATIONS:
+ //
+ // ITERATION 1
+ // propagate tracks to fNominalPos vertex
+ // preselect them:
+ // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
+ // else reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
+ // ITERATION 2
+ // propagate tracks to best between initVertex and fCurrentVertex
+ // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
+ // between initVertex and fCurrentVertex)
+ for(Int_t iter=0; iter<2; iter++) {
+ if(fOnlyFitter && iter==0) continue;
+ Int_t nTrksPrep = PrepareTracks(*trkTree,iter+1);
+ if(fDebug) printf(" tracks prepared - iteration %d: %d\n",iter+1,nTrksPrep);
+ if(nTrksPrep < fMinTracks) {
+ if(fDebug) printf("TooFewTracks\n");
+ TooFewTracks(esdEvent);
+ if(fDebug) fCurrentVertex->PrintStatus();
+ fTrkArray.Delete();
+ delete trkTree;
+ return fCurrentVertex;
+ }
+
+ // vertex finder
+ if(!fOnlyFitter) {
+ if(nTrksPrep==1){
+ if(fDebug) printf("Just one track\n");
+ OneTrackVertFinder();
+ }else{
+ switch (fAlgo) {
+ 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;
+ }
+ }
+ if(fDebug) printf(" Vertex finding completed\n");
+ }
+
+ // vertex fitter
+ VertexFitter(fConstraint);
+ if(fDebug) printf(" Vertex fit completed\n");
+ if(iter==0) fTrkArray.Delete();
+ } // end loop on the two iterations
+
+
+ // take true pos from SPD vertex in ESD and write it in tracks' vertex
+ Double_t tp[3];
+ esdEvent->GetVertex()->GetTruePos(tp);
+ fCurrentVertex->SetTruePos(tp);
+
+ if(fConstraint) {
+ if(fOnlyFitter) {
+ fCurrentVertex->SetTitle("VertexerTracksWithConstraintOnlyFitter");
+ } else {
+ fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
+ }
+ } else {
+ fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
+ }
+
+
+ // set indices of used tracks
+ UShort_t *indices = 0;
+ AliESDtrack *ett = 0;
+ if(fCurrentVertex->GetNContributors()>0) {
+ indices = new UShort_t[fCurrentVertex->GetNContributors()];
+ for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
+ ett = (AliESDtrack*)fTrkArray.At(jj);
+ indices[jj] = (UShort_t)ett->GetID();
+ }
+ fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
+ }
+ delete [] indices;
+
+ delete trkTree;
+
+ fTrkArray.Delete();
+
+ if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
+
+
+ if(fDebug) fCurrentVertex->PrintStatus();
+ if(fDebug) fCurrentVertex->PrintIndices();
+
+
+ return fCurrentVertex;
+}
+//------------------------------------------------------------------------
+Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
+{
+ //
+ Double_t det=matr[0][0]*matr[1][1]*matr[2][2]-matr[0][0]*matr[1][2]*matr[2][1]-matr[0][1]*matr[1][0]*matr[2][2]+matr[0][1]*matr[1][2]*matr[2][0]+matr[0][2]*matr[1][0]*matr[2][1]-matr[0][2]*matr[1][1]*matr[2][0];
+ return det;
+}
+//-------------------------------------------------------------------------
+void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
+{
+ //
+ Double_t x12=p0[0]-p1[0];
+ Double_t y12=p0[1]-p1[1];
+ Double_t z12=p0[2]-p1[2];
+ Double_t kk=x12*x12+y12*y12+z12*z12;
+ m[0][0]=2-2/kk*x12*x12;
+ m[0][1]=-2/kk*x12*y12;
+ m[0][2]=-2/kk*x12*z12;
+ m[1][0]=-2/kk*x12*y12;
+ 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][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;
+ d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
+
+}
+//--------------------------------------------------------------------------
+void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
+{
+ //
+ Double_t x12=p1[0]-p0[0];
+ Double_t y12=p1[1]-p0[1];
+ Double_t z12=p1[2]-p0[2];
+
+ Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
+
+ Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
+
+ Double_t cc[3];
+ cc[0]=-x12/sigmasq[0];
+ cc[1]=-y12/sigmasq[1];
+ cc[2]=-z12/sigmasq[2];
+
+ Double_t ww=(-p0[0]*x12*sigmasq[1]*sigmasq[2]-p0[1]*y12*sigmasq[0]*sigmasq[2]-p0[2]*z12*sigmasq[0]*sigmasq[1])/den;
+
+ Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
+
+ Double_t aa[3];
+ aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
+ aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
+ aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
+
+ m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
+ m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
+ m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
+
+ m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
+ m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
+ m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
+
+ m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
+ m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
+ m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
+
+ d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
+ d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
+ d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
+
+ }
+//--------------------------------------------------------------------------
+Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
+{
+ //
+ Double_t x12=p0[0]-p1[0];
+ Double_t y12=p0[1]-p1[1];
+ Double_t z12=p0[2]-p1[2];
+ 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);
+}
+//---------------------------------------------------------------------------
+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",fTrkArray.GetEntries());
+ AliESDtrack *track1;
+ track1 = (AliESDtrack*)fTrkArray.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);
+ AliStrLine *line1 = new AliStrLine(pos,dir);
+ Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
+ Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
+ AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
+ Double_t crosspoint[3]={0.,0.,0.};
+ Double_t sigma=999.;
+ Int_t nContrib=-1;
+ Int_t retcode = zeta->Cross(line1,crosspoint);
+ if(retcode>=0){
+ sigma=line1->GetDistFromPoint(crosspoint);
+ nContrib=1;
+ }
+ delete zeta;
+ delete line1;
+ fVert.SetXYZ(crosspoint);
+ fVert.SetDispersion(sigma);
+ fVert.SetNContributors(nContrib);
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::HelixVertexFinder()
+{
+ // Get estimate of vertex position in (x,y) from tracks DCA
+
+
+ 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)fTrkArray.GetEntriesFast();
+
+ Double_t aver[3]={0.,0.,0.};
+ Double_t averquad[3]={0.,0.,0.};
+ Double_t sigmaquad[3]={0.,0.,0.};
+ Double_t sigma=0;
+ Int_t ncombi = 0;
+ AliESDtrack *track1;
+ AliESDtrack *track2;
+ Double_t distCA;
+ Double_t x, par[5];
+ Double_t alpha, cs, sn;
+ Double_t crosspoint[3];
+ for(Int_t i=0; i<nacc; i++){
+ track1 = (AliESDtrack*)fTrkArray.At(i);
+
+
+ for(Int_t j=i+1; j<nacc; j++){
+ track2 = (AliESDtrack*)fTrkArray.At(j);
+
+ distCA=track2->PropagateToDCA(track1,field);
+ if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
+ track1->GetExternalParameters(x,par);
+ alpha=track1->GetAlpha();
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t x1=x*cs - par[0]*sn;
+ Double_t y1=x*sn + par[0]*cs;
+ Double_t z1=par[1];
+ Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
+ track2->GetExternalParameters(x,par);
+ alpha=track2->GetAlpha();
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t x2=x*cs - par[0]*sn;
+ Double_t y2=x*sn + par[0]*cs;
+ Double_t z2=par[1];
+ Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
+ Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
+ Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+ Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+ Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+ crosspoint[0]=wx1*x1 + wx2*x2;
+ crosspoint[1]=wy1*y1 + wy2*y2;
+ crosspoint[2]=wz1*z1 + wz2*z2;
+
+ ncombi++;
+ for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
+ for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
+ }
+ }
+
+ }
+ if(ncombi>0){
+ for(Int_t jj=0;jj<3;jj++){
+ initPos[jj] = aver[jj]/ncombi;
+ averquad[jj]/=ncombi;
+ sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
+ sigma+=sigmaquad[jj];
+ }
+ sigma=TMath::Sqrt(TMath::Abs(sigma));
+ }
+ else {
+ Warning("HelixVertexFinder","Finder did not succed");
+ sigma=999;
+ }
+ fVert.SetXYZ(initPos);
+ fVert.SetDispersion(sigma);
+ fVert.SetNContributors(ncombi);