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