[u/mrichter/AliRoot.git] / ITS / AliITSVertexerTracks.cxx

/**************************************************************************
 * 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 <TKey.h>
#include <TFile.h>
#include <TTree.h>
#include <TVector3.h>
#include <TMatrixD.h>
#include <TRandom.h>
//---- AliRoot headers -----
#include <AliRun.h>
#include "AliKalmanTrack.h"
#include "AliITSStrLine.h"
#include "AliITStrackV2.h"
#include "AliESDVertex.h"
#include "AliITSVertexerTracks.h"
#include "AliESD.h"
#include "AliESDtrack.h"


ClassImp(AliITSVertexerTracks)


//----------------------------------------------------------------------------
AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
//
// Default constructor
//
  fInFile  = 0;
  fOutFile = 0;
  SetVtxStart();
  SetMinTracks();
  fTrksToSkip = 0;
  fNTrksToSkip = 0;
  for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
}
//----------------------------------------------------------------------------
AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
                                           Double_t field,
                                           Int_t fEv,Int_t lEv,
                                           Double_t xStart,Double_t yStart) {
//
// Standard constructor
//
  fCurrentVertex = 0;
  fInFile  = inFile;
  fOutFile = outFile;
  SetFirstEvent(fEv);
  SetLastEvent(lEv);
  SetField(field);
  SetVtxStart(xStart,yStart);
  SetMinTracks();
  fTrksToSkip = 0;
  fNTrksToSkip = 0;
  for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
  SetDebug();
}
//----------------------------------------------------------------------------
AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn,
					   Double_t xStart,Double_t yStart)
                                          :AliITSVertexer(fn) {
//
// Alternative constructor
//
  fInFile  = 0;
  fOutFile = 0;
  SetField(field);
  SetVtxStart(xStart,yStart);
  SetMinTracks();
  fTrksToSkip = 0;
  fNTrksToSkip = 0;
  for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
}
//-----------------------------------------------------------------------------
AliITSVertexerTracks::~AliITSVertexerTracks() {
  // Default Destructor
  // The objects poited by the following pointers are not owned
  // by this class and are not deleted

  fCurrentVertex = 0;
  fInFile        = 0;
  fOutFile       = 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;

  TDirectory *curdir = 0;

  // 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();

    // write vertex to file
    TString vtxName = "Vertex_";
    vtxName += ev;
    fCurrentVertex->SetName(vtxName.Data()); 
    fCurrentVertex->SetTitle("VertexerTracks");
    //WriteCurrentVertex();
    curdir = gDirectory;
    fOutFile->cd();
    fCurrentVertex->Write();
    curdir->cd();
    fCurrentVertex = 0;
  } // loop over events

  return;
}
//---------------------------------------------------------------------------
void AliITSVertexerTracks::FindVerticesESD() {
//
// Vertices for all events from fFirstEvent to fLastEvent
//

  // Check if the conv. const. has been set
  if(!CheckField()) return;

  TDirectory *curdir = 0;

  fInFile->cd();
  TKey *key=0;
  TIter next(fInFile->GetListOfKeys());
  // loop on events in file
  while ((key=(TKey*)next())!=0) {
    AliESD *esdEvent=(AliESD*)key->ReadObj();
    if(!esdEvent) { 
      printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n"); 
      return; 
    }
    Int_t ev = (Int_t)esdEvent->GetEventNumber();
    if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; }
    if(ev % 100 == 0 || fDebug) 
      printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);

    FindVertexForCurrentEvent(esdEvent);

    if(!fCurrentVertex) {
      printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
      continue;
    }

    cout<<"VERTICE TROVATO\n";
    if(fDebug) fCurrentVertex->PrintStatus();

    // write the ESD to file
    curdir = gDirectory;
    Char_t ename[100];
    sprintf(ename,"%d",ev);
    fOutFile->cd();
    esdEvent->Dump();
    esdEvent->Write(ename,TObject::kOverwrite);
    curdir->cd();
    fCurrentVertex = 0;
  } // 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);

  return;
}
//----------------------------------------------------------------------------
AliESDVertex* 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);
  VertexFitter();
  if(fDebug) printf(" vertex fit completed\n");

  return fCurrentVertex;
}
//----------------------------------------------------------------------------
AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
{
//
// Vertex for current ESD event
//
  fCurrentVertex = 0;
  Double_t vtx[3],cvtx[6];

  // put tracks reco in ITS in a tree
  Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
  TTree *trkTree = new TTree("TreeT_ITS","its tracks");
  AliITStrackV2 *itstrack = 0;
  trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0);

  for(Int_t i=0; i<entr; i++) {
    AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
    if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
      { delete esdTrack; continue; }
    itstrack = new AliITStrackV2(*esdTrack);
    trkTree->Fill();
    itstrack = 0;
    delete esdTrack; 
  }
  delete itstrack;

  // preselect 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; }

  // Set initial vertex position from ESD
  esdEvent->GetVertex()->GetXYZ(vtx);
  SetVtxStart(vtx[0],vtx[1]);

  // VERTEX FINDER
  VertexFinder();

  // VERTEX FITTER
  ComputeMaxChi2PerTrack(nTrks);
  VertexFitter();
  if(fDebug) printf(" vertex fit completed\n");

  // store vertex information in ESD
  fCurrentVertex->GetXYZ(vtx);
  fCurrentVertex->GetCovMatrix(cvtx);
  esdEvent->SetVertex(fCurrentVertex);

  cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
  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; 
}
//---------------------------------------------------------------------------
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 AliESDVertex(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;
      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 AliESDVertex(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;
}
//----------------------------------------------------------------------------
AliESDVertex *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);
  VertexFitter();
  if(fDebug) printf(" vertex fit completed\n");

  return fCurrentVertex;
}
//----------------------------------------------------------------------------
Commit	Line	Data
cab6ff9b	1	/**************************************************************************
	2	* Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------
	17	// Implementation of the vertexer from tracks
	18	//
	19	// Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it
	20	// M.Masera, Torino, massimo.masera@to.infn.it
	21	//-----------------------------------------------------------------
	22
	23	//---- standard headers ----
	24	#include <Riostream.h>
	25	//---- Root headers --------
11ba84a4	26	#include <TKey.h>
cab6ff9b	27	#include <TFile.h>
	28	#include <TTree.h>
	29	#include <TVector3.h>
	30	#include <TMatrixD.h>
11ba84a4	31	#include <TRandom.h>
cab6ff9b	32	//---- AliRoot headers -----
	33	#include <AliRun.h>
	34	#include "AliKalmanTrack.h"
	35	#include "AliITSStrLine.h"
	36	#include "AliITStrackV2.h"
d681bb2d	37	#include "AliESDVertex.h"
cab6ff9b	38	#include "AliITSVertexerTracks.h"
11ba84a4	39	#include "AliESD.h"
11ba84a4	40	#include "AliESDtrack.h"
cab6ff9b	41
	42
	43	ClassImp(AliITSVertexerTracks)
	44
	45
	46	//----------------------------------------------------------------------------
	47	AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
	48	//
	49	// Default constructor
	50	//
11ba84a4	51	fInFile = 0;
11ba84a4	52	fOutFile = 0;
cab6ff9b	53	SetVtxStart();
cab6ff9b	54	SetMinTracks();
cab6ff9b	55	fTrksToSkip = 0;
	56	fNTrksToSkip = 0;
	57	for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
	58	}
cab6ff9b	59	//----------------------------------------------------------------------------
11ba84a4	60	AliITSVertexerTracks::AliITSVertexerTracks(TFile inFile,TFile outFile,
	61	Double_t field,
	62	Int_t fEv,Int_t lEv,
	63	Double_t xStart,Double_t yStart) {
cab6ff9b	64	//
	65	// Standard constructor
	66	//
11ba84a4	67	fCurrentVertex = 0;
	68	fInFile = inFile;
	69	fOutFile = outFile;
	70	SetFirstEvent(fEv);
	71	SetLastEvent(lEv);
cab6ff9b	72	SetField(field);
	73	SetVtxStart(xStart,yStart);
	74	SetMinTracks();
cab6ff9b	75	fTrksToSkip = 0;
cab6ff9b	76	fNTrksToSkip = 0;
11ba84a4	77	for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
11ba84a4	78	SetDebug();
cab6ff9b	79	}
cab6ff9b	80	//----------------------------------------------------------------------------
11ba84a4	81	AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn,
	82	Double_t xStart,Double_t yStart)
	83	:AliITSVertexer(fn) {
	84	//
	85	// Alternative constructor
	86	//
	87	fInFile = 0;
	88	fOutFile = 0;
	89	SetField(field);
	90	SetVtxStart(xStart,yStart);
	91	SetMinTracks();
	92	fTrksToSkip = 0;
	93	fNTrksToSkip = 0;
	94	for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
	95	}
	96	//-----------------------------------------------------------------------------
	97	AliITSVertexerTracks::~AliITSVertexerTracks() {
	98	// Default Destructor
	99	// The objects poited by the following pointers are not owned
	100	// by this class and are not deleted
	101
	102	fCurrentVertex = 0;
	103	fInFile = 0;
	104	fOutFile = 0;
	105	}
	106	//-----------------------------------------------------------------------------
cab6ff9b	107	Bool_t AliITSVertexerTracks::CheckField() const {
	108	//
	109	// Check if the conv. const. has been set
	110	//
	111	AliITStrackV2 t;
	112	Double_t cc = t.GetConvConst();
	113	Double_t field = 100./0.299792458/cc;
	114
	115	if(field<0.1 \|\| field>0.6) {
	116	printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
	117	return kFALSE;
	118	}
	119	printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field);
	120	return kTRUE;
	121	}
	122	//---------------------------------------------------------------------------
	123	void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
	124	//
	125	// Max. contr. to the chi2 has been tuned as a function of multiplicity
	126	//
	127	if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
	128	} else { fMaxChi2PerTrack = 100.; }
	129
	130	return;
	131	}
	132	//---------------------------------------------------------------------------
	133	void AliITSVertexerTracks::FindVertices() {
	134	//
	135	// Vertices for all events from fFirstEvent to fLastEvent
	136	//
	137
	138	// Check if the conv. const. has been set
	139	if(!CheckField()) return;
	140
11ba84a4	141	TDirectory *curdir = 0;
cab6ff9b	142
	143	// loop over events
	144	for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
	145	if(ev % 100 == 0 \|\| fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
	146
	147	FindVertexForCurrentEvent(ev);
	148
	149	if(!fCurrentVertex) {
	150	printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
	151	continue;
	152	}
	153
	154	if(fDebug) fCurrentVertex->PrintStatus();
11ba84a4	155
11ba84a4	156	// write vertex to file
d4221964	157	TString vtxName = "Vertex_";
d4221964	158	vtxName += ev;
11ba84a4	159	fCurrentVertex->SetName(vtxName.Data());
d4221964	160	fCurrentVertex->SetTitle("VertexerTracks");
11ba84a4	161	//WriteCurrentVertex();
	162	curdir = gDirectory;
	163	fOutFile->cd();
	164	fCurrentVertex->Write();
	165	curdir->cd();
	166	fCurrentVertex = 0;
	167	} // loop over events
	168
	169	return;
	170	}
	171	//---------------------------------------------------------------------------
	172	void AliITSVertexerTracks::FindVerticesESD() {
	173	//
	174	// Vertices for all events from fFirstEvent to fLastEvent
	175	//
	176
	177	// Check if the conv. const. has been set
	178	if(!CheckField()) return;
	179
	180	TDirectory *curdir = 0;
	181
	182	fInFile->cd();
	183	TKey *key=0;
	184	TIter next(fInFile->GetListOfKeys());
	185	// loop on events in file
	186	while ((key=(TKey*)next())!=0) {
	187	AliESD esdEvent=(AliESD)key->ReadObj();
	188	if(!esdEvent) {
	189	printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
	190	return;
	191	}
	192	Int_t ev = (Int_t)esdEvent->GetEventNumber();
	193	if(ev<fFirstEvent \|\| ev>fLastEvent) { delete esdEvent; continue; }
	194	if(ev % 100 == 0 \|\| fDebug)
	195	printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
	196
	197	FindVertexForCurrentEvent(esdEvent);
	198
	199	if(!fCurrentVertex) {
	200	printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
	201	continue;
	202	}
	203
	204	cout<<"VERTICE TROVATO\n";
	205	if(fDebug) fCurrentVertex->PrintStatus();
	206
	207	// write the ESD to file
	208	curdir = gDirectory;
	209	Char_t ename[100];
	210	sprintf(ename,"%d",ev);
	211	fOutFile->cd();
	212	esdEvent->Dump();
	213	esdEvent->Write(ename,TObject::kOverwrite);
	214	curdir->cd();
	215	fCurrentVertex = 0;
cab6ff9b	216	} // loop over events
	217
	218	return;
	219	}
	220	//----------------------------------------------------------------------------
	221	Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
	222	//
	223	// Propagate tracks to initial vertex position and store them in a TObjArray
	224	//
	225	Double_t maxd0rphi = 3.;
	226	Double_t alpha,xlStart,d0rphi;
	227	Int_t nTrks = 0;
	228	Bool_t skipThis;
	229
	230	Int_t nEntries = (Int_t)trkTree.GetEntries();
	231
	232	if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
	233	fTrkArray.Expand(nEntries);
	234
	235	if(fDebug) {
	236	printf(" PrepareTracks()\n");
	237	trkTree.Print();
	238	}
	239
	240	for(Int_t i=0; i<nEntries; i++) {
	241	// check tracks to skip
	242	skipThis = kFALSE;
	243	for(Int_t j=0; j<fNTrksToSkip; j++) {
	244	if(i==fTrksToSkip[j]) {
	245	if(fDebug) printf("skipping track: %d\n",i);
	246	skipThis = kTRUE;
	247	}
	248	}
	249	if(skipThis) continue;
	250
	251	AliITStrackV2 *itstrack = new AliITStrackV2;
	252	trkTree.SetBranchAddress("tracks",&itstrack);
	253	trkTree.GetEvent(i);
	254
	255
	256	// propagate track to vtxSeed
	257	alpha = itstrack->GetAlpha();
	258	xlStart = fNominalPos[0]TMath::Cos(alpha)+fNominalPos[1]TMath::Sin(alpha);
	259	itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
	260	itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed
	261
	262	// select tracks with d0rphi < maxd0rphi
	263	d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
	264	if(d0rphi > maxd0rphi) { delete itstrack; continue; }
	265
	266	fTrkArray.AddLast(itstrack);
	267
	268	nTrks++;
	269	}
	270
	271	if(fTrksToSkip) delete [] fTrksToSkip;
	272
	273	return nTrks;
	274	}
	275	//----------------------------------------------------------------------------
	276	void AliITSVertexerTracks::PrintStatus() const {
	277	//
	278	// Print status
	279	//
280	printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
281	printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
282	printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
283	printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
cab6ff9b	284
	285	return;
	286	}
	287	//----------------------------------------------------------------------------
d681bb2d	288	AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
cab6ff9b	289	//
	290	// Vertex for current event
	291	//
	292	fCurrentVertex = 0;
	293
	294	// get tree with tracks from input file
	295	TString treeName = "TreeT_ITS_";
	296	treeName += evnumb;
11ba84a4	297	TTree trkTree=(TTree)fInFile->Get(treeName.Data());
cab6ff9b	298	if(!trkTree) return fCurrentVertex;
	299
	300
	301	// get tracks and propagate them to initial vertex position
	302	Int_t nTrks = PrepareTracks(*trkTree);
	303	delete trkTree;
	304	if(fDebug) printf(" tracks prepared: %d\n",nTrks);
	305	if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
	306
	307	// VERTEX FINDER
	308	VertexFinder();
	309
	310	// VERTEX FITTER
	311	ComputeMaxChi2PerTrack(nTrks);
cab6ff9b	312	VertexFitter();
	313	if(fDebug) printf(" vertex fit completed\n");
	314
cab6ff9b	315	return fCurrentVertex;
cab6ff9b	316	}
cab6ff9b	317	//----------------------------------------------------------------------------
d681bb2d	318	AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
11ba84a4	319	{
cab6ff9b	320	//
11ba84a4	321	// Vertex for current ESD event
cab6ff9b	322	//
11ba84a4	323	fCurrentVertex = 0;
	324	Double_t vtx[3],cvtx[6];
	325
	326	// put tracks reco in ITS in a tree
	327	Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
	328	TTree *trkTree = new TTree("TreeT_ITS","its tracks");
	329	AliITStrackV2 *itstrack = 0;
	330	trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0);
	331
	332	for(Int_t i=0; i<entr; i++) {
	333	AliESDtrack esdTrack = (AliESDtrack)esdEvent->GetTrack(i);
	334	if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
	335	{ delete esdTrack; continue; }
	336	itstrack = new AliITStrackV2(*esdTrack);
	337	trkTree->Fill();
	338	itstrack = 0;
	339	delete esdTrack;
	340	}
	341	delete itstrack;
	342
	343	// preselect tracks and propagate them to initial vertex position
	344	Int_t nTrks = PrepareTracks(*trkTree);
	345	delete trkTree;
	346	if(fDebug) printf(" tracks prepared: %d\n",nTrks);
	347	if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
cab6ff9b	348
11ba84a4	349	// Set initial vertex position from ESD
2257f27e	350	esdEvent->GetVertex()->GetXYZ(vtx);
11ba84a4	351	SetVtxStart(vtx[0],vtx[1]);
cab6ff9b	352
11ba84a4	353	// VERTEX FINDER
11ba84a4	354	VertexFinder();
cab6ff9b	355
11ba84a4	356	// VERTEX FITTER
	357	ComputeMaxChi2PerTrack(nTrks);
	358	VertexFitter();
	359	if(fDebug) printf(" vertex fit completed\n");
cab6ff9b	360
11ba84a4	361	// store vertex information in ESD
	362	fCurrentVertex->GetXYZ(vtx);
	363	fCurrentVertex->GetCovMatrix(cvtx);
2257f27e	364	esdEvent->SetVertex(fCurrentVertex);
cab6ff9b	365
11ba84a4	366	cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
11ba84a4	367	return fCurrentVertex;
cab6ff9b	368	}
cab6ff9b	369	//---------------------------------------------------------------------------
11ba84a4	370	void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
cab6ff9b	371	//
11ba84a4	372	// Mark the tracks not ot be used in the vertex finding
cab6ff9b	373	//
11ba84a4	374	fNTrksToSkip = n;
	375	fTrksToSkip = new Int_t[n];
	376	for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
	377	return;
cab6ff9b	378	}
	379	//---------------------------------------------------------------------------
	380	void AliITSVertexerTracks::TooFewTracks() {
	381	//
	382	// When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
	383	// and the number of tracks to -1
	384	//
d681bb2d	385	fCurrentVertex = new AliESDVertex(0.,0.,-1);
cab6ff9b	386	return;
	387	}
	388	//---------------------------------------------------------------------------
	389	void AliITSVertexerTracks::VertexFinder() {
	390
	391	// Get estimate of vertex position in (x,y) from tracks DCA
	392	// Then this estimate is stored to the data member fInitPos
	393	// (previous values are overwritten)
	394
	395
	396
	397	/*
	398	***** TEMPORARY!!! FOR TEST ONLY!!! ********************************
	399
	400	fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
	401	fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
	402	*/
	403
	404	fInitPos[2] = 0.;
	405	for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
	406
	407	Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
	408
	409	Double_t aver[3]={0.,0.,0.};
	410	Int_t ncombi = 0;
	411	AliITStrackV2 *track1;
	412	AliITStrackV2 *track2;
	413	for(Int_t i=0; i<nacc; i++){
	414	track1 = (AliITStrackV2*)fTrkArray.At(i);
	415	if(fDebug>5){
	416	Double_t xv,par[5];
	417	track1->GetExternalParameters(xv,par);
	418	cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
	419	for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
	420	cout<<endl;
	421	}
	422	Double_t mom1[3];
	423	Double_t alpha = track1->GetAlpha();
	424	Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
	425	Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
	426	mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
	427	mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
	428	mom1[2] = TMath::Cos(theta);
	429
	430	Double_t pos1[3];
	431	Double_t mindist = TMath::Cos(alpha)fNominalPos[0]+TMath::Sin(alpha)fNominalPos[1];
	432	track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
	433	AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
	434	for(Int_t j=i+1; j<nacc; j++){
	435	track2 = (AliITStrackV2*)fTrkArray.At(j);
	436	Double_t mom2[3];
	437	alpha = track2->GetAlpha();
	438	azim = TMath::ASin(track2->GetSnp())+alpha;
	439	theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
	440	mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
	441	mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
	442	mom2[2] = TMath::Cos(theta);
	443	Double_t pos2[3];
	444	mindist = TMath::Cos(alpha)fNominalPos[0]+TMath::Sin(alpha)fNominalPos[1];
	445	track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
	446	AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
	447	Double_t crosspoint[3];
	448	Int_t retcode = line2->Cross(line1,crosspoint);
	449	if(retcode<0){
450	if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
451	if(fDebug>10)cout<<"bad intersection\n";
452	line1->PrintStatus();
453	line2->PrintStatus();
454	}
455	else {
456	ncombi++;
457	for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
458	if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
459	if(fDebug>10)cout<<"\n Cross point: ";
460	if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
461	}
462	delete line2;
463	}
464	delete line1;
465	}
466	if(ncombi>0){
467	for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
468	}
469	else {
470	Warning("VertexFinder","Finder did not succed");
471	}
472
473
474	//************************************************************************
475	return;
476
477	}
cab6ff9b	478	//---------------------------------------------------------------------------
	479	void AliITSVertexerTracks::VertexFitter() {
	480	//
	481	// The optimal estimate of the vertex position is given by a "weighted
	482	// average of tracks positions"
	483	// Original method: CMS Note 97/0051
	484	//
	485	if(fDebug) {
	486	printf(" VertexFitter(): start\n");
	487	PrintStatus();
	488	}
	489
	490
	491	Int_t i,j,k,step=0;
	492	TMatrixD rv(3,1);
	493	TMatrixD V(3,3);
	494	rv(0,0) = fInitPos[0];
	495	rv(1,0) = fInitPos[1];
	496	rv(2,0) = 0.;
	497	Double_t xlStart,alpha;
	498	Double_t rotAngle;
	499	Double_t cosRot,sinRot;
	500	Double_t cc[15];
	501	Int_t nUsedTrks;
	502	Double_t chi2,chi2i;
	503	Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
	504	AliITStrackV2 *t = 0;
	505	Int_t failed = 0;
	506
	507	Int_t *skipTrack = new Int_t[arrEntries];
	508	for(i=0; i<arrEntries; i++) skipTrack[i]=0;
	509
	510	// 3 steps:
	511	// 1st - first estimate of vtx using all tracks
	512	// 2nd - apply cut on chi2 max per track
	513	// 3rd - estimate of global chi2
	514	for(step=0; step<3; step++) {
	515	if(fDebug) printf(" step = %d\n",step);
	516	chi2 = 0.;
	517	nUsedTrks = 0;
	518
	519	TMatrixD SumWiri(3,1);
	520	TMatrixD SumWi(3,3);
	521	for(i=0; i<3; i++) {
	522	SumWiri(i,0) = 0.;
	523	for(j=0; j<3; j++) SumWi(j,i) = 0.;
	524	}
	525
	526	// loop on tracks
	527	for(k=0; k<arrEntries; k++) {
	528	if(skipTrack[k]) continue;
	529	// get track from track array
	530	t = (AliITStrackV2*)fTrkArray.At(k);
	531	alpha = t->GetAlpha();
	532	xlStart = fInitPos[0]TMath::Cos(alpha)+fInitPos[1]TMath::Sin(alpha);
	533	t->PropagateTo(xlStart,0.,0.); // to vtxSeed
11ba84a4	534	rotAngle = alpha;
cab6ff9b	535	if(alpha<0.) rotAngle += 2.*TMath::Pi();
	536	cosRot = TMath::Cos(rotAngle);
	537	sinRot = TMath::Sin(rotAngle);
	538
	539	// vector of track global coordinates
	540	TMatrixD ri(3,1);
	541	ri(0,0) = t->GetX()cosRot-t->GetY()sinRot;
	542	ri(1,0) = t->GetX()sinRot+t->GetY()cosRot;
	543	ri(2,0) = t->GetZ();
	544
	545	// matrix to go from global (x,y,z) to local (y,z);
	546	TMatrixD Qi(2,3);
	547	Qi(0,0) = -sinRot;
	548	Qi(0,1) = cosRot;
	549	Qi(0,2) = 0.;
	550	Qi(1,0) = 0.;
	551	Qi(1,1) = 0.;
	552	Qi(1,2) = 1.;
	553
	554	// covariance matrix of local (y,z) - inverted
	555	TMatrixD Ui(2,2);
	556	t->GetExternalCovariance(cc);
	557	Ui(0,0) = cc[0];
	558	Ui(0,1) = cc[1];
	559	Ui(1,0) = cc[1];
	560	Ui(1,1) = cc[2];
	561
	562	// weights matrix: Wi = QiT * UiInv * Qi
	563	if(Ui.Determinant() <= 0.) continue;
	564	TMatrixD UiInv(TMatrixD::kInverted,Ui);
	565	TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi);
	566	TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi);
	567
	568	// track chi2
	569	TMatrixD deltar = rv; deltar -= ri;
	570	TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar);
	571	chi2i = deltar(0,0)*Wideltar(0,0)+
	572	deltar(1,0)*Wideltar(1,0)+
	573	deltar(2,0)*Wideltar(2,0);
	574
	575
	576	if(step==1 && chi2i > fMaxChi2PerTrack) {
	577	skipTrack[k] = 1;
	578	continue;
	579	}
	580
	581	// add to total chi2
	582	chi2 += chi2i;
	583
	584	TMatrixD Wiri(Wi,TMatrixD::kMult,ri);
	585
	586	SumWiri += Wiri;
	587	SumWi += Wi;
	588
	589	nUsedTrks++;
	590	} // end loop on tracks
	591
	592	if(nUsedTrks < fMinTracks) {
	593	failed=1;
	594	continue;
	595	}
	596
	597	Double_t determinant = SumWi.Determinant();
	598	//cerr<<" determinant: "<<determinant<<endl;
