meerkat is hosted by Hepforge, IPPP Durham
Meerkat  v1r3
Multidimensional kernel density estimation package
ExtendedDalitzPhaseSpace.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <vector>
3 #include <stdlib.h>
4 
5 #include "TMath.h"
6 
8 
9 #include "Logger.hh"
10 
12  Double_t mdmin,
13  Double_t mdmax,
14  Double_t ma,
15  Double_t mb,
16  Double_t mc) : AbsPhaseSpace(phspName) {
17 
18  m_a = ma;
19  m_b = mb;
20  m_c = mc;
21  m_d_min = mdmin;
22  m_d_max = mdmax;
23 
24  if (mdmin >= mdmax) {
25  Logger::print(2,"%20.20s ERROR: Minimum mother mass is larger than maximum\n", m_name);
26  abort();
27  }
28 
29  Logger::print(0,"%20.20s INFO: Creating extended Dalitz phase space for mother mass in range (%f,%f) and daughters of masses %f, %f, %f\n",
31 
32  m_a2 = m_a*m_a;
33  m_b2 = m_b*m_b;
34  m_c2 = m_c*m_c;
35  m_SqSum = m_a2 + m_b2 + m_c2;
36 
37  m_MinAB = TMath::Power(m_a + m_b, 2);
38  m_MaxAB = TMath::Power(m_d_max - m_c, 2);
39  m_MinBC = TMath::Power(m_b + m_c, 2);
40  m_MaxBC = TMath::Power(m_d_max - m_a, 2);
41 
42 }
43 
45 
46 }
47 
48 Bool_t ExtendedDalitzPhaseSpace::withinLimits(std::vector<Double_t> &x) {
49 
50  Double_t m_d = x[0];
51  Double_t m2ab = x[1];
52  Double_t m2bc = x[2];
53 
54  if (m_d < m_d_min || m_d > m_d_max) return 0;
55 
56  Double_t m_d2 = m_d*m_d;
57 // Double_t sqSum = m_SqSum + m_d2;
58 
59  Double_t minAB = TMath::Power(m_a + m_b, 2);
60  Double_t maxAB = TMath::Power(m_d - m_c, 2);
61  Double_t minBC = TMath::Power(m_b + m_c, 2);
62  Double_t maxBC = TMath::Power(m_d - m_a, 2);
63 
64 // Logger::print(0,"%f %f %f %f %f %f\n", m2ab, m2bc, m_MinAB, m_MinBC, m_MaxAB, m_MaxBC);
65 
66  if (m2ab < minAB || m2ab > maxAB || m2bc < minBC || m2bc > maxBC) return 0;
67 
68 // Logger::print(0,"1\n");
69 
70  Double_t eb = (m2ab - m_a2 + m_b2)/2./TMath::Sqrt(m2ab);
71  Double_t ec = (m_d2 - m2ab - m_c2)/2./TMath::Sqrt(m2ab);
72 
73  Double_t p2b = TMath::Power(eb, 2) - m_b2;
74  Double_t p2c = TMath::Power(ec, 2) - m_c2;
75 
76  if (p2c<0 || p2b<0) return 0;
77 
78 // Logger::print(0,"2\n");
79 
80  Double_t m2bc_max = TMath::Power(eb+ec, 2) - TMath::Power(sqrt(p2b) - sqrt(p2c), 2);
81  Double_t m2bc_min = TMath::Power(eb+ec, 2) - TMath::Power(sqrt(p2b) + sqrt(p2c), 2);
82 
83  if (m2bc < m2bc_min || m2bc > m2bc_max) return 0;
84 
85 // Logger::print(0,"3\n");
86 
87  return 1;
88 }
89 
91  switch(var) {
92  case 0: return m_d_min;
93  case 1: return m_MinAB;
94  case 2: return m_MinBC;
95  default:
96  Logger::print(2,"%20.20s ERROR: var=%d for lower limit of 2D phase space\n", m_name, var);
97  abort();
98  }
99 }
100 
102  switch(var) {
103  case 0: return m_d_max;
104  case 1: return m_MaxAB;
105  case 2: return m_MaxBC;
106  default:
107  Logger::print(2,"%20.20s ERROR: var=%d for upper limit of 2D phase space\n", m_name, var);
108  abort();
109  }
110 }
111 
112 Bool_t ExtendedDalitzPhaseSpace::limits(UInt_t var, std::vector<Double_t> &x,
113  Double_t* lower, Double_t* upper) {
114 
115  if (var == 0) {
116  *upper = m_d_max;
117  *lower = m_d_min;
118  }
119  else if (var == 1) { // AB
120  Double_t m_d = x[0];
121  Double_t m_d2 = m_d*m_d;
122  Double_t m2bc = x[2];
123 
124 // Double_t minAB = TMath::Power(m_a + m_b, 2);
125 // Double_t maxAB = TMath::Power(m_d - m_c, 2);
126  Double_t minBC = TMath::Power(m_b + m_c, 2);
127  Double_t maxBC = TMath::Power(m_d - m_a, 2);
128 
129  if (m2bc < minBC || m2bc > maxBC) return 0;
130 
131  Double_t eb = (m2bc + m_b2 - m_c2)/2./TMath::Sqrt(m2bc);
132  Double_t ea = (m_d2 - m2bc - m_a2)/2./TMath::Sqrt(m2bc);
133 
134  Double_t p2b = TMath::Power(eb, 2) - m_b2;
135  Double_t p2a = TMath::Power(ea, 2) - m_a2;
136 
137  if (p2a<0 || p2b<0) return 0;
138 
139 // Logger::print(0,"2\n");
140 
141  *upper = TMath::Power(ea+eb, 2) - TMath::Power(sqrt(p2a) - sqrt(p2b), 2);
142  *lower = TMath::Power(ea+eb, 2) - TMath::Power(sqrt(p2a) + sqrt(p2b), 2);
143 
144  }
145  else if (var == 2) { // BC
146 
147  Double_t m_d = x[0];
148  Double_t m_d2 = m_d*m_d;
149  Double_t m2ab = x[1];
150 
151  Double_t minAB = TMath::Power(m_a + m_b, 2);
152  Double_t maxAB = TMath::Power(m_d - m_c, 2);
153 // Double_t minBC = TMath::Power(m_b + m_c, 2);
154 // Double_t maxBC = TMath::Power(m_d - m_a, 2);
155 
156  if (m2ab < minAB || m2ab > maxAB) return 0;
157 
158  Double_t eb = (m2ab - m_a2 + m_b2)/2./TMath::Sqrt(m2ab);
159  Double_t ec = (m_d2 - m2ab - m_c2)/2./TMath::Sqrt(m2ab);
160 
161  Double_t p2b = TMath::Power(eb, 2) - m_b2;
162  Double_t p2c = TMath::Power(ec, 2) - m_c2;
163 
164  if (p2c<0 || p2b<0) return 0;
165 
166 // Logger::print(0,"2\n");
167 
168  *upper = TMath::Power(eb+ec, 2) - TMath::Power(sqrt(p2b) - sqrt(p2c), 2);
169  *lower = TMath::Power(eb+ec, 2) - TMath::Power(sqrt(p2b) + sqrt(p2c), 2);
170 
171  } else {
172  Logger::print(2,"%20.20s ERROR: var=%d for limits of 2D phase space\n", m_name, var);
173  abort();
174  }
175 
176  return 1;
177 
178 }
Double_t m_a
Mass of the particle A.
char m_name[256]
Phase space name.
Double_t m_MinAB
Lower limit of AB invariant mass values.
Bool_t limits(UInt_t var, std::vector< Double_t > &x, Double_t *lower, Double_t *upper)
Return limits (lower and upper) for the variable at the certain point of the phase space...
Double_t m_b
Mass of the particle B.
void print(int level, const char *format,...)
Definition: Logger.cpp:27
Double_t m_a2
Squared mass of the particle A.
Double_t m_b2
Squared mass of the particle B.
virtual ~ExtendedDalitzPhaseSpace()
Destructor.
Double_t m_SqSum
Sum of squared masses of daughter particles.
Double_t m_MinBC
Lower limit of BC invariant mass values.
Abstract class which defines phase space interface.
Double_t m_MaxAB
Upper limit of AB invariant mass values.
Double_t m_MaxBC
Upper limit of BC invariant mass values.
Double_t m_d_max
Maximum mass of the mother particle.
ExtendedDalitzPhaseSpace(const char *phaseSpaceName, Double_t mDmin, Double_t mDmax, Double_t mA, Double_t mB, Double_t mC)
Constructor.
Double_t m_c
Mass of the particle C.
Bool_t withinLimits(std::vector< Double_t > &x)
Check if the point is within the phase space limits.
Double_t m_c2
Squared mass of the particle C.
Double_t lowerLimit(UInt_t var)
Get lower limit.
Double_t upperLimit(UInt_t var)
Get upper limit.
Double_t m_d_min
Minimum mass of the mother particle.