STEER/AliKFVertex.cxx

   1 //----------------------------------------------------------------------------
   2 // Implementation of the AliKFVertex class
   3 // .
   4 // @author  S.Gorbunov, I.Kisel
   5 // @version 1.0
   6 // @since   13.05.07
   7 //
   8 // Class to reconstruct and store primary and secondary vertices
   9 // The method is described in CBM-SOFT note 2007-003,
  10 // ``Reconstruction of decayed particles based on the Kalman filter'',
  11 // http://www.gsi.de/documents/DOC-2007-May-14-1.pdf
  12 //
  13 // This class is ALICE interface to general mathematics in AliKFParticleCore
  14 //
  15 //  -= Copyright &copy ALICE HLT Group =-
  16 //____________________________________________________________________________
  17
  18
  19 #include "AliKFVertex.h"
  20
  21 ClassImp(AliKFVertex)
  22
  23
  24 AliKFVertex::AliKFVertex( const AliVVertex &vertex ): fIsConstrained(0)
  25 {
  26   // Constructor from ALICE VVertex
  27
  28   vertex.GetXYZ( fP );
  29   vertex.GetCovarianceMatrix( fC );
  30   fChi2 = vertex.GetChi2();
  31   fNDF = 2*vertex.GetNContributors() - 3;
  32   fQ = 0;
  33   fAtProductionVertex = 0;
  34   fIsLinearized = 0;
  35   fSFromDecay = 0;
  36 }
  37
  38 /*
  39 void     AliKFVertex::Print(Option_t* ) const
  40 {
  41   cout<<"AliKFVertex position:    "<<GetX()<<" "<<GetY()<<" "<<GetZ()<<endl;
  42   cout<<"AliKFVertex cov. matrix: "<<GetCovariance(0)<<endl;
  43   cout<<"                         "<<GetCovariance(1)<<" "<<GetCovariance(2)<<endl;
  44   cout<<"                         "<<GetCovariance(3)<<" "<<GetCovariance(4)<<" "<<GetCovariance(5)<<endl;
  45 }
  46   */
  47
  48 void AliKFVertex::SetBeamConstraint( Double_t x, Double_t y, Double_t z,
  49                                      Double_t errX, Double_t errY, Double_t errZ )
  50 {
  51   // Set beam constraint to the vertex
  52   fP[0] = x;
  53   fP[1] = y;
  54   fP[2] = z;
  55   fC[0] = errX*errX;
  56   fC[1] = 0;
  57   fC[2] = errY*errY;
  58   fC[3] = 0;
  59   fC[4] = 0;
  60   fC[5] = errZ*errZ;
  61   fIsConstrained = 1;
  62 }
  63
  64 void AliKFVertex::SetBeamConstraintOff()
  65 {
  66   fIsConstrained = 0;
  67 }
  68
  69 void AliKFVertex::ConstructPrimaryVertex( const AliKFParticle *vDaughters[],
  70                                           int NDaughters, Bool_t vtxFlag[],
  71                                           Double_t ChiCut  )
  72 {
  73   //* Primary vertex finder with simple rejection of outliers
  74
  75   if( NDaughters<2 ) return;
  76   double constrP[3]={fP[0], fP[1], fP[2]};
  77   double constrC[6]={fC[0], fC[1], fC[2], fC[3], fC[4], fC[5]};
  78
  79   Construct( vDaughters, NDaughters, 0, -1, fIsConstrained );
  80
  81   SetVtxGuess( fVtxGuess[0], fVtxGuess[1], fVtxGuess[2] );
  82
  83   for( int i=0; i<NDaughters; i++ ) vtxFlag[i] = 1;
  84
  85   Int_t nRest = NDaughters;
  86   while( nRest>2 )
  87     {
  88       Double_t worstChi = 0.;
  89       Int_t worstDaughter = 0;
  90       for( Int_t it=0; it<NDaughters; it++ ){
  91         if( !vtxFlag[it] ) continue;
  92         const AliKFParticle &p = *(vDaughters[it]);
  93         AliKFVertex tmp = *this - p;
  94         Double_t chi = p.GetDeviationFromVertex( tmp );
  95         if( worstChi < chi ){
  96           worstChi = chi;
  97           worstDaughter = it;
  98         }
  99       }
 100       if( worstChi < ChiCut ) break;
 101
 102       vtxFlag[worstDaughter] = 0;
 103       *this -= *(vDaughters[worstDaughter]);
 104       nRest--;
 105     }
 106
 107   if( nRest>=2 ){// final refit
 108     SetVtxGuess( fP[0], fP[1], fP[2] );
 109     if( fIsConstrained ){
 110       fP[0] = constrP[0];
 111       fP[1] = constrP[1];
 112       fP[2] = constrP[2];
 113       for( int i=0; i<6; i++ ) fC[i] = constrC[i];
 114     }
 115     int nDaughtersNew=0;
 116     const AliKFParticle **vDaughtersNew=new const AliKFParticle *[NDaughters];
 117     for( int i=0; i<NDaughters; i++ ){
 118       if( vtxFlag[i] )  vDaughtersNew[nDaughtersNew++] = vDaughters[i];
 119     }
 120     Construct( vDaughtersNew, nDaughtersNew, 0, -1, fIsConstrained );
 121     delete[] vDaughtersNew;
 122   }
 123
 124   if( nRest<=2 && GetChi2()>ChiCut*ChiCut*GetNDF() ){
 125     for( int i=0; i<NDaughters; i++ ) vtxFlag[i] = 0;
 126     fNDF = -3;
 127     fChi2 = 0;
 128   }
 129 }