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