LCFIVertex  0.7.2
vertexmass.cpp
1 #include <inc/vertex.h>
2 #include <inc/decaychain.h>
3 #include <inc/track.h>
4 #include <inc/jet.h>
5 #include <util/inc/matrix.h>
6 #include <util/inc/projection.h>
7 #include <algo/inc/vertexmass.h>
8 #include <vector>
9 #include <string>
10 #include <math.h>
11 #include <util/inc/string.h>
12 
13 //author Erik Devetak
14 
15 namespace vertex_lcfi
16 {
17 
19  {
20  _MaxMomentumAngle = 3;
21  _MaxKinematicCorrectionSigma = 2;
22  _MaxMomentumCorrection =2 ;
23 
24  _ParameterNames.push_back("MaxMomentumAngle");
25  _ParameterNames.push_back("MaxKinematicCorrectionSigma");
26  _ParameterNames.push_back("MaxMomentumCorrection");
27 
28  _ParameterValues.push_back(makeString(_MaxMomentumAngle));
29  _ParameterValues.push_back(makeString(_MaxKinematicCorrectionSigma));
30  _ParameterValues.push_back(makeString(_MaxMomentumCorrection));
31 
32  }
33  string VertexMass::name() const
34  {
35  return _Name;
36  }
37 
38  std::vector<string> VertexMass::parameterNames() const
39  {
40  return _ParameterNames;
41  }
42 
43  std::vector<string> VertexMass::parameterValues() const
44  {
45  _ParameterValues.clear();
46  _ParameterValues.push_back(makeString(_MaxMomentumAngle));
47  _ParameterValues.push_back(makeString(_MaxKinematicCorrectionSigma));
48  _ParameterValues.push_back(makeString(_MaxMomentumCorrection));
49  return _ParameterValues;
50  }
51 
52  void VertexMass::setStringParameter(const string & Parameter, const string & Value)
53  {
54  this->badParameter(Parameter);
55  }
56 
57  void VertexMass::setDoubleParameter(const string & Parameter, const double Value)
58  {
59  if (Parameter == "MaxMomentumAngle")
60  {
61  _MaxMomentumAngle = Value;
62  return;
63  }
64  if (Parameter == "MaxKinematicCorrectionSigma")
65  {
66  _MaxKinematicCorrectionSigma = Value;
67  return;
68  }
69  if (Parameter == "MaxMomentumCorrection")
70  {
71  _MaxMomentumCorrection = Value;
72  return;
73  }
74  this->badParameter(Parameter);
75  }
76 
77  void VertexMass::setPointerParameter(const string & Parameter, void * Value)
78  {
79  this->badParameter(Parameter);
80  }
81 
82  double VertexMass::ptvertexcons( Vector3* totalmomentum, Vertex* IPVertex , Vertex* TheVertex, double numberoftracks ) const
83  {
84 
85  double ptv = totalmomentum->mag2() - pow((totalmomentum->dot(TheVertex->position() - IPVertex->position())/(TheVertex->position().subtract(IPVertex->position())).mag()),2);
86 
87  if( numberoftracks < 1 ||ptv <0 )
88  return 0;
89 
90  return sqrt( ptv );
91  }
92 
93 
94 
95 
96  double VertexMass::calculateFor(DecayChain* MyDecayChain) const
97  {
98  int totaltracks = 0;
99  double totalenergy = 0;
100  double Rmass = 0;
101  int vertexcounter = 0;
102  int tempvertex =0;
103  int numberoftracks =0;
104  double MomentumCorr = 0;
105  double CorrectMass = 0;
106  double PtVertexCons = 0;
107 
108 
109  Vector3 totalmom(0,0,0);
110 
111  //Count tracks in vertices and decide the seed this is redundant if we have the track attach class algorithm preceding it,
112  //but it is still here in case we use something else.
113  for (std::vector<Vertex*>::const_iterator iVertex = (MyDecayChain->vertices().begin()); iVertex != MyDecayChain->vertices().end() ;++iVertex)
114  {
115 
116  numberoftracks = (*iVertex)->tracks( ).size();
117  if ( numberoftracks > 1 || vertexcounter > 0 )
118  {
119  tempvertex = vertexcounter;
120  }
121  vertexcounter++;
122  }
123 
124 
125  // calculate the total energy and momentum using all tracks ( note that in this case we are effectively using the track attachment before,
126  // which makes it so that these are tracks in the seed vertex + attached tracks
127  for(std::vector<Track*>::const_iterator iTrack = (MyDecayChain->allTracks().begin()); iTrack != MyDecayChain->allTracks().end() ;++iTrack)
128  {
129 
130  totalmom = totalmom.add( (*iTrack)->momentum());
131  totalenergy += sqrt((*iTrack)->momentum().mag2() + 0.139567*0.139567); //momentum + pion mass
132 
133  totaltracks++;
134  }
135 
136  // Calculate the mass
137  Rmass = totalenergy*totalenergy - totalmom.mag2();
138 
139  // This calculates p^2*(1-cos^2(theta)) where theta is the angle between the vertex axis and the momentum wrt. the origin)
140  PtVertexCons = ptvertexcons( &totalmom, MyDecayChain->vertices()[0] , MyDecayChain->vertices()[tempvertex], numberoftracks );
141 
142  if( tempvertex > 0 )
143  MomentumCorr = Ptcalc( MyDecayChain->vertices()[0], MyDecayChain->vertices()[tempvertex], &totalmom, _MaxKinematicCorrectionSigma );
144 
145 
146 
147  //Note here we are putting the real mass = to 0 if themomentum is not in the same direction as the devay and the momentum is high
148  if (PtVertexCons > _MaxMomentumAngle * Rmass)
149  Rmass = 0;
150 
151  if ( Rmass < 0.0 || totaltracks < 1)
152  Rmass = 0;
153  else
154  {
155  Rmass = sqrt(Rmass);
156  }
157 
158  CorrectMass = sqrt(Rmass*Rmass + MomentumCorr*MomentumCorr) + MomentumCorr;
159 
160 
161  //and Mind this call if Rmass = 0 we return 0
162  if (CorrectMass > _MaxMomentumCorrection * Rmass )
163  CorrectMass = _MaxMomentumCorrection * Rmass;
164 
165 
166  return CorrectMass;
167 
168  }
169 
170 
171  double VertexMass::Ptcalc( Vertex* IPVertex , Vertex* TheVertex, Vector3* momentum , float sigmax )
172  {
173 
174  //I should really look into symmetric matrices (eventually)!!!
175 
176  double totalmom;
177  double x0, y0, lambda, c1, c2, mu, sigma = 0;
178  double ptmin0, ptmax0, theta, dir, step;
179  int imm, n;
180  double cosi, sine, aa, bb, cc, det2;
181  double smin, amin, amax, temp, xmin, xmax, ymin, ymax ;
182  bool iterate = true ;
183  double p3norm =0;
184  double ptmin = 0;
185  double ptmax = 0;
186 
187  Vector3 location = TheVertex->position()- IPVertex->position();
188  Vector3 length(0,0,0);
189  Vector3 yvector(0,0,0);
190  Vector3 normom(0,0,0);
191  Vector3 normy0(0,0,0);
192 
193  ptmin = 1001;
194  ptmax = -1;
195 
196  Matrix3x3 errormatrix = TheVertex->positionError()+ IPVertex->positionError();
197  Matrix3x3 result;
198  Matrix3x3 rotation;
199  Matrix3x3 inverse;
200 
201  result.clear();
202 
203  totalmom = momentum->mag();
204 
205  lambda = momentum->dot(location)/(totalmom*totalmom);
206  length = momentum->mult(lambda);
207  x0 = length.mag();
208  yvector =location.subtract(length);
209  y0 = yvector.mag();
210 
211  normom = momentum->unit();
212  normy0 = yvector.unit();
213 
214  rotation(0,0)= normom.x();
215  rotation(0,1)= normom.y();
216  rotation(0,2)= normom.z();
217  rotation(1,0)= normy0.x();
218  rotation(1,1)= normy0.y();
219  rotation(1,2)= normy0.z();
220  rotation(2,0)= normom.y()*normy0.z() - normom.z()*normy0.y();
221  rotation(2,1)= normom.z()*normy0.x() - normom.x()*normy0.z();
222  rotation(2,2)= normom.x()*normy0.y() - normom.y()*normy0.x();
223 
224  errormatrix = prod( errormatrix,trans( rotation ));
225  errormatrix = prod( rotation , errormatrix );
226 
227  if( determinant(errormatrix) !=0 )
228  {
229  inverse = InvertMatrix ( errormatrix );
230 
231  }
232  else
233  {
234 
235  std::cerr << "Error matrix is singular - vertexmass.cpp:236" << std::endl;
236  return -1;
237  }
238 
239 
240  lambda = x0 + (y0*inverse(0,1))/inverse(0,0);
241  sigma = (inverse(0,0)*(lambda-x0) * (lambda-x0)) - (2*inverse(0,1)*(lambda-x0)*y0) + (y0*y0*inverse(1,1));
242 
243  if( sigma > 0 ) sigma = sqrt(sigma);
244  else sigma = 0;
245 
246  p3norm = sigma;
247 
248  if(sigma <= sigmax) ptmin = 0;
249 
250  c1 = inverse(0,1)*x0 + inverse(1,1)*y0;
251  c2 = inverse(0,2)*x0 + inverse(2,1)*y0;
252  lambda = ( c1*inverse(2,2) -c2*inverse(2,1) ) / ( inverse(1,1)*inverse(2,2)- inverse(0,2)*inverse(0,2) );
253  mu = ( c2 - inverse(2,1)*lambda ) / inverse(2,2);
254 
255  sigma = inverse(0,0) * x0 * x0 - 2 * inverse(0,1)* x0 * (lambda-y0)
256  + inverse(1,1) * (lambda -y0) * (lambda-y0) - 2 * inverse(0,2) * x0 * mu
257  + 2 * inverse(2,1) * (lambda-y0) * mu + inverse(2,2) * mu * mu;
258 
259  if( sigma > 0 ) sigma = sqrt(sigma);
260  else sigma = 0;
261 
262  if(sigma <= sigmax) ptmax = totalmom;
263 
264 
265  for( imm = 0; imm<2; imm++ )
266  {
267  if ( ptmin <= 0.00001 && imm == 0 ) continue;
268  if ( fabs( ptmax-totalmom ) <= 0.00001 && imm == 1 ) continue;
269 
270  //initializations for each loop
271  ptmin0 = 1000;
272  ptmax0 =0;
273  n = -1;
274  theta =0;
275  dir = 1;
276 
277  //determine step sizes for the kinematic correction of the decay axis
278  // this is an iterative process
279 
280  if( ( 0.001 *sigmax / y0 ) < 0.1 )
281  {
282  step = 0.001 * sigmax / y0;
283  smin = 0.00005 * sigmax /y0;
284  }
285  else
286  {
287  step =0.1;
288  smin =0.005;
289  }
290 
291  step = step * 5;
292 
293 
294  // iteration in correcting the the vertex axis with the aim to find the pt min. At the same time we also calculate the maximum pt
295  // correction using the same error matrices. However this value is not used in the future.
296  do
297  {
298  cosi = cos ( theta );
299  sine = sin ( theta );
300  iterate = true;
301 
302  aa = inverse(0,1) * inverse(0,1)* x0 * x0 - inverse(0,0) * inverse(1,1)* x0 * x0 + inverse(1,1) * (sigmax) * (sigmax);
303  bb = 2* inverse(0,0) * inverse(1,1) *x0 * y0 - 2 * inverse(0,1) * inverse(0,1) * x0 * y0 + 2 * inverse(0,1) * (sigmax) * (sigmax);
304  cc = inverse(0,1) * inverse(0,1) * y0 * y0 - inverse(0,0) * inverse(1,1) * y0 * y0 + inverse(0,0) * (sigmax) *(sigmax);
305 
306  //plane rotation
307 
308  aa *= cosi*cosi;
309 
310  aa += sine * ( inverse(0,2) * inverse(0,2) *x0 *x0 * sine + inverse(1,2) * inverse(1,2)*y0*y0* sine
311  + 2 * inverse(0,1) * inverse(0,2) * x0 * x0 * cosi + 2 * inverse(1,1) * inverse(0,2) * y0 * x0 * cosi
312  + 2 * inverse(0,2) * inverse(2,1) * x0 * y0 * sine - 2 * inverse(0,0) * inverse(2,1) * x0 * x0 * cosi
313  - inverse(0,0) * inverse(2,2) * x0 * x0 * sine - 2 * inverse(0,1) * inverse(2,1) * x0 * y0 * cosi
314  - 2 * inverse(0,1) * inverse (2,2) * x0 * y0 * sine - inverse(1,1) * inverse(2,2) * y0 * y0 * sine
315  + (2 * inverse(2,1) * cosi + inverse(2,2) * sine) * (sigmax) * (sigmax) );
316 
317 
318  bb *= cosi;
319  bb += 2 * sine * ( inverse(0,0)* inverse(1,2)* x0 * y0 + inverse(0,1) * inverse (1,2) *y0 * y0
320  - inverse(0,1)* inverse(0,2) * x0 * y0 - inverse(1,1) * inverse(0,2) *y0 * y0 + inverse(0,2) * (sigmax) * (sigmax) );
321 
322  det2 = bb*bb - 4*aa*cc;
323 
324  if( det2 > 0.0 )
325  {
326  amin = - (bb - sqrt(det2))/(2*aa);
327  amax = - (bb + sqrt(det2))/(2*aa);
328 
329  if(fabs(amin)> fabs(amax))
330  {
331  temp = amin;
332  amin = amax;
333  amax = temp;
334  }
335 
336  //coordinate calculation and mometnum corrections. Imm controls whether we are calculating the minimum or maximum correction
337  if(imm ==0 )
338  {
339  xmin = (inverse(0,0) * x0 + inverse(0,1) * y0 + ( inverse(0,1) * x0 + inverse(1,1) * y0) * amin * cosi
340  + ( inverse(0,2) * x0 + inverse(1,2) * y0 ) * amin * sine )/ ( inverse(0,0) + 2 * inverse(0,1) * amin * cosi
341  + inverse(1,1) * amin * amin * cosi * cosi + 2 * inverse (0,2) * amin * sine
342  + 2 * inverse (2,1) * amin * amin * cosi * sine + inverse(2,2) * amin * amin * sine *sine);
343 
344  ymin = amin * xmin;
345  ptmin = totalmom * fabs( ymin ) / sqrt(xmin * xmin+ ymin * ymin);
346  }
347  else
348  {
349 
350  xmax = (inverse(0,0) * x0 + inverse(0,1) * y0 + ( inverse(0,1) * x0 + inverse(1,1) * y0) * amax * cosi
351  + ( inverse(0,2) * x0 + inverse(1,2) * y0 ) * amax * sine )/ ( inverse(0,0) + 2 * inverse(0,1) * amax * cosi
352  + inverse(1,1) * amax * amax * cosi * cosi + 2 * inverse (0,2) * amax * sine
353  + 2 * inverse (2,1) * amax * amax * cosi * sine + inverse(2,2) * amax * amax * sine *sine);
354  ymax = amax*xmax;
355  ptmax = totalmom * fabs( ymax ) / sqrt(xmax * xmax+ ymax * ymax);
356  }
357 
358  }
359  else
360  {
361  if( imm == 0 ) ptmin = 1001;
362  else ptmax = -1;
363  }
364 
365  //here all the control processes.
366  //effectively stop when the step size is too small and keep the smalest pmin and the biggest pmax.
367  //if the new calculation of pt min is not smaller then change direction/angle of correction
368  //similar thinking for pt max
369  //notice that the step size decreases
370  if( ( imm == 0 && ptmin< ptmin0 ) || ( imm == 1 && ptmax > ptmax0 ) )
371  {
372  n++;
373  ptmin0 = ptmin;
374  ptmax0 = ptmax;
375  theta += dir * step;
376  }
377  else if( n == 0 )
378  {
379  n++;
380  theta -= 2* dir * step;
381  dir = - dir;
382 
383  }
384  else
385  {
386  theta -= dir * step;
387  if ( step < smin ) iterate = false;
388  else
389  {
390  step = step/5;
391  theta+= dir * step;
392  n = 0;
393  }
394  }
395 
396  }while(iterate == true );
397  }
398  if( ptmin > 999 )
399  {
400  ptmin = -1;
401  ptmax = -1;
402  }
403 
404  return ptmin;
405 
406  }
407 
408 }
void setDoubleParameter(const string &Parameter, const double Value)
Set Double Parameter.
Definition: vertexmass.cpp:57
VertexMass()
Default Constructor.
Definition: vertexmass.cpp:18
static double Ptcalc(Vertex *, Vertex *, Vector3 *, float)
Calculate vertex-constrained Pt correction.
Definition: vertexmass.cpp:171
const std::vector< Track * > & allTracks() const
All tracks contained in DecayChain.
Definition: decaychain.cpp:52
string name() const
Name.
Definition: vertexmass.cpp:33
double calculateFor(DecayChain *MyDecayChain) const
Run the algorithm on a jet.
Definition: vertexmass.cpp:96
void setStringParameter(const string &Parameter, const string &Value)
Set String Parameter.
Definition: vertexmass.cpp:52
std::vector< string > parameterValues() const
Parameter Values.
Definition: vertexmass.cpp:43
const std::vector< Vertex * > & vertices() const
Vertices contained in DecayChain.
Definition: decaychain.cpp:80
const Vector3 & position() const
Position.
std::vector< string > parameterNames() const
Parameter Names.
Definition: vertexmass.cpp:38
const SymMatrix3x3 & positionError() const
Position.
void setPointerParameter(const string &Parameter, void *Value)
Set Pointer Parameter.
Definition: vertexmass.cpp:77