laura is hosted by Hepforge, IPPP Durham
Laura++  v3r2
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauKinematics.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <iostream>
16 
17 #include "TF2.h"
18 #include "TMath.h"
19 #include "TRandom.h"
20 
21 #include "LauConstants.hh"
22 #include "LauKinematics.hh"
23 #include "LauRandom.hh"
24 
26 
27 
28 LauKinematics::LauKinematics(const Double_t m1, const Double_t m2, const Double_t m3, const Double_t mParent, const Bool_t calcSquareDPCoords, const Bool_t symmetricalDP, const Bool_t fullySymmetricDP) :
29  symmetricalDP_(symmetricalDP),
30  fullySymmetricDP_(fullySymmetricDP),
31  m1_(m1), m2_(m2), m3_(m3), mParent_(mParent),
32  m1Sq_(m1*m1), m2Sq_(m2*m2), m3Sq_(m3*m3), mParentSq_(mParent*mParent),
33  mDTot_(m1 + m2 + m3),
34  massInt_(mParent - (m1+m2+m3)),
35  mSqDTot_(m1*m1 + m2*m2 + m3*m3),
36  m12_(0.0), m23_(0.0), m13_(0.0),
37  m12Sq_(0.0), m23Sq_(0.0), m13Sq_(0.0),
38  c12_(0.0), c23_(0.0), c13_(0.0),
39  mPrime_(0.0), thetaPrime_(0.0),
40  qi_(0.0), qk_(0.0),
41  p1_12_(0.0), p3_12_(0.0),
42  p2_23_(0.0), p1_23_(0.0),
43  p1_13_(0.0), p2_13_(0.0),
44  p1_Parent_(0.0), p2_Parent_(0.0), p3_Parent_(0.0),
45  squareDP_(calcSquareDPCoords), warnings_(kTRUE)
46 {
47  // Constructor
48  mass_.clear(); mMin_.clear(); mMax_.clear(); mSqMin_.clear(); mSqMax_.clear();
49  mSq_.clear();
50  mSqDiff_.clear();
51  mDiff_.clear();
52 
53  mass_.push_back(m1_);
54  mass_.push_back(m2_);
55  mass_.push_back(m3_);
56 
57  mSq_.push_back(m1Sq_);
58  mSq_.push_back(m2Sq_);
59  mSq_.push_back(m3Sq_);
60 
61  // DP max and min kinematic boundary for circumscribed box
62  // (see figure in PDG book).
63  for (Int_t i = 0; i < 3; i++) {
64  mMin_.push_back(mDTot_ - mass_[i]);
65  mMax_.push_back(mParent_ - mass_[i]);
66  mSqMin_.push_back(mMin_[i]*mMin_[i]);
67  mSqMax_.push_back(mMax_[i]*mMax_[i]);
68  mSqDiff_.push_back(mSqMax_[i] - mSqMin_[i]);
69  mDiff_.push_back(mMax_[i] - mMin_[i]);
70  }
71 
72  if (this->squareDP()) {
73  std::cout<<"INFO in LauKinematics::LauKinematics : squareDP = kTRUE"<<std::endl;
74  } else {
75  std::cout<<"INFO in LauKinematics::LauKinematics : squareDP = kFALSE"<<std::endl;
76  }
77 
78  // add covariant factor calculation
79 }
80 
82 {
83  // Destructor
84 }
85 
86 void LauKinematics::updateKinematics(const Double_t m13Sq, const Double_t m23Sq)
87 {
88  // Sets the internal private data members
89  // m13Sq and m23Sq, as well as m13 and m23.
90  // Also calculate m12Sq and m12, given mParent defined in the constructor.
91  // Lastly, calculate the helicity cosines.
92 
93  // Update the 3 mass-squares
94  this->updateMassSquares(m13Sq, m23Sq);
95 
96  // Now update the helicity cosines
97  this->calcHelicities();
98 
99  // Also calculate the m', theta' variables
100  if (squareDP_) {this->calcSqDPVars();}
101 }
102 
103 void LauKinematics::updateSqDPKinematics(const Double_t mPrime, const Double_t thetaPrime)
104 {
105  // From square DP co-ordinates, calculate remaining kinematic variables
106  this->updateSqDPMassSquares(mPrime, thetaPrime);
107 
108  // Finally calculate the helicities and track momenta
109  this->calcHelicities();
110 }
111 
113 {
114  // For given m_12 and cos(theta_12) values, calculate m' and theta' for the square Dalitz plot
115  Double_t value = (2.0*(m12_ - mMin_[2])/mDiff_[2]) - 1.0;
116  mPrime_ = LauConstants::invPi*TMath::ACos(value);
117  thetaPrime_ = LauConstants::invPi*TMath::ACos(c12_);
118  // Sometimes events are assigned exactly thetaPrime = 0.0 or 1.0
119  // which gives problems with efficiency and other histograms
120  if (thetaPrime_ == 0.0)
121  {
122  thetaPrime_ += 1.0e-10;
123  } else if (thetaPrime_ == 1.0)
124  {
125  thetaPrime_ -= 1.0e-10;
126  }
127 }
128 
130 {
131  return this->calcSqDPJacobian(mPrime_,thetaPrime_);
132 }
133 
134 Double_t LauKinematics::calcSqDPJacobian(const Double_t mPrime, const Double_t thPrime) const
135 {
136  // Calculate the Jacobian for the transformation
137  // m23^2, m13^2 -> m', theta' (square DP)
138  const Double_t m12 = 0.5*mDiff_[2]*(1.0 + TMath::Cos(LauConstants::pi*mPrime)) + mMin_[2];
139  const Double_t m12Sq = m12*m12;
140 
141  const Double_t e1Cms12 = (m12Sq - m2Sq_ + m1Sq_)/(2.0*m12);
142  const Double_t e3Cms12 = (mParentSq_ - m12Sq - m3Sq_)/(2.0*m12);
143 
144  const Double_t p1Cms12 = this->pCalc(e1Cms12, m1Sq_);
145  const Double_t p3Cms12 = this->pCalc(e3Cms12, m3Sq_);
146 
147  const Double_t deriv1 = LauConstants::piBy2*mDiff_[2]*TMath::Sin(LauConstants::pi*mPrime);
148  const Double_t deriv2 = LauConstants::pi*TMath::Sin(LauConstants::pi*thPrime);
149 
150  const Double_t jacobian = 4.0*p1Cms12*p3Cms12*m12*deriv1*deriv2;
151 
152  return jacobian;
153 }
154 
155 void LauKinematics::updateMassSquares(const Double_t m13Sq, const Double_t m23Sq)
156 {
157  m13Sq_ = m13Sq;
158  if (m13Sq_ > 0.0) {
159  m13_ = TMath::Sqrt(m13Sq_);
160  } else {
161  m13_ = 0.0;
162  }
163 
164  m23Sq_ = m23Sq;
165  if (m23Sq_ > 0.0) {
166  m23_ = TMath::Sqrt(m23Sq_);
167  } else {
168  m23_ = 0.0;
169  }
170 
171  // Now calculate m12Sq and m12.
172  this->calcm12Sq();
173 
174  // Calculate momenta of tracks in parent (B, D etc.) rest-frame
175  this->calcParentFrameMomenta();
176 }
177 
178 void LauKinematics::updateSqDPMassSquares(const Double_t mPrime, const Double_t thetaPrime)
179 {
180  // From square DP co-ordinates, calculate only the mass-squares
181 
182  // First set the square-DP variables
183  mPrime_ = mPrime; thetaPrime_ = thetaPrime;
184 
185  // Next calculate m12 and c12
186  Double_t m12 = 0.5*mDiff_[2]*(1.0 + TMath::Cos(LauConstants::pi*mPrime)) + mMin_[2];
187  Double_t c12 = TMath::Cos(LauConstants::pi*thetaPrime);
188 
189  // From these calculate m13Sq and m23Sq
190  this->updateMassSq_m12(m12, c12);
191 
192  // Calculate momenta of tracks in parent (B, D etc.) rest-frame
193  this->calcParentFrameMomenta();
194 }
195 
197 {
198  // Calculate m12Sq from m13Sq and m23Sq.
200 
201  // If m12Sq is too low, return lower limit,
202  // and modify m13Sq accordingly.
203  if (m12Sq_ < mSqMin_[2]) {
204  m12Sq_ = mSqMin_[2] + 1.0e-3;
206  }
207 
208  if (m12Sq_ > 0.0) {
209  m12_ = TMath::Sqrt(m12Sq_);
210  } else {
211  m12_ = 0.0;
212  }
213 }
214 
216 {
217  Double_t e1 = (mParentSq_ + m1Sq_ - m23Sq_) / (2.0*mParent_);
218  Double_t e2 = (mParentSq_ + m2Sq_ - m13Sq_) / (2.0*mParent_);
219  Double_t e3 = (mParentSq_ + m3Sq_ - m12Sq_) / (2.0*mParent_);
220 
221  p1_Parent_ = this->pCalc(e1, m1Sq_);
222  p2_Parent_ = this->pCalc(e2, m2Sq_);
223  p3_Parent_ = this->pCalc(e3, m3Sq_);
224 }
225 
227 {
228  // Calculate helicity angle cosines, given m12Sq, m13Sq and m23Sq.
229  // The notation is confusing for the helicity angle cosines:
230  // cij is helicity angle for the pair which is made from tracks i and j.
231  // It is the angle between tracks i and k in the ij rest frame
232  // Make sure that setMassSqVars is run first.
233  Int_t zero(0), one(1), two(2);
234 
235  c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
236  p1_12_ = qi_; p3_12_ = qk_; // i, j = 12 (rest frame), k = 3
237  c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
238  p2_23_ = qi_; p1_23_ = qk_; // i, j = 23 (rest frame), k = 1
239  c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
240  p1_13_ = qi_; p2_13_ = qk_; // i, j = 31 (rest frame), k = 2
241 
242 }
243 
244 Double_t LauKinematics::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
245 {
246  // Routine to calculate the cos(helicity) variables from the masses of the particles.
247  // cij is helicity angle for the pair which is made from tracks i and j.
248  // It is the angle between tracks i and k in the ij rest frame.
249  Double_t EiCmsij = (mijSq - mSq_[j] + mSq_[i])/(2.0*mij);
250  Double_t EkCmsij = (mParentSq_ - mijSq - mSq_[k])/(2.0*mij);
251 
252  if (EiCmsij < mass_[i]) {
253  if (warnings_) {
254  std::cerr<<"WARNING in LauKinematics::cFromM : EiCmsij = "<<EiCmsij<<" too small < mass_["<<i<<"] = "<<mass_[i]<<" in cFromM, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
255  std::cerr<<" : mijSq = "<<mijSq<<"; mij = "<<mij<<"; mSq_["<<j<<"] = "<<mSq_[j]<<"; mSq_["<<i<<"] = "<<mSq_[i]<<std::endl;
256  }
257  return 0.0;
258  }
259  if (EkCmsij < mass_[k]) {
260  if (warnings_) {
261  std::cerr<<"WARNING in LauKinematics::cFromM : EkCmsij = "<<EkCmsij<<" too small < mass_["<<k<<"] = "<<mass_[k]<<" in cFromM, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
262  std::cerr<<" : mijSq = "<<mijSq<<"; mij = "<<mij<<"; mSq_["<<j<<"] = "<<mSq_[j]<<"; mSq_["<<i<<"] = "<<mSq_[i]<<std::endl;
263  }
264  return 0.0;
265  }
266 
267  // Find track i and k momenta in ij rest frame
268  qi_ = this->pCalc(EiCmsij, mSq_[i]);
269  qk_ = this->pCalc(EkCmsij, mSq_[k]);
270 
271  // Find ij helicity angle
272  Double_t cosHel = -(mikSq - mSq_[i] - mSq_[k] - 2.0*EiCmsij*EkCmsij)/(2.0*qi_*qk_);
273 
274  if (cosHel > 1.0) {
275  cosHel = 1.0;
276  } else if (cosHel < -1.0) {
277  cosHel = -1.0;
278  }
279 
280  if (i == 1) {cosHel *= -1.0;}
281 
282  return cosHel;
283 }
284 
285 Double_t LauKinematics::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
286 {
287  // Returns the mass-squared for a pair of particles, i,j, given their
288  // invariant mass (squared) and the helicity angle.
289  // cij is helicity angle for the pair which is made from tracks i and j.
290  // It is the angle between tracks i and k in the ij rest frame.
291 
292  Double_t cosHel( cij );
293  if (i == 1) {cosHel *= -1.0;}
294 
295  Double_t EiCmsij = (mijSq - mSq_[j] + mSq_[i])/(2.0*mij);
296  Double_t EkCmsij = (mParentSq_ - mijSq - mSq_[k])/(2.0*mij);
297 
298  if (TMath::Abs(EiCmsij - mass_[i]) > 1e-6 && EiCmsij < mass_[i]) {
299  if (warnings_) {
300  std::cerr<<"WARNING in LauKinematics::mFromC : EiCmsij = "<<EiCmsij<<" < "<<mass_[i]<<" in mFromC, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
301  }
302  return 0.0;
303  }
304  if (TMath::Abs(EkCmsij - mass_[k]) > 1e-6 && EkCmsij < mass_[k]) {
305  if (warnings_) {
306  std::cerr<<"WARNING in LauKinematics::mFromC : EkCmsij = "<<EkCmsij<<" < "<<mass_[k]<<" in mFromC, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
307  }
308  return 0.0;
309  }
310 
311  // Find track i and k momenta in ij rest fram
312  Double_t qi = this->pCalc(EiCmsij, mSq_[i]);
313  Double_t qk = this->pCalc(EkCmsij, mSq_[k]);
314 
315  // Find mikSq
316  Double_t massSq = mSq_[i] + mSq_[k] + 2.0*EiCmsij*EkCmsij - 2.0*qi*qk*cosHel;
317 
318  if (massSq < mSqMin_[j]) {
319  if (warnings_) {
320  std::cerr<<"WARNING in LauKinematics::mFromC : mFromC below bound: i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
321  }
322  massSq = mSqMin_[j];
323  }
324 
325  return massSq;
326 }
327 
328 void LauKinematics::genFlatPhaseSpace(Double_t& m13Sq, Double_t& m23Sq) const
329 {
330  // Routine to generate flat phase-space events.
331  // DP max kinematic boundaries in circumscribed box
332  // See DP figure in PDG book.
333  // m13max=mbrec-mass(2)
334  // m13sqmax=m13max*m13max
335  // m23max=mbrec-mass(1)
336  // m23sqmax=m23max*m23max
337 
338  // Generate m13Sq and m23Sq flat within DP overall rectangular box
339  // Loop if the generated point is not within the DP kinematic boundary
340 
341  do {
342  m13Sq = mSqMin_[1] + LauRandom::randomFun()->Rndm()*mSqDiff_[1];
343  m23Sq = mSqMin_[0] + LauRandom::randomFun()->Rndm()*mSqDiff_[0];
344 
345  } while ( ! this->withinDPLimits( m13Sq, m23Sq ) );
346 }
347 
348 void LauKinematics::genFlatSqDP(Double_t& mPrime, Double_t& thetaPrime) const
349 {
350  // Generate random event in the square Dalitz plot
351  mPrime = LauRandom::randomFun()->Rndm();
352  thetaPrime = LauRandom::randomFun()->Rndm();
353 }
354 
355 Bool_t LauKinematics::withinDPLimits(const Double_t m13Sq, const Double_t m23Sq) const
356 {
357  // Find out whether the point (m13Sq,m23Sq) is within the limits of the
358  // Dalitz plot. The limits are specified by the invariant masses
359  // of the parent (e.g. B) and its three daughters that were
360  // defined in the constructor of this class. Here
361  // m_13Sq = square of invariant mass of daughters 1 and 3
362  // m_23Sq = square of invariant mass of daughters 2 and 3.
363 
364  Bool_t withinDP = kFALSE;
365 
366  // First check that m13Sq is within its absolute min and max
367  if (!((m13Sq > mSqMin_[1]) && (m13Sq < mSqMax_[1]))) {
368  return kFALSE;
369  }
370 
371  // Now for the given value of m13Sq calculate the local min and max of m23Sq
372  Double_t m13 = TMath::Sqrt(m13Sq);
373 
374  Double_t e3Cms13 = (m13Sq - m1Sq_ + m3Sq_)/(2.0*m13);
375  Double_t p3Cms13 = this->pCalc(e3Cms13, m3Sq_);
376 
377  Double_t e2Cms13 = (mParentSq_ - m13Sq - m2Sq_)/(2.0*m13);
378  Double_t p2Cms13 = this->pCalc(e2Cms13, m2Sq_);
379 
380  Double_t term = 2.0*e2Cms13*e3Cms13 + m2Sq_ + m3Sq_;
381 
382  Double_t m23SqLocMin = term - 2.0*p2Cms13*p3Cms13;
383  Double_t m23SqLocMax = term + 2.0*p2Cms13*p3Cms13;
384 
385  // Check whether the given value of m23Sq satisfies these bounds
386  if (m23Sq > m23SqLocMin && m23Sq < m23SqLocMax) {
387  withinDP = kTRUE;
388  }
389 
390  return withinDP;
391 }
392 
393 Bool_t LauKinematics::withinDPLimits2(const Double_t m13Sq, const Double_t m23Sq) const
394 {
395  // Same as withinDPLimits, but this time testing whether the m13Sq
396  // variable is within the kinematic boundary for the given m23Sq value
397 
398  Bool_t withinDP = kFALSE;
399 
400  // First check that m13Sq is within its absolute min and max
401  if (!((m23Sq > mSqMin_[0]) && (m23Sq < mSqMax_[0]))) {
402  return kFALSE;
403  }
404 
405  // Now for the given value of m13Sq calculate the local min and max of m23Sq
406  Double_t m23 = TMath::Sqrt(m23Sq);
407 
408  Double_t e3Cms23 = (m23Sq - m2Sq_ + m3Sq_)/(2.0*m23);
409  Double_t p3Cms23 = this->pCalc(e3Cms23, m3Sq_);
410 
411  Double_t e1Cms23 = (mParentSq_ - m23Sq - m1Sq_)/(2.0*m23);
412  Double_t p1Cms23 = this->pCalc(e1Cms23, m1Sq_);
413 
414  Double_t term = 2.0*e3Cms23*e1Cms23 + m1Sq_ + m3Sq_;
415 
416  Double_t m13SqLocMin = term - 2.0*p3Cms23*p1Cms23;
417  Double_t m13SqLocMax = term + 2.0*p3Cms23*p1Cms23;
418 
419  // Check whether the given value of m23Sq satisfies these bounds
420  if (m13Sq > m13SqLocMin && m13Sq < m13SqLocMax) {
421  withinDP = kTRUE;
422  }
423 
424  return withinDP;
425 }
426 
427 Bool_t LauKinematics::withinSqDPLimits(const Double_t mPrime, const Double_t thetaPrime) const
428 {
429  // Check whether m' and theta' are between 0 and 1
430  Bool_t withinDP(kFALSE);
431  if (mPrime > 0.0 && mPrime < 1.0 &&
432  thetaPrime > 0.0 && thetaPrime < 1.0) {
433  withinDP = kTRUE;
434  }
435 
436  return withinDP;
437 }
438 
439 Double_t LauKinematics::calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const
440 {
441  // Calculate one massSq from the other two
442  return mParentSq_ + mSqDTot_ - firstMassSq - secondMassSq;
443 }
444 
446 {
447  return this->distanceFromDPCentre(m13Sq_,m23Sq_);
448 }
449 
450 Double_t LauKinematics::distanceFromDPCentre(const Double_t m13Sq, const Double_t m23Sq) const
451 {
452  // DP centre is defined as the point where m12=m13=m23
453 
454  // First find the m^2_ij value for the centre
455  Double_t centreMijSq = (mParentSq_ + m1Sq_ + m2Sq_ + m3Sq_)/3.0;
456 
457  // Then find the difference between this and the two provided co-ords
458  Double_t diff13 = m13Sq - centreMijSq;
459  Double_t diff23 = m23Sq - centreMijSq;
460 
461  // Calculate the total distance
462  Double_t distance = TMath::Sqrt(diff13*diff13 + diff23*diff23);
463  return distance;
464 }
465 
466 Double_t LauKinematics::pCalc(const Double_t energy, const Double_t massSq) const
467 {
468  // Routine to calculate the momentum of a particle, given its energy and
469  // invariant mass (squared).
470  Double_t arg = energy*energy - massSq;
471 
472  if (arg < 0.0) {
473  //if (warnings_) {
474  //std::cerr<<"WARNING in LauKinematics::pCalc : arg < 0.0: "<<arg<<" for e = "<<energy<<", mSq = "<<massSq<<std::endl;
475  //}
476  arg = 0.0;
477  }
478  Double_t pCalcVal = TMath::Sqrt(arg);
479  return pCalcVal;
480 
481 }
482 
484 {
485  // Flips the DP variables m13^2 <-> m23^2.
486  // Used in the case of symmetrical Dalitz plots (i.e. when two of the
487  // daughter tracks are the same type) within the
488  // LauIsobarDynamics::resAmp function.
490 }
491 
493 {
494  // Cyclically rotates the DP variables (m12 -> m23, m23 -> m13, m13 -> m12)
495  // Used in the case of fully symmetric Dalitz plots (i.e. when all
496  // three of the daughter tracks are the same type) within the
497  // LauIsobarDynamics::resAmp function.
499 }
500 
501 void LauKinematics::updateMassSq_m23(const Double_t m23, const Double_t c23)
502 {
503  // Update the variables m12Sq_ and m13Sq_ given m23 and c23.
504  m23_ = m23; m23Sq_ = m23*m23; c23_ = c23;
505 
506  const Int_t zero(0), one(1), two(2);
507  m12Sq_ = this->mFromC(m23Sq_, c23_, m23_, one, two, zero);
509 
510  m12_ = TMath::Sqrt(m12Sq_);
511  m13_ = TMath::Sqrt(m13Sq_);
512 }
513 
514 void LauKinematics::updateMassSq_m13(const Double_t m13, const Double_t c13)
515 {
516  // Update the variables m12Sq_ and m23Sq_ given m13 and c13.
517  m13_ = m13; m13Sq_ = m13*m13; c13_ = c13;
518 
519  const Int_t zero(0), one(1), two(2);
520  m23Sq_ = this->mFromC(m13Sq_, c13_, m13_, two, zero, one);
522 
523  m23_ = TMath::Sqrt(m23Sq_);
524  m12_ = TMath::Sqrt(m12Sq_);
525 }
526 
527 void LauKinematics::updateMassSq_m12(const Double_t m12, const Double_t c12)
528 {
529  // Update the variables m23Sq_ and m13Sq_ given m12 and c12.
530  m12_ = m12; m12Sq_ = m12*m12; c12_ = c12;
531 
532  const Int_t zero(0), one(1), two(2);
533  m13Sq_ = this->mFromC(m12Sq_, c12_, m12_, zero, one, two);
535  m13_ = TMath::Sqrt(m13Sq_);
536  m23_ = TMath::Sqrt(m23Sq_);
537 }
538 
539 void LauKinematics::updateKinematicsFrom23(const Double_t m23, const Double_t c23)
540 {
541  // Calculate the other mass squares
542  this->updateMassSq_m23(m23,c23);
543 
544  // Calculate the remaining helicity angles
545  this->calcHelicities();
546 
547  // Calculate momenta of tracks in parent (B, D etc.) rest-frame
548  this->calcParentFrameMomenta();
549 
550  // Also calculate the m', theta' variables
551  if (squareDP_) {this->calcSqDPVars();}
552 }
553 
554 void LauKinematics::updateKinematicsFrom13(const Double_t m13, const Double_t c13)
555 {
556  // Calculate the other mass squares
557  this->updateMassSq_m13(m13,c13);
558 
559  // Calculate the remaining helicity angles
560  this->calcHelicities();
561 
562  // Calculate momenta of tracks in parent (B, D etc.) rest-frame
563  this->calcParentFrameMomenta();
564 
565  // Also calculate the m', theta' variables
566  if (squareDP_) {this->calcSqDPVars();}
567 }
568 
569 void LauKinematics::updateKinematicsFrom12(const Double_t m12, const Double_t c12)
570 {
571  // Calculate the other mass squares
572  this->updateMassSq_m12(m12,c12);
573 
574  // Calculate the remaining helicity angles
575  this->calcHelicities();
576 
577  // Calculate momenta of tracks in parent (B, D etc.) rest-frame
578  this->calcParentFrameMomenta();
579 
580  // Also calculate the m', theta' variables
581  if (squareDP_) {this->calcSqDPVars();}
582 }
583 
584 Double_t LauKinematics::genm13Sq() const
585 {
586  Double_t m13Sq = mSqMin_[1] + LauRandom::randomFun()->Rndm()*mSqDiff_[1];
587  return m13Sq;
588 }
589 
590 Double_t LauKinematics::genm23Sq() const
591 {
592  Double_t m23Sq = mSqMin_[0] + LauRandom::randomFun()->Rndm()*mSqDiff_[0];
593  return m23Sq;
594 }
595 
596 Double_t LauKinematics::genm12Sq() const
597 {
598  Double_t m12Sq = mSqMin_[2] + LauRandom::randomFun()->Rndm()*mSqDiff_[2];
599  return m12Sq;
600 }
601 
602 
603 void LauKinematics::drawDPContour(Int_t orientation, Int_t nbins)
604 {
605  // orientation -
606  // 1323 : x = m13, y = m23
607  // etc.
608 
609  Double_t m1 = this->getm1();
610  Double_t m2 = this->getm2();
611  Double_t m3 = this->getm3();
612  Double_t mParent = this->getmParent();
613 
614  Double_t m13SqMin = this->getm13SqMin();
615  Double_t m23SqMin = this->getm23SqMin();
616  Double_t m12SqMin = this->getm12SqMin();
617  Double_t m13SqMax = this->getm13SqMax();
618  Double_t m23SqMax = this->getm23SqMax();
619  Double_t m12SqMax = this->getm12SqMax();
620 
621  Double_t xMin(0.0);
622  Double_t xMax(mParent*mParent);
623  Double_t yMin(0.0);
624  Double_t yMax(mParent*mParent);
625  if (orientation == 1323) {
626  xMin = m13SqMin-1.0; xMax = m13SqMax+1.0;
627  yMin = m23SqMin-1.0; yMax = m23SqMax+1.0;
628  } else if (orientation == 2313) {
629  xMin = m23SqMin-1.0; xMax = m23SqMax+1.0;
630  yMin = m13SqMin-1.0; yMax = m13SqMax+1.0;
631  } else if (orientation == 1213) {
632  xMin = m12SqMin-1.0; xMax = m12SqMax+1.0;
633  yMin = m13SqMin-1.0; yMax = m13SqMax+1.0;
634  } else if (orientation == 1312) {
635  xMin = m13SqMin-1.0; xMax = m13SqMax+1.0;
636  yMin = m12SqMin-1.0; yMax = m12SqMax+1.0;
637  } else if (orientation == 1223) {
638  xMin = m12SqMin-1.0; xMax = m12SqMax+1.0;
639  yMin = m23SqMin-1.0; yMax = m23SqMax+1.0;
640  } else if (orientation == 2312) {
641  xMin = m23SqMin-1.0; xMax = m23SqMax+1.0;
642  yMin = m12SqMin-1.0; yMax = m12SqMax+1.0;
643  } else {
644  std::cerr<<"ERROR in LauKinematics::drawDPContour : Unrecognised orientation, not drawing contour."<<std::endl;
645  return;
646  }
647 
648  Int_t npar = 4;
649  TF2 * f2 = new TF2("contour", &dal, xMin, xMax, yMin, yMax, npar);
650 
651  // Set the parameters
652  f2->SetParameter(0,mParent);
653  if (orientation == 1323) {
654  f2->SetParameter(1,m1);
655  f2->SetParameter(2,m2);
656  f2->SetParameter(3,m3);
657  } else if (orientation == 2313) {
658  f2->SetParameter(1,m2);
659  f2->SetParameter(2,m1);
660  f2->SetParameter(3,m3);
661  } else if (orientation == 1213) {
662  f2->SetParameter(1,m2);
663  f2->SetParameter(2,m3);
664  f2->SetParameter(3,m1);
665  } else if (orientation == 1312) {
666  f2->SetParameter(1,m3);
667  f2->SetParameter(2,m2);
668  f2->SetParameter(3,m1);
669  } else if (orientation == 1223) {
670  f2->SetParameter(1,m1);
671  f2->SetParameter(2,m3);
672  f2->SetParameter(3,m2);
673  } else if (orientation == 2312) {
674  f2->SetParameter(1,m3);
675  f2->SetParameter(2,m1);
676  f2->SetParameter(3,m2);
677  }
678 
679  // Set up the contour to be drawn when the value of the function == 1.0
680  Double_t b[]={1.0};
681  f2->SetContour(1,b);
682 
683  // Set the number of bins for the contour to be sampled over
684  f2->SetNpx(nbins);
685  f2->SetNpy(nbins);
686  // and the line style
687  f2->SetLineWidth(3);
688  f2->SetLineStyle(kSolid);
689 
690  // Draw the contour on top of the histo that should already have been drawn
691  f2->DrawCopy("same");
692 
693  delete f2;
694 }
695 
696 Double_t dal(Double_t* x, Double_t* par)
697 {
698  Double_t mParent = par[0];
699  Double_t mi = par[1];
700  Double_t mj = par[2];
701  Double_t mk = par[3];
702 
703  Double_t mikSq=x[0];
704  Double_t mjkSq=x[1];
705  Double_t mik = TMath::Sqrt(mikSq);
706  Double_t mjk = TMath::Sqrt(mjkSq);
707 
708  Double_t ejcmsik = (mParent*mParent-mj*mj-mik*mik)/(2.0*mik);
709  Double_t ekcmsik = (mik*mik+mk*mk-mi*mi)/(2.0*mik);
710  if (ekcmsik<mk || ejcmsik<mj) return 2.0;
711 
712  Double_t pj = TMath::Sqrt(ejcmsik*ejcmsik-mj*mj);
713  Double_t pk = TMath::Sqrt(ekcmsik*ekcmsik-mk*mk);
714  Double_t coshelik = (mjk*mjk - mk*mk - mj*mj - 2.0*ejcmsik*ekcmsik)/(2.0*pj*pk);
715 
716  Double_t coshelikSq = coshelik*coshelik;
717  return coshelikSq;
718 }
719 
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 getm13SqMin() const
Get the m13Sq minimum, (m1 + m3)^2.
const Double_t mParentSq_
Mass of parent particle squared.
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
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.
std::vector< Double_t > mSqMin_
Vector of the minimum mijSq values.
ClassImp(LauAbsCoeffSet)
Double_t m23_
Invariant mass m23.
Double_t p2_Parent_
Momentum of track 2 in parent rest frame.
std::vector< Double_t > mass_
Vector of daughter particles masses.
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 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.
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 getm13SqMax() const
Get the m13Sq maximum, (mParent - m2)^2.
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 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.
File containing declaration of LauKinematics class.
Double_t genm23Sq() const
Randomly generate the invariant mass squared m23Sq.
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.
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) ...
const Double_t m3Sq_
Mass of particle 3 squared.
void flipAndUpdateKinematics()
Flips the DP variables m13^2 &lt;-&gt; m23^2 and recalculates all kinematic quantities. ...
const Double_t invPi
One over Pi.
const Double_t pi
Pi.
Definition: LauConstants.hh:89
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.
File containing LauRandom namespace.
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...
Double_t genm13Sq() const
Randomly generate the invariant mass squared m13Sq.
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 ...
Double_t dal(Double_t *x, Double_t *par)
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.
std::vector< Double_t > mDiff_
Vector of the difference between the mMax and mMin.
Bool_t squareDP_
Should we calculate the square DP co-ordinates or not?
File containing LauConstants namespace.
Double_t m12_
Invariant mass m12.
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 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.
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
Double_t genm12Sq() const
Randomly generate the invariant mass squared m12Sq.
const Double_t piBy2
Pi divided by two.
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...
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.