Number of sigma pedestal cut increased to 4
[u/mrichter/AliRoot.git] / 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 }