1 #include <inc/vertex.h>
2 #include <inc/decaychain.h>
5 #include <util/inc/matrix.h>
6 #include <util/inc/projection.h>
7 #include <algo/inc/vertexmass.h>
11 #include <util/inc/string.h>
20 _MaxMomentumAngle = 3;
21 _MaxKinematicCorrectionSigma = 2;
22 _MaxMomentumCorrection =2 ;
24 _ParameterNames.push_back(
"MaxMomentumAngle");
25 _ParameterNames.push_back(
"MaxKinematicCorrectionSigma");
26 _ParameterNames.push_back(
"MaxMomentumCorrection");
28 _ParameterValues.push_back(makeString(_MaxMomentumAngle));
29 _ParameterValues.push_back(makeString(_MaxKinematicCorrectionSigma));
30 _ParameterValues.push_back(makeString(_MaxMomentumCorrection));
40 return _ParameterNames;
45 _ParameterValues.clear();
46 _ParameterValues.push_back(makeString(_MaxMomentumAngle));
47 _ParameterValues.push_back(makeString(_MaxKinematicCorrectionSigma));
48 _ParameterValues.push_back(makeString(_MaxMomentumCorrection));
49 return _ParameterValues;
54 this->badParameter(Parameter);
59 if (Parameter ==
"MaxMomentumAngle")
61 _MaxMomentumAngle = Value;
64 if (Parameter ==
"MaxKinematicCorrectionSigma")
66 _MaxKinematicCorrectionSigma = Value;
69 if (Parameter ==
"MaxMomentumCorrection")
71 _MaxMomentumCorrection = Value;
74 this->badParameter(Parameter);
79 this->badParameter(Parameter);
82 double VertexMass::ptvertexcons(
Vector3* totalmomentum,
Vertex* IPVertex ,
Vertex* TheVertex,
double numberoftracks )
const
85 double ptv = totalmomentum->mag2() - pow((totalmomentum->dot(TheVertex->
position() - IPVertex->
position())/(TheVertex->
position().subtract(IPVertex->
position())).mag()),2);
87 if( numberoftracks < 1 ||ptv <0 )
99 double totalenergy = 0;
101 int vertexcounter = 0;
103 int numberoftracks =0;
104 double MomentumCorr = 0;
105 double CorrectMass = 0;
106 double PtVertexCons = 0;
113 for (std::vector<Vertex*>::const_iterator iVertex = (MyDecayChain->
vertices().begin()); iVertex != MyDecayChain->
vertices().end() ;++iVertex)
116 numberoftracks = (*iVertex)->tracks( ).size();
117 if ( numberoftracks > 1 || vertexcounter > 0 )
119 tempvertex = vertexcounter;
127 for(std::vector<Track*>::const_iterator iTrack = (MyDecayChain->
allTracks().begin()); iTrack != MyDecayChain->
allTracks().end() ;++iTrack)
130 totalmom = totalmom.add( (*iTrack)->momentum());
131 totalenergy += sqrt((*iTrack)->momentum().mag2() + 0.139567*0.139567);
137 Rmass = totalenergy*totalenergy - totalmom.mag2();
140 PtVertexCons = ptvertexcons( &totalmom, MyDecayChain->
vertices()[0] , MyDecayChain->
vertices()[tempvertex], numberoftracks );
143 MomentumCorr =
Ptcalc( MyDecayChain->
vertices()[0], MyDecayChain->
vertices()[tempvertex], &totalmom, _MaxKinematicCorrectionSigma );
148 if (PtVertexCons > _MaxMomentumAngle * Rmass)
151 if ( Rmass < 0.0 || totaltracks < 1)
158 CorrectMass = sqrt(Rmass*Rmass + MomentumCorr*MomentumCorr) + MomentumCorr;
162 if (CorrectMass > _MaxMomentumCorrection * Rmass )
163 CorrectMass = _MaxMomentumCorrection * Rmass;
177 double x0, y0, lambda, c1, c2, mu, sigma = 0;
178 double ptmin0, ptmax0, theta, dir, step;
180 double cosi, sine, aa, bb, cc, det2;
181 double smin, amin, amax, temp, xmin, xmax, ymin, ymax ;
182 bool iterate = true ;
203 totalmom = momentum->mag();
205 lambda = momentum->dot(location)/(totalmom*totalmom);
206 length = momentum->mult(lambda);
208 yvector =location.subtract(length);
211 normom = momentum->unit();
212 normy0 = yvector.unit();
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();
224 errormatrix = prod( errormatrix,trans( rotation ));
225 errormatrix = prod( rotation , errormatrix );
227 if( determinant(errormatrix) !=0 )
229 inverse = InvertMatrix ( errormatrix );
235 std::cerr <<
"Error matrix is singular - vertexmass.cpp:236" << std::endl;
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));
243 if( sigma > 0 ) sigma = sqrt(sigma);
248 if(sigma <= sigmax) ptmin = 0;
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);
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;
259 if( sigma > 0 ) sigma = sqrt(sigma);
262 if(sigma <= sigmax) ptmax = totalmom;
265 for( imm = 0; imm<2; imm++ )
267 if ( ptmin <= 0.00001 && imm == 0 )
continue;
268 if ( fabs( ptmax-totalmom ) <= 0.00001 && imm == 1 )
continue;
280 if( ( 0.001 *sigmax / y0 ) < 0.1 )
282 step = 0.001 * sigmax / y0;
283 smin = 0.00005 * sigmax /y0;
298 cosi = cos ( theta );
299 sine = sin ( theta );
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);
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) );
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) );
322 det2 = bb*bb - 4*aa*cc;
326 amin = - (bb - sqrt(det2))/(2*aa);
327 amax = - (bb + sqrt(det2))/(2*aa);
329 if(fabs(amin)> fabs(amax))
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);
345 ptmin = totalmom * fabs( ymin ) / sqrt(xmin * xmin+ ymin * ymin);
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);
355 ptmax = totalmom * fabs( ymax ) / sqrt(xmax * xmax+ ymax * ymax);
361 if( imm == 0 ) ptmin = 1001;
370 if( ( imm == 0 && ptmin< ptmin0 ) || ( imm == 1 && ptmax > ptmax0 ) )
380 theta -= 2* dir * step;
387 if ( step < smin ) iterate =
false;
396 }
while(iterate ==
true );
void setDoubleParameter(const string &Parameter, const double Value)
Set Double Parameter.
VertexMass()
Default Constructor.
static double Ptcalc(Vertex *, Vertex *, Vector3 *, float)
Calculate vertex-constrained Pt correction.
const std::vector< Track * > & allTracks() const
All tracks contained in DecayChain.
string name() const
Name.
double calculateFor(DecayChain *MyDecayChain) const
Run the algorithm on a jet.
void setStringParameter(const string &Parameter, const string &Value)
Set String Parameter.
std::vector< string > parameterValues() const
Parameter Values.
const std::vector< Vertex * > & vertices() const
Vertices contained in DecayChain.
const Vector3 & position() const
Position.
std::vector< string > parameterNames() const
Parameter Names.
const SymMatrix3x3 & positionError() const
Position.
void setPointerParameter(const string &Parameter, void *Value)
Set Pointer Parameter.