599	if(determinant < 100.) {
600	printf("det(V) = 0\n");
601	failed=1;
602	continue;
603	}
604
605	// inverted of weights matrix
606	TMatrixD InvSumWi(TMatrixD::kInverted,SumWi);
607	V = InvSumWi;
608
609	// position of primary vertex
610	rv.Mult(V,SumWiri);
611
612	} // end loop on the 3 steps
613
614	delete [] skipTrack;
615	delete t;
616
617	if(failed) {
618	TooFewTracks();
619	return;
620	}
621
622	Double_t position[3];
623	position[0] = rv(0,0);
624	position[1] = rv(1,0);
625	position[2] = rv(2,0);
626	Double_t covmatrix[6];
627	covmatrix[0] = V(0,0);
628	covmatrix[1] = V(0,1);
629	covmatrix[2] = V(1,1);
630	covmatrix[3] = V(0,2);
631	covmatrix[4] = V(1,2);
632	covmatrix[5] = V(2,2);
633
634	// store data in the vertex object
d681bb2d	635	fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
cab6ff9b	636
	637	if(fDebug) {
	638	printf(" VertexFitter(): finish\n");
	639	printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
	640	fCurrentVertex->PrintStatus();
	641	}
	642
	643	return;
	644	}
	645	//----------------------------------------------------------------------------
d681bb2d	646	AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
cab6ff9b	647	//
	648	// Return vertex from tracks in trkTree
	649	//
	650	if(fCurrentVertex) fCurrentVertex = 0;
	651
	652	// get tracks and propagate them to initial vertex position
	653	Int_t nTrks = PrepareTracks(*(&trkTree));
	654	if(fDebug) printf(" tracks prepared: %d\n",nTrks);
	655	if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
	656
	657	// VERTEX FINDER
	658	VertexFinder();
	659
	660	// VERTEX FITTER
	661	ComputeMaxChi2PerTrack(nTrks);
cab6ff9b	662	VertexFitter();
	663	if(fDebug) printf(" vertex fit completed\n");
	664
	665	return fCurrentVertex;
	666	}
	667	//----------------------------------------------------------------------------
	668
	669
	670