New class to find primary vertex coordinates using reconstructed tracks (V2)
authormasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Feb 2003 16:10:55 +0000 (16:10 +0000)
committermasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Feb 2003 16:10:55 +0000 (16:10 +0000)
ITS/AliITSVertexerTracks.cxx [new file with mode: 0644]
ITS/AliITSVertexerTracks.h [new file with mode: 0644]

diff --git a/ITS/AliITSVertexerTracks.cxx b/ITS/AliITSVertexerTracks.cxx
new file mode 100644 (file)
index 0000000..fb8d7db
--- /dev/null
@@ -0,0 +1,644 @@
+/**************************************************************************
+ * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//    Implementation of the vertexer from tracks
+//
+// Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it
+//         M.Masera,  Torino, massimo.masera@to.infn.it
+//-----------------------------------------------------------------
+
+//---- standard headers ----
+#include <Riostream.h>
+//---- Root headers --------
+#include <TFile.h>
+#include <TTree.h>
+#include <TVector3.h>
+#include <TMatrixD.h>
+#include <TObjArray.h>
+#include <TRandom.h> // TEMPORARY !!!!!!!
+//---- AliRoot headers -----
+#include <AliRun.h>
+#include "AliKalmanTrack.h"
+#include "AliITSStrLine.h"
+#include "AliITStrackV2.h"
+#include "AliITSVertex.h"
+#include "AliITSVertexerTracks.h"
+
+
+ClassImp(AliITSVertexerTracks)
+
+
+//----------------------------------------------------------------------------
+AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
+//
+// Default constructor
+//
+
+  SetVtxStart();
+  SetMinTracks();
+  SetUseThrustFrame();
+  SetPhiThrust();
+  fTrksToSkip = 0;
+  fNTrksToSkip = 0;
+  for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
+}
+
+//----------------------------------------------------------------------------
+AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
+                                          Double_t field,
+                                          Double_t xStart,Double_t yStart,
+                                          Int_t useThFr)
+                                          :AliITSVertexer(inFile,outFile) {
+//
+// Standard constructor
+//
+
+  SetField(field);
+  SetVtxStart(xStart,yStart);
+  SetMinTracks();
+  SetUseThrustFrame(useThFr);
+  SetPhiThrust();
+  fTrksToSkip = 0;
+  fNTrksToSkip = 0;
+  for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
+}
+//----------------------------------------------------------------------------
+Bool_t AliITSVertexerTracks::CheckField() const {
+//
+// Check if the conv. const. has been set
+//
+  AliITStrackV2 t;
+  Double_t cc    = t.GetConvConst();
+  Double_t field = 100./0.299792458/cc;
+
+  if(field<0.1 || field>0.6) {
+    printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
+    return kFALSE;
+  }
+  printf("AliITSVertexerTracks::CheckField():  Using B = %3.1f T\n",field);
+  return kTRUE;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
+//
+// Max. contr. to the chi2 has been tuned as a function of multiplicity
+//
+  if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
+  } else { fMaxChi2PerTrack = 100.; }
+
+  return;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::FindVertices() {
+//
+// Vertices for all events from fFirstEvent to fLastEvent
+//
+
+  // Check if the conv. const. has been set
+  if(!CheckField()) return;
+
+
+  // loop over events
+  for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
+    if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
+
+    FindVertexForCurrentEvent(ev);
+
+    if(!fCurrentVertex) {
+      printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
+      continue;
+    }
+
+    if(fDebug) fCurrentVertex->PrintStatus();
+
+    WriteCurrentVertex();
+  } // loop over events
+
+  return;
+}
+//----------------------------------------------------------------------------
+Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
+//
+// Propagate tracks to initial vertex position and store them in a TObjArray
+//
+  Double_t maxd0rphi = 3.;
+  Double_t alpha,xlStart,d0rphi;
+  Int_t    nTrks    = 0;
+  Bool_t   skipThis;
+
+  Int_t    nEntries = (Int_t)trkTree.GetEntries();
+
+  if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
+  fTrkArray.Expand(nEntries);
+
+  if(fDebug) {
+    printf(" PrepareTracks()\n");
+    trkTree.Print();
+  }
+
+  for(Int_t i=0; i<nEntries; i++) {
+    // check tracks to skip
+    skipThis = kFALSE;
+    for(Int_t j=0; j<fNTrksToSkip; j++) { 
+      if(i==fTrksToSkip[j]) {
+       if(fDebug) printf("skipping track: %d\n",i);
+       skipThis = kTRUE;
+      }
+    }
+    if(skipThis) continue;
+
+    AliITStrackV2 *itstrack = new AliITStrackV2; 
+    trkTree.SetBranchAddress("tracks",&itstrack);
+    trkTree.GetEvent(i);
+
+
+    // propagate track to vtxSeed
+    alpha  = itstrack->GetAlpha();
+    xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
+    itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be) 
+    itstrack->PropagateTo(xlStart,0.,0.);   // to vtxSeed
+
+    // select tracks with d0rphi < maxd0rphi
+    d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
+    if(d0rphi > maxd0rphi) { delete itstrack; continue; }
+   
+    fTrkArray.AddLast(itstrack);
+
+    nTrks++; 
+  }
+
+  if(fTrksToSkip) delete [] fTrksToSkip;
+
+  return nTrks;
+} 
+//----------------------------------------------------------------------------
+void AliITSVertexerTracks::PrintStatus() const {
+//
+// Print status
+//
+  printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
+  printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
+  printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
+  printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
+  printf(" Using Thrust Frame: %d    fPhiThrust = %f\n",fUseThrustFrame,fPhiThrust);
+
+  return;
+}
+//----------------------------------------------------------------------------
+AliITSVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
+//
+// Vertex for current event
+//
+  fCurrentVertex = 0;
+
+  // get tree with tracks from input file
+  TString treeName = "TreeT_ITS_";
+  treeName += evnumb;
+  TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
+  if(!trkTree) return fCurrentVertex;
+
+
+  // get tracks and propagate them to initial vertex position
+  Int_t nTrks = PrepareTracks(*trkTree);
+  delete trkTree;
+  if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+  if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
+
+  // VERTEX FINDER
+  VertexFinder();
+
+  // VERTEX FITTER
+  ComputeMaxChi2PerTrack(nTrks);
+  if(fUseThrustFrame) ThrustFinderXY();
+  if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
+  VertexFitter();
+  if(fDebug) printf(" vertex fit completed\n");
+
+  TString vtxName;
+  vtxName = "VertexTracks_";
+  vtxName += evnumb;
+  fCurrentVertex->SetName(vtxName.Data());
+  return fCurrentVertex;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
+//
+// Mark the tracks not ot be used in the vertex finding
+//
+  fNTrksToSkip = n;
+  fTrksToSkip = new Int_t[n]; 
+  for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
+  return; 
+}
+//----------------------------------------------------------------------------
+Double_t AliITSVertexerTracks::SumPl(TTree &momTree,Double_t phi) const {
+//
+// Function to be maximized for thrust determination
+//
+  TVector3 u(1.,1.,0);
+  u.SetMag(1.);
+  u.SetPhi(phi);
+  Double_t pl;
+
+  Double_t sum = 0.;
+
+  TVector3* mom=0;
+  momTree.SetBranchAddress("momenta",&mom);
+  Int_t entries = (Int_t)momTree.GetEntries();
+
+  for(Int_t i=0; i<entries; i++) {
+    momTree.GetEvent(i);
+    pl = mom->Dot(u);
+    if(pl>0.) sum += pl;
+  } 
+
+  delete mom;
+
+  return sum;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::ThrustFinderXY() {
+//
+// This function looks for the thrust direction, \vec{u}, in the (x,y) plane.
+// The following function is maximized:
+// \Sum_{\vec{p}\cdot\vec{u}} \vec{p}\cdot\vec{u} / \Sum |\vec{p}| 
+// where \vec{p} = (p_x,p_y)
+//
+  Double_t pt,alpha,phi;
+  Double_t totPt;
+
+  // tree for thrust determination
+  TVector3 *ioMom = new TVector3;
+  TTree *t = new TTree("Tree_Momenta","Tree with momenta");  
+  t->Branch("momenta","TVector3",&ioMom);
+  totPt     = 0.;
+
+  AliITStrackV2 *itstrack = 0; 
+  Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
+
+  // loop on tracks
+  for(Int_t i=0; i<arrEntries; i++) {
+    itstrack = (AliITStrackV2*)fTrkArray.At(i);
+    // momentum of the track at the vertex
+    pt = 1./TMath::Abs(itstrack->Get1Pt());
+    alpha = itstrack->GetAlpha();
+    phi = alpha+TMath::ASin(itstrack->GetSnp());
+    ioMom->SetX(pt*TMath::Cos(phi)); 
+    ioMom->SetY(pt*TMath::Sin(phi));
+    ioMom->SetZ(0.);
+
+    totPt   += ioMom->Pt();
+    t->Fill();
+  } // end loop on tracks
+
+  Double_t tValue=0.,tPhi=0.;
+  Double_t maxSumPl = 0.;
+  Double_t thisSumPl;
+  Double_t dPhi;
+  Int_t nSteps,iStep;
+
+  phi = 0.;
+  nSteps = 100;
+  dPhi = 2.*TMath::Pi()/(Double_t)nSteps;
+
+  for(iStep=0; iStep<nSteps; iStep++) {
+    phi += dPhi;
+    thisSumPl = SumPl(*t,phi);
+    if(thisSumPl > maxSumPl) {
+      maxSumPl = thisSumPl;
+      tPhi = phi;
+    }
+  }
+
+  phi = tPhi-dPhi;
+  nSteps = 10;
+  dPhi /= 5.;
+
+  for(iStep=0; iStep<nSteps; iStep++) {
+    phi += dPhi;
+    thisSumPl = SumPl(*t,phi);
+    if(thisSumPl > maxSumPl) {
+      maxSumPl = thisSumPl;
+      tPhi = phi;
+    }
+  }
+
+  tValue = 2.*maxSumPl/totPt;
+  if(tPhi<0.) tPhi += 2.*TMath::Pi();
+  if(tPhi>2.*TMath::Pi()) tPhi -= 2.*TMath::Pi();
+
+  SetPhiThrust(tPhi);
+
+  delete t;
+  delete ioMom;
+
+  return;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::TooFewTracks() {
+//
+// When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
+// and the number of tracks to -1
+//
+  fCurrentVertex = new AliITSVertex(0.,0.,-1);
+  return;
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::VertexFinder() {
+
+  // Get estimate of vertex position in (x,y) from tracks DCA
+  // Then this estimate is stored to the data member fInitPos   
+  // (previous values are overwritten)
+
+
+  /*
+******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
+
+fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
+fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
+  */
+
+  fInitPos[2] = 0.;
+  for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
+
+  Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
+
+  Double_t aver[3]={0.,0.,0.};
+  Int_t ncombi = 0;
+  AliITStrackV2 *track1;
+  AliITStrackV2 *track2;
+  for(Int_t i=0; i<nacc; i++){
+    track1 = (AliITStrackV2*)fTrkArray.At(i);
+    if(fDebug>5){
+      Double_t xv,par[5];
+      track1->GetExternalParameters(xv,par);
+      cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
+      for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
+      cout<<endl;
+    }
+    Double_t mom1[3];
+    Double_t alpha = track1->GetAlpha();
+    Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
+    Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
+    mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
+    mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
+    mom1[2] = TMath::Cos(theta);
+
+    Double_t pos1[3];
+    Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+    track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
+    AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
+    for(Int_t j=i+1; j<nacc; j++){
+      track2 = (AliITStrackV2*)fTrkArray.At(j);
+      Double_t mom2[3];
+      alpha = track2->GetAlpha();
+      azim = TMath::ASin(track2->GetSnp())+alpha;
+      theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
+      mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
+      mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
+      mom2[2] = TMath::Cos(theta);
+      Double_t pos2[3];
+      mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+      track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
+      AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
+      Double_t crosspoint[3];
+      Int_t retcode = line2->Cross(line1,crosspoint);
+      if(retcode<0){
+       if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
+       if(fDebug>10)cout<<"bad intersection\n";
+       line1->PrintStatus();
+       line2->PrintStatus();
+      }
+      else {
+       ncombi++;
+       for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
+       if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
+       if(fDebug>10)cout<<"\n Cross point: ";
+       if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
+      }
+      delete line2;
+    }
+    delete line1;
+  }
+  if(ncombi>0){
+    for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
+  }
+  else {
+    Warning("VertexFinder","Finder did not succed");
+  }
+
+
+  //************************************************************************
+  return;
+
+}
+
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::VertexFitter() {
+//
+// The optimal estimate of the vertex position is given by a "weighted 
+// average of tracks positions"
+// Original method: CMS Note 97/0051
+//
+  if(fDebug) { 
+    printf(" VertexFitter(): start\n");
+    PrintStatus();
+  }
+
+
+  Int_t i,j,k,step=0;
+  TMatrixD rv(3,1);
+  TMatrixD V(3,3);
+  rv(0,0) = fInitPos[0];
+  rv(1,0) = fInitPos[1];
+  rv(2,0) = 0.;
+  Double_t xlStart,alpha;
+  Double_t rotAngle;
+  Double_t cosRot,sinRot;
+  Double_t cc[15];
+  Int_t nUsedTrks;
+  Double_t chi2,chi2i;
+  Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
+  AliITStrackV2 *t = 0;
+  Int_t failed = 0;
+
+  Int_t *skipTrack = new Int_t[arrEntries];
+  for(i=0; i<arrEntries; i++) skipTrack[i]=0;
+
+  // 3 steps:
+  // 1st - first estimate of vtx using all tracks
+  // 2nd - apply cut on chi2 max per track
+  // 3rd - estimate of global chi2
+  for(step=0; step<3; step++) {
+    if(fDebug) printf(" step = %d\n",step);
+    chi2 = 0.;
+    nUsedTrks = 0;
+
+    TMatrixD SumWiri(3,1);
+    TMatrixD SumWi(3,3);
+    for(i=0; i<3; i++) {
+      SumWiri(i,0) = 0.;
+      for(j=0; j<3; j++) SumWi(j,i) = 0.;
+    }
+
+    // loop on tracks  
+    for(k=0; k<arrEntries; k++) {
+      if(skipTrack[k]) continue;
+      // get track from track array
+      t = (AliITStrackV2*)fTrkArray.At(k);
+      alpha = t->GetAlpha();
+      xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
+      t->PropagateTo(xlStart,0.,0.);   // to vtxSeed
+      rotAngle = alpha-fPhiThrust;
+      if(alpha<0.) rotAngle += 2.*TMath::Pi();
+      cosRot = TMath::Cos(rotAngle);
+      sinRot = TMath::Sin(rotAngle);
+      
+      // vector of track global coordinates
+      TMatrixD ri(3,1);
+      ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
+      ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
+      ri(2,0) = t->GetZ();
+
+      // matrix to go from global (x,y,z) to local (y,z);
+      TMatrixD Qi(2,3);
+      Qi(0,0) = -sinRot;
+      Qi(0,1) = cosRot;
+      Qi(0,2) = 0.;
+      Qi(1,0) = 0.;
+      Qi(1,1) = 0.;
+      Qi(1,2) = 1.;
+
+      // covariance matrix of local (y,z) - inverted
+      TMatrixD Ui(2,2);
+      t->GetExternalCovariance(cc);
+      Ui(0,0) = cc[0];
+      Ui(0,1) = cc[1];
+      Ui(1,0) = cc[1];
+      Ui(1,1) = cc[2];
+
+      // weights matrix: Wi = QiT * UiInv * Qi
+      if(Ui.Determinant() <= 0.) continue;
+      TMatrixD UiInv(TMatrixD::kInverted,Ui);
+      TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi);
+      TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi);
+
+      // track chi2
+      TMatrixD deltar = rv; deltar -= ri;
+      TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar);
+      chi2i = deltar(0,0)*Wideltar(0,0)+
+              deltar(1,0)*Wideltar(1,0)+
+             deltar(2,0)*Wideltar(2,0);
+
+
+      if(step==1 && chi2i > fMaxChi2PerTrack) {
+       skipTrack[k] = 1;
+       continue;
+      }
+
+      // add to total chi2
+      chi2 += chi2i;
+
+      TMatrixD Wiri(Wi,TMatrixD::kMult,ri); 
+
+      SumWiri += Wiri;
+      SumWi   += Wi;
+
+      nUsedTrks++;
+    } // end loop on tracks
+
+    if(nUsedTrks < fMinTracks) {
+      failed=1;
+      continue;
+    }
+
+    Double_t determinant = SumWi.Determinant();
+    //cerr<<" determinant: "<<determinant<<endl;
+    if(determinant < 100.)  { 
+      printf("det(V) = 0\n");       
+      failed=1;
+      continue;
+    }
+
+    // inverted of weights matrix
+    TMatrixD InvSumWi(TMatrixD::kInverted,SumWi);
+    V = InvSumWi;
+     
+    // position of primary vertex
+    rv.Mult(V,SumWiri);
+
+  } // end loop on the 3 steps
+
+  delete [] skipTrack;
+  delete t;
+
+  if(failed) { 
+    TooFewTracks(); 
+    return; 
+  }
+
+  Double_t position[3];
+  position[0] = rv(0,0);
+  position[1] = rv(1,0);
+  position[2] = rv(2,0);
+  Double_t covmatrix[6];
+  covmatrix[0] = V(0,0);
+  covmatrix[1] = V(0,1);
+  covmatrix[2] = V(1,1);
+  covmatrix[3] = V(0,2);
+  covmatrix[4] = V(1,2);
+  covmatrix[5] = V(2,2);
+  
+  // store data in the vertex object
+  fCurrentVertex = new AliITSVertex(fPhiThrust,position,covmatrix,chi2,nUsedTrks);
+
+  if(fDebug) {
+    printf(" VertexFitter(): finish\n");
+    printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
+    fCurrentVertex->PrintStatus();
+  }
+
+  return;
+}
+//----------------------------------------------------------------------------
+AliITSVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
+//
+// Return vertex from tracks in trkTree
+//
+  if(fCurrentVertex) fCurrentVertex = 0;
+
+  // get tracks and propagate them to initial vertex position
+  Int_t nTrks = PrepareTracks(*(&trkTree));
+  if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+  if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
+
+  // VERTEX FINDER
+  VertexFinder();
+
+  // VERTEX FITTER
+  ComputeMaxChi2PerTrack(nTrks);
+  if(fUseThrustFrame) ThrustFinderXY();
+  if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
+  VertexFitter();
+  if(fDebug) printf(" vertex fit completed\n");
+
+  return fCurrentVertex;
+}
+//----------------------------------------------------------------------------
+
+
+
diff --git a/ITS/AliITSVertexerTracks.h b/ITS/AliITSVertexerTracks.h
new file mode 100644 (file)
index 0000000..54a1192
--- /dev/null
@@ -0,0 +1,89 @@
+#ifndef ALIITSVERTEXERTRACKS_H
+#define ALIITSVERTEXERTRACKS_H
+/* Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+//-------------------------------------------------------
+// Class for primary vertex determination with ITS tracks
+//
+//   Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it
+//           M.Masera,  Torino, massimo.masera@to.infn.it 
+//-------------------------------------------------------
+
+/*****************************************************************************
+ *                                                                           *
+ * This class determines the primary vertex position using ITS tracks.       *
+ * This is done in two steps:                                                *
+ * 1) Vertex Finding: a reasonable estimate of the vertex position is        *
+ *    obtained from a mean of "points of closest approach" between all       *
+ *    possible pairs of tracks.                                              *
+ * 2) Vertex Fitting: once tracks are propagated to the position given by    *
+ *    first step, the optimal estimate of the position of vertex is obtained *
+ *    from a weighted average of the track positions. A covariance           *
+ *    matrix and a chi2 for the vertex are given.                            *
+ *                                                                           *
+ *****************************************************************************/
+
+#include "AliKalmanTrack.h"
+#include "AliITSVertexer.h"
+
+class TTree; 
+class TFile; 
+class TObjArray; 
+class AliITSVertex; 
+
+class AliITSVertexerTracks : public AliITSVertexer {
+  
+ public:
+  // default constructor
+  AliITSVertexerTracks();  
+  // standard constructor     
+  AliITSVertexerTracks(TFile *inFile,TFile *outFile,Double_t field,
+                      Double_t xStart=0,Double_t yStart=0,Int_t useThFr=0); 
+  // destructor
+  virtual ~AliITSVertexerTracks() {}
+  // return vertex from the set of tracks in the tree
+  AliITSVertex *VertexOnTheFly(TTree &trkTree);
+  // computes the vertex for the current event
+  virtual AliITSVertex*    FindVertexForCurrentEvent(Int_t evnumb);
+  // computes the vertex for each event and stores it on fOutFile
+  virtual void  FindVertices();
+  virtual void  PrintStatus() const;
+  // computes the vertex for the current event    
+  void  SetField(Double_t field) const
+    { AliKalmanTrack::SetConvConst(100./0.299792458/field); return; }
+  void  SetMinTracks(Int_t n=3) { fMinTracks = n; return; }
+  void  SetSkipTracks(Int_t n,Int_t *skipped);
+  void  SetUseThrustFrame(Int_t utf=0) { fUseThrustFrame = utf; return; } 
+  void  SetVtxStart(Double_t x=0,Double_t y=0) 
+    { fNominalPos[0]=x; fNominalPos[1]=y; return; }
+  
+ private:
+  Double_t  fInitPos[3];      // vertex position after vertex finder
+  Double_t  fNominalPos[2];   // initial knowledge on vertex position
+  Int_t     fMinTracks;       // minimum number of tracks
+  Double_t  fMaxChi2PerTrack; // maximum contribition to the chi2 
+  Int_t     fUseThrustFrame;  // if !=0 vertex is given in thrust ref. frame 
+  Double_t  fPhiThrust;       // thrust direction
+  TObjArray fTrkArray;        // array with tracks to be processed
+  Int_t     *fTrksToSkip;     // tracks to be skipped for find and fit 
+  Int_t     fNTrksToSkip;     // number of tracks to be skipped 
+
+  Bool_t   CheckField() const; 
+  void     ComputeMaxChi2PerTrack(Int_t nTracks);
+  Int_t    PrepareTracks(TTree &trkTree);
+  void     SetPhiThrust(Double_t phi=0.) { fPhiThrust=phi; return; } 
+  Double_t SumPl(TTree &momTree,Double_t phi) const;
+  void     ThrustFinderXY();
+  void     TooFewTracks();
+  void     VertexFinder();
+  void     VertexFitter();
+
+  ClassDef(AliITSVertexerTracks,1) // 3D Vertexing with ITS tracks 
+};
+
+#endif
+
+
+