laura is hosted by Hepforge, IPPP Durham
Laura++  v3r4
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauKinematics.hh
Go to the documentation of this file.
1 
2 /*
3 Copyright 2004 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
37 #ifndef LAU_KINEMATICS
38 #define LAU_KINEMATICS
39 
40 #include <vector>
41 
42 #include "TMath.h"
43 
44 
46 
47  public:
49 
58  LauKinematics(const Double_t m1, const Double_t m2, const Double_t m3, const Double_t mParent,
59  const Bool_t calcSquareDPCoords = kFALSE, const Bool_t symmetricalDP = kFALSE, const Bool_t fullySymmetricDP = kFALSE);
60 
62  virtual ~LauKinematics();
63 
65 
68  inline void squareDP( const Bool_t calcSquareDPCoords ) { squareDP_ = calcSquareDPCoords; }
69 
71 
74  inline Bool_t squareDP() const { return squareDP_; }
75 
77 
80  inline Bool_t gotSymmetricalDP() const { return (symmetricalDP_ && ! fullySymmetricDP_); }
81 
83 
86  inline Bool_t gotFullySymmetricDP() const { return fullySymmetricDP_; }
87 
89  inline void warningMessages(const Bool_t boolean) { warnings_ = boolean; }
90 
92 
98  void updateKinematics(const Double_t m13Sq, const Double_t m23Sq);
99 
101 
107  void updateSqDPKinematics(const Double_t mPrime, const Double_t thetaPrime);
108 
110 
114  void updateKinematicsFrom23(const Double_t m23, const Double_t c23);
115 
117 
121  void updateKinematicsFrom13(const Double_t m13, const Double_t c13);
122 
124 
128  void updateKinematicsFrom12(const Double_t m12, const Double_t c12);
129 
131 
136  Double_t calcSqDPJacobian(const Double_t mPrime, const Double_t thPrime) const;
137 
139 
142  Double_t calcSqDPJacobian() const;
143 
145 
149  void genFlatPhaseSpace(Double_t& m13Sq, Double_t& m23Sq) const;
150 
152 
156  void genFlatSqDP(Double_t& mPrime, Double_t& thetaPrime) const;
157 
159 
169  Bool_t withinDPLimits(const Double_t m13Sq, const Double_t m23Sq) const;
170 
172 
182  Bool_t withinDPLimits2(const Double_t m13Sq, const Double_t m23Sq) const;
183 
185 
190  Bool_t withinSqDPLimits(const Double_t mPrime, const Double_t thetaPrime) const;
191 
193 
198  Double_t calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const;
199 
201 
204  Double_t distanceFromDPCentre() const;
205 
207 
210  Double_t distanceFromDPCentre(const Double_t m13Sq, const Double_t m23Sq) const;
211 
213 
216  inline Double_t getm12() const {return m12_;}
218 
221  inline Double_t getm23() const {return m23_;}
223 
226  inline Double_t getm13() const {return m13_;}
227 
229 
232  inline Double_t getm12Sq() const {return m12Sq_;}
234 
237  inline Double_t getm23Sq() const {return m23Sq_;}
239 
242  inline Double_t getm13Sq() const {return m13Sq_;}
243 
245 
249  inline Double_t getc12() const {return c12_;}
251 
255  inline Double_t getc23() const {return c23_;}
257 
261  inline Double_t getc13() const {return c13_;}
262 
264 
267  inline Double_t getmPrime() const {return mPrime_;}
269 
272  inline Double_t getThetaPrime() const {return thetaPrime_;}
273 
275 
278  inline Double_t getmParent() const {return mParent_;}
280 
283  inline Double_t getmParentSq() const {return mParentSq_;}
284 
286 
292  inline Double_t getDPBoxArea() const {return (mSqMax_[1] - mSqMin_[1])*(mSqMax_[0] - mSqMin_[0]);}
293 
295 
299 
301 
305 
307 
310  inline Double_t getm1() const {return m1_;}
312 
315  inline Double_t getm2() const {return m2_;}
317 
320  inline Double_t getm3() const {return m3_;}
321 
323 
326  inline Double_t getm23Min() const {return TMath::Sqrt(mSqMin_[0]);}
328 
331  inline Double_t getm13Min() const {return TMath::Sqrt(mSqMin_[1]);}
333 
336  inline Double_t getm12Min() const {return TMath::Sqrt(mSqMin_[2]);}
337 
339 
342  inline Double_t getm23Max() const {return TMath::Sqrt(mSqMax_[0]);}
344 
347  inline Double_t getm13Max() const {return TMath::Sqrt(mSqMax_[1]);}
349 
352  inline Double_t getm12Max() const {return TMath::Sqrt(mSqMax_[2]);}
353 
355 
358  inline Double_t getm23SqMin() const {return mSqMin_[0];}
360 
363  inline Double_t getm13SqMin() const {return mSqMin_[1];}
365 
368  inline Double_t getm12SqMin() const {return mSqMin_[2];}
369 
371 
374  inline Double_t getm23SqMax() const {return mSqMax_[0];}
376 
379  inline Double_t getm13SqMax() const {return mSqMax_[1];}
381 
384  inline Double_t getm12SqMax() const {return mSqMax_[2];}
385 
387 
390  inline Double_t getp1_12() const {return p1_12_;}
392 
395  inline Double_t getp3_12() const {return p3_12_;}
397 
400  inline Double_t getp2_23() const {return p2_23_;}
402 
405  inline Double_t getp1_23() const {return p1_23_;}
407 
410  inline Double_t getp1_13() const {return p1_13_;}
412 
415  inline Double_t getp2_13() const {return p2_13_;}
416 
418 
421  inline Double_t getp1_Parent() const {return p1_Parent_;}
423 
426  inline Double_t getp2_Parent() const {return p2_Parent_;}
428 
431  inline Double_t getp3_Parent() const {return p3_Parent_;}
432 
434 
438  void drawDPContour(const Int_t orientation = 1323, const Int_t nbins = 100);
439 
441 
444  inline Double_t getcov12() const {return (mParentSq_ + m12Sq_ - m3Sq_)/(2.*mParent_*m12_);}
445 
447 
450  inline Double_t getcov13() const {return (mParentSq_ + m13Sq_ - m2Sq_)/(2.*mParent_*m13_);}
451 
453 
456  inline Double_t getcov23() const {return (mParentSq_ + m23Sq_ - m1Sq_)/(2.*mParent_*m23_);}
457 
458  protected:
460 
464  void updateMassSq_m12(const Double_t m12, const Double_t c12);
465 
467 
471  void updateMassSq_m23(const Double_t m23, const Double_t c23);
472 
474 
478  void updateMassSq_m13(const Double_t m13, const Double_t c13);
479 
481 
487  void updateMassSquares(const Double_t m13Sq, const Double_t m23Sq);
488 
490 
496  void updateSqDPMassSquares(const Double_t mPrime, const Double_t thetaPrime);
497 
499 
508  Double_t cFromM(const Double_t mijSq, const Double_t mikSq, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const;
509 
511 
520  Double_t mFromC(const Double_t mijSq, const Double_t cij, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const;
521 
523 
528  Double_t pCalc(const Double_t energy, const Double_t massSq) const;
529 
531 
534  Double_t genm13Sq() const;
536 
539  Double_t genm23Sq() const;
541 
544  Double_t genm12Sq() const;
545 
547  void calcm12Sq();
548 
550  void calcHelicities();
551 
553  void calcSqDPVars();
554 
556  void calcParentFrameMomenta();
557 
558  private:
560  LauKinematics(const LauKinematics& rhs);
561 
564 
566  const Bool_t symmetricalDP_;
568  const Bool_t fullySymmetricDP_;
569 
571  const Double_t m1_;
573  const Double_t m2_;
575  const Double_t m3_;
577  const Double_t mParent_;
578 
580  const Double_t m1Sq_;
582  const Double_t m2Sq_;
584  const Double_t m3Sq_;
586  const Double_t mParentSq_;
587 
589  std::vector<Double_t> mass_;
591  std::vector<Double_t> mMin_;
593  std::vector<Double_t> mMax_;
595  std::vector<Double_t> mDiff_;
596 
598  std::vector<Double_t> mSq_;
600  std::vector<Double_t> mSqMin_;
602  std::vector<Double_t> mSqMax_;
604  std::vector<Double_t> mSqDiff_;
605 
607  const Double_t mDTot_;
609  const Double_t massInt_;
611  const Double_t mSqDTot_;
612 
614  Double_t m12_;
616  Double_t m23_;
618  Double_t m13_;
619 
621  Double_t m12Sq_;
623  Double_t m23Sq_;
625  Double_t m13Sq_;
626 
628  Double_t c12_;
630  Double_t c23_;
632  Double_t c13_;
633 
635  Double_t mPrime_;
637  Double_t thetaPrime_;
638 
640  mutable Double_t qi_;
642  mutable Double_t qk_;
643 
645  Double_t p1_12_;
647  Double_t p3_12_;
649  Double_t p2_23_;
651  Double_t p1_23_;
653  Double_t p1_13_;
655  Double_t p2_13_;
656 
658  Double_t p1_Parent_;
660  Double_t p2_Parent_;
662  Double_t p3_Parent_;
663 
665  Bool_t squareDP_;
667  Bool_t warnings_;
668 
669  ClassDef(LauKinematics,0)
670 };
671 
672 Double_t dal(Double_t* x, Double_t* par);
673 
674 #endif
Bool_t withinDPLimits(const Double_t m13Sq, const Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot...
void updateMassSq_m12(const Double_t m12, const Double_t c12)
Update the variables m23Sq_ and m13Sq_ given the invariant mass m12 and the cosine of the helicity an...
Double_t p1_23_
Momentum of track 1 in 2-3 rest frame.
Double_t getm23Max() const
Get the m23 maximum defined as (mParent - m1)
Double_t getm13SqMin() const
Get the m13Sq minimum, (m1 + m3)^2.
const Double_t mParentSq_
Mass of parent particle squared.
Bool_t squareDP() const
Are the square Dalitz plot co-ordinates being calculated?
Bool_t gotSymmetricalDP() const
Is the DP symmetric?
Double_t getc23() const
Get the cosine of the helicity angle theta23.
Double_t getm23SqMin() const
Get the m23Sq minimum, (m2 + m3)^2.
void updateMassSq_m13(const Double_t m13, const Double_t c13)
Update the variables m12Sq_ and m23Sq_ given the invariant mass m13 and the cosine of the helicity an...
Double_t qi_
Momentum q of particle i.
void warningMessages(const Bool_t boolean)
Enable/disable warning messages.
Double_t getp2_Parent() const
Get the momentum of the track 2 in parent rest frame.
std::vector< Double_t > mSqMin_
Vector of the minimum mijSq values.
LauKinematics & operator=(const LauKinematics &rhs)
Copy assignment operator (not implemented)
Double_t getcov12() const
Get covariant factor in 12 axis.
void squareDP(const Bool_t calcSquareDPCoords)
Enable/disable the calculation of square Dalitz plot co-ordinates.
Double_t getm23() const
Get the m23 invariant mass.
Double_t m23_
Invariant mass m23.
Double_t getcov23() const
Get covariant factor in 23 axis.
Double_t p2_Parent_
Momentum of track 2 in parent rest frame.
Double_t getmPrime() const
Get m&#39; value.
Bool_t gotFullySymmetricDP() const
Is the DP fully symmetric?
std::vector< Double_t > mass_
Vector of daughter particles masses.
Double_t getc13() const
Get the cosine of the helicity angle theta13.
const Double_t m1Sq_
Mass of particle 1 squared.
Double_t getm3() const
Get the m3 mass.
Double_t cFromM(const Double_t mijSq, const Double_t mikSq, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const
General method to calculate the cos(helicity) variables from the masses of the particles.
Double_t getp2_13() const
Get the momentum of the track 2 in 13 rest frame.
Double_t getm1() const
Get the m1 mass.
Double_t getm23SqMax() const
Get the m23Sq maximum, (mParent - m1)^2.
Double_t m13_
Invariant mass m13.
Double_t p1_13_
Momentum of track 1 in 1-3 rest frame.
Bool_t warnings_
Enable/disable warning messages.
Double_t getp1_12() const
Get the momentum of the track 1 in 12 rest frame.
void updateSqDPKinematics(const Double_t mPrime, const Double_t thetaPrime)
Update all kinematic quantities based on the square DP co-ordinates m&#39; and theta&#39;.
Double_t getm13Min() const
Get the m13 minimum defined as (m1 + m3)
Double_t getm13SqMax() const
Get the m13Sq maximum, (mParent - m2)^2.
Double_t getm23Sq() const
Get the m23 invariant mass square.
Double_t distanceFromDPCentre() const
Calculate the distance from the currently set (m13Sq, m23Sq) point to the centre of the Dalitz plot (...
void updateMassSq_m23(const Double_t m23, const Double_t c23)
Update the variables m12Sq_ and m13Sq_ given the invariant mass m23 and the cosine of the helicity an...
const Double_t mDTot_
Sum of the daughter masses.
const Double_t mSqDTot_
Sum of the squares of the daughter masses.
Double_t getm12SqMin() const
Get the m12Sq minimum, (m1 + m2)^2.
Double_t mFromC(const Double_t mijSq, const Double_t cij, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const
General method to calculate mikSq given mijSq and cosHel_ij.
std::vector< Double_t > mMin_
Vector of the minimum mij values.
Double_t getm13Max() const
Get the m13 maximum defined as (mParent - m2)
Double_t genm23Sq() const
Randomly generate the invariant mass squared m23Sq.
Double_t getp1_Parent() const
Get the momentum of the track 1 in parent rest frame.
Double_t getm12Min() const
Get the m12 minimum defined as (m1 + m2)
void calcm12Sq()
Calculate m12Sq from m13Sq and m23Sq.
Double_t p2_23_
Momentum of track 2 in 2-3 rest frame.
Double_t thetaPrime_
theta&#39; co-ordinate
std::vector< Double_t > mSqDiff_
Vector of the difference between the mSqMax and mSqMin.
Double_t calcSqDPJacobian() const
Calculate the Jacobian for the transformation m23^2, m13^2 -&gt; m&#39;, theta&#39; (square DP) at the currently ...
void updateKinematicsFrom23(const Double_t m23, const Double_t c23)
Update all kinematic quantities based on m23 and cos(theta23)
Double_t getm2() const
Get the m2 mass.
Double_t getcov13() const
Get covariant factor in 13 axis.
Bool_t withinDPLimits2(const Double_t m13Sq, const Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot (alternat...
Double_t p2_13_
Momentum of track 2 in 1-3 rest frame.
void genFlatSqDP(Double_t &mPrime, Double_t &thetaPrime) const
Routine to generate events flat in the square Dalitz plot.
Double_t p1_Parent_
Momentum of track 1 in parent rest frame.
Double_t calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const
Calculate the third invariant mass square from the two provided (e.g. mjkSq from mijSq and mikSq) ...
Double_t getm12Max() const
Get the m12 maximum defined as (mParent - m3)
Double_t getp1_23() const
Get the momentum of the track 1 in 23 rest frame.
const Double_t m3Sq_
Mass of particle 3 squared.
Double_t getc12() const
Get the cosine of the helicity angle theta12.
Double_t getm13() const
Get the m13 invariant mass.
void flipAndUpdateKinematics()
Flips the DP variables m13^2 &lt;-&gt; m23^2 and recalculates all kinematic quantities. ...
Double_t getp1_13() const
Get the momentum of the track 1 in 13 rest frame.
Double_t getm23Min() const
Get the m23 minimum defined as (m2 + m3)
std::vector< Double_t > mMax_
Vector of the maximum mij values.
Double_t getmParent() const
Get parent mass.
Double_t m13Sq_
Invariant mass m13 squared.
void updateKinematics(const Double_t m13Sq, const Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
Double_t m12Sq_
Invariant mass m12 squared.
void updateSqDPMassSquares(const Double_t mPrime, const Double_t thetaPrime)
Update some kinematic quantities based on the square DP co-ordinates m&#39; and theta&#39;.
Double_t qk_
Momentum q of particle k.
Double_t p3_Parent_
Momentum of track 3 in parent rest frame.
void rotateAndUpdateKinematics()
Cyclically rotates the DP variables (m12 -&gt; m23, m23 -&gt; m13, m13 -&gt; m12) and recalculates all kinemat...
const Double_t m3_
Mass of particle 3.
Double_t genm13Sq() const
Randomly generate the invariant mass squared m13Sq.
Double_t getm12Sq() const
Get the m12 invariant mass square.
Double_t getp3_12() const
Get the momentum of the track 3 in 12 rest frame.
std::vector< Double_t > mSqMax_
Vector of the maximum mijSq values.
Double_t c23_
Cosine of the helicity angle theta23, which is defined as the angle between 1&amp;2 in the rest frame of ...
const Double_t massInt_
Mass difference between the parent particle and the sum of the daughter particles.
Double_t dal(Double_t *x, Double_t *par)
Double_t getm13Sq() const
Get the m13 invariant mass square.
void calcParentFrameMomenta()
Calculate the momenta of each daughter in the parent rest frame.
void updateMassSquares(const Double_t m13Sq, const Double_t m23Sq)
Update some kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
const Bool_t fullySymmetricDP_
Fully-symmetrical DP.
std::vector< Double_t > mDiff_
Vector of the difference between the mMax and mMin.
Double_t getm12() const
Get the m12 invariant mass.
Double_t getDPBoxArea() const
Get the box area defined from the kinematic bounds.
Bool_t squareDP_
Should we calculate the square DP co-ordinates or not?
Double_t m12_
Invariant mass m12.
Double_t getp2_23() const
Get the momentum of the track 2 in 23 rest frame.
Double_t c12_
Cosine of the helicity angle theta12, which is defined as the angle between 1&amp;3 in the rest frame of ...
const Double_t mParent_
Mass of parent particle.
void drawDPContour(const Int_t orientation=1323, const Int_t nbins=100)
Method to draw the Dalitz plot contours on the top of the histo previously drawn. ...
Double_t p3_12_
Momentum of track 3 in 1-2 rest frame.
Double_t getThetaPrime() const
Get theta&#39; value.
Double_t m23Sq_
Invariant mass m23 squared.
Double_t getm12SqMax() const
Get the m12Sq maximum, (mParent - m3)^2.
void genFlatPhaseSpace(Double_t &m13Sq, Double_t &m23Sq) const
Routine to generate events flat in phase-space.
virtual ~LauKinematics()
Destructor.
const Bool_t symmetricalDP_
Symmetrical DP.
Class for calculating 3-body kinematic quantities.
LauKinematics(const Double_t m1, const Double_t m2, const Double_t m3, const Double_t mParent, const Bool_t calcSquareDPCoords=kFALSE, const Bool_t symmetricalDP=kFALSE, const Bool_t fullySymmetricDP=kFALSE)
Constructor.
Double_t genm12Sq() const
Randomly generate the invariant mass squared m12Sq.
const Double_t m2_
Mass of particle 2.
Double_t c13_
Cosine of the helicity angle theta13, which is defined as the angle between 1&amp;2 in the rest frame of ...
void updateKinematicsFrom13(const Double_t m13, const Double_t c13)
Update all kinematic quantities based on m13 and cos(theta13)
void calcHelicities()
Calculate cosines of the helicity angles, momenta of daughters and bachelor in various ij rest frames...
const Double_t m2Sq_
Mass of particle 2 squared.
void updateKinematicsFrom12(const Double_t m12, const Double_t c12)
Update all kinematic quantities based on m12 and cos(theta12)
Double_t p1_12_
Momentum of track 1 in 1-2 rest frame.
Double_t pCalc(const Double_t energy, const Double_t massSq) const
General method to calculate the momentum of a particle, given its energy and invariant mass squared...
Double_t mPrime_
m&#39; co-ordinate
Bool_t withinSqDPLimits(const Double_t mPrime, const Double_t thetaPrime) const
Check whether a given (m&#39;,theta&#39;) point is within the kinematic limits of the Dalitz plot...
Double_t getmParentSq() const
Get parent mass squared.
Double_t getp3_Parent() const
Get the momentum of the track 3 in parent rest frame.
std::vector< Double_t > mSq_
Vector of daughter particles masses squared.
void calcSqDPVars()
Calculate the m&#39; and theta&#39; variables for the square Dalitz plot.
const Double_t m1_
Mass of particle 1.