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