[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtDalitzPlot.cpp

//-----------------------------------------------------------------------
// File and Version Information: 
//      $Id: EvtDalitzPlot.cpp,v 1.3 2009-03-16 15:53:27 robbep Exp $
// 
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations. If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information:
//      Copyright (C) 1998 Caltech, UCSB
//
// Module creator:
//      Alexei Dvoretskii, Caltech, 2001-2002.
//-----------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"

// Global 3-body Dalitz decay kinematics as defined by the mass
// of the mother and the daughters. Spins are not considered.

#include <math.h>
#include <assert.h>
#include <stdio.h>
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtDalitzPlot.hh"
#include "EvtGenBase/EvtTwoBodyVertex.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtDecayMode.hh"

using namespace EvtCyclic3;


EvtDalitzPlot::EvtDalitzPlot()
  : _mA(0.), _mB(0.), _mC(0.), _bigM(0.),
    _ldel(0.), _rdel(0.)
{}


EvtDalitzPlot::EvtDalitzPlot(double mA, double mB, double mC, double bigM,
			     double ldel, double rdel)
  : _mA(mA), _mB(mB), _mC(mC), _bigM(bigM),
    _ldel(ldel), _rdel(rdel) 
{
  sanityCheck();
}

EvtDalitzPlot::EvtDalitzPlot(const EvtDecayMode& mode, double ldel, double rdel )
{
  _mA = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(A)));
  _mB = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(B)));
  _mC = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(C)));
  _bigM = EvtPDL::getMeanMass(EvtPDL::getId(mode.mother()));

  _ldel = ldel;
  _rdel = rdel;
  
  sanityCheck();
}


EvtDalitzPlot::EvtDalitzPlot(const EvtDalitzPlot& other) 
  : _mA(other._mA), _mB(other._mB), _mC(other._mC), _bigM(other._bigM),
    _ldel(other._ldel), _rdel(other._rdel)
{}


EvtDalitzPlot::~EvtDalitzPlot()
{}


bool EvtDalitzPlot::operator==(const EvtDalitzPlot& other) const
{
  bool ret = false;
  if(_mA == other._mA && 
     _mB == other._mB &&
     _mC == other._mC &&
     _bigM == other._bigM) ret = true;

  return ret;
}


const EvtDalitzPlot* EvtDalitzPlot::clone() const
{
  return new EvtDalitzPlot(*this);
}


void EvtDalitzPlot::sanityCheck() const
{
  if(_mA < 0 || _mB < 0 || _mC < 0 || _bigM <= 0 || _bigM - _mA - _mB - _mC < 0.) {
    
    printf("Invalid Dalitz plot %f %f %f %f\n",_mA,_mB,_mC,_bigM);
    assert(0);
  }
  assert(_ldel <= 0.);
  assert(_rdel >= 0.);
}


double EvtDalitzPlot::m(Index i) const {

  double m = _mA;
  if(i == B) m = _mB;
  else
    if(i == C) m = _mC;

  return m;
}


double EvtDalitzPlot::sum() const
{
  return _mA*_mA + _mB*_mB + _mC*_mC + _bigM*_bigM;
}


double EvtDalitzPlot::qAbsMin(Pair i) const
{
  Index j = first(i);
  Index k = second(i);

  return (m(j) + m(k))*(m(j) + m(k));
}


double EvtDalitzPlot::qAbsMax(Pair i) const
{
  Index j = other(i);
  return (_bigM-m(j))*(_bigM-m(j));
}


double EvtDalitzPlot::qResAbsMin(EvtCyclic3::Pair i) const
{
  return qAbsMin(i) - sum()/3.;
}

double EvtDalitzPlot::qResAbsMax(EvtCyclic3::Pair i) const
{
  return  qAbsMax(i) - sum()/3.;
}

double EvtDalitzPlot::qHelAbsMin(EvtCyclic3::Pair i) const
{
  Pair j = next(i);
  Pair k = prev(i);
  return  (qAbsMin(j) - qAbsMax(k))/2.;
}

double EvtDalitzPlot::qHelAbsMax(EvtCyclic3::Pair i) const
{
  Pair j = next(i);
  Pair k = prev(i);
  return  (qAbsMax(j) - qAbsMin(k))/2.;
}


double EvtDalitzPlot::mAbsMin(Pair i) const
{
  return sqrt(qAbsMin(i));
}


double EvtDalitzPlot::mAbsMax(Pair i) const
{
  return sqrt(qAbsMax(i));
}


// parallel

double EvtDalitzPlot::qMin(Pair i, Pair j, double q) const
{
  if(i == j) return q;

  else {

    // Particle pair j defines the rest-frame
    // 0 - particle common to r.f. and angle calculations
    // 1 - particle belonging to r.f. but not angle
    // 2 - particle not belonging to r.f.

    Index k0 = common(i,j);
    Index k2 = other(j);
    Index k1 = other(k0,k2);

    // Energy, momentum of particle common to rest-frame and angle
    EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q)); 
    double pk = jpair.p();
    double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);


    // Energy and momentum of the other particle
    EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
    double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
    double pj = mother.p(EvtTwoBodyKine::A);


    // See PDG 34.4.3.1
    return (ek+ej)*(ek+ej) - (pk+pj)*(pk+pj);
  }
}


// antiparallel

double EvtDalitzPlot::qMax(Pair i, Pair j, double q) const
{

  if(i == j) return q;
  else {

    // Particle pair j defines the rest-frame
    // 0 - particle common to r.f. and angle calculations
    // 1 - particle belonging to r.f. but not angle
    // 2 - particle not belonging to r.f.
    
    Index k0 = common(i,j);
    Index k2 = other(j);
    Index k1 = other(k0,k2); 

    // Energy, momentum of particle common to rest-frame and angle
    EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q)); 
    double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
    double pk = jpair.p();

    // Energy and momentum of the other particle
    EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
    double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
    double pj = mother.p(EvtTwoBodyKine::A);

    
    // See PDG 34.4.3.1
    return (ek+ej)*(ek+ej) - (pk-pj)*(pk-pj);
  }
}


double EvtDalitzPlot::getArea(int N, Pair i, Pair j) const
{
  // Trapezoidal integral over qi. qj can be calculated.
  // The first and the last point are zero, so they are not counted

  double dh = (qAbsMax(i) - qAbsMin(i))/((double) N);
  double sum = 0;

  int ii;
  for(ii=1;ii<N;ii++) {

    double x = qAbsMin(i) + ii*dh;
    double dy = qMax(j,i,x) - qMin(j,i,x);
    sum += dy;
  }

  return sum * dh;
}


double EvtDalitzPlot::cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
{
  if(i1 == i2) return 1.;
  
  double qmax = qMax(i1,i2,q2);
  double qmin = qMin(i1,i2,q2);
  
  double cos = (qmax + qmin - 2*q1)/(qmax - qmin);
  
  return cos;
}


double EvtDalitzPlot::e(Index i, Pair j, double q) const
{
  if(i == other(j)) {
 
    // i does not belong to pair j

    return (bigM()*bigM()-q-m(i)*m(i))/2/sqrt(q);
  }
  else {
    
    // i and k make pair j

    Index k;
    if(first(j) == i) k = second(j);
    else k = first(j); 

    double e = (q + m(i)*m(i) - m(k)*m(k))/2/sqrt(q);	       
    return e;
  }
}


double EvtDalitzPlot::p(Index i, Pair j, double q) const
{
  double en = e(i,j,q);
  double p2 = en*en - m(i)*m(i);
  
  if(p2 < 0) {
    printf("Bad value of p2 %f %d %d %f %f\n",p2,i,j,en,m(i));
    assert(0);
  }

  return sqrt(p2);  
}


double EvtDalitzPlot::q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const
{
  if(i1 == i2) return q2;

  EvtCyclic3::Index f = first(i1);
  EvtCyclic3::Index s = second(i1);
  return m(f)*m(f) + m(s)*m(s) + 2*e(f,i2,q2)*e(s,i2,q2) - 2*p(f,i2,q2)*p(s,i2,q2)*cosTh;
}


double EvtDalitzPlot::jacobian(EvtCyclic3::Pair i, double q) const
{
  return 2*p(first(i),i,q)*p(other(i),i,q);  // J(BC) = 2pA*pB = 2pA*pC
}


EvtTwoBodyVertex EvtDalitzPlot::vD(Pair iRes, double m0, int L) const
{
  return EvtTwoBodyVertex(m(first(iRes)),
			  m(second(iRes)),m0,L);
}


EvtTwoBodyVertex EvtDalitzPlot::vB(Pair iRes, double m0, int L) const
{
  return EvtTwoBodyVertex(m0,m(other(iRes)),bigM(),L);
}


void EvtDalitzPlot::print() const
{
  // masses
  printf("Mass  M    %f\n",bigM());
  printf("Mass mA    %f\n",_mA);
  printf("Mass mB    %f\n",_mB);
  printf("Mass mC    %f\n",_mC);
  // limits
  printf("Limits qAB %f : %f\n",qAbsMin(AB),qAbsMax(AB));
  printf("Limits qBC %f : %f\n",qAbsMin(BC),qAbsMax(BC));
  printf("Limits qCA %f : %f\n",qAbsMin(CA),qAbsMax(CA));
  printf("Sum q       %f\n",sum());
  printf("Limit qsum  %f : %f\n",qSumMin(),qSumMax());
}
Commit	Line	Data
da0e9ce3	1	//-----------------------------------------------------------------------
da0e9ce3	2	// File and Version Information:
0ca57c2f	3	// $Id: EvtDalitzPlot.cpp,v 1.3 2009-03-16 15:53:27 robbep Exp $
da0e9ce3	4	//
	5	// Environment:
	6	// This software is part of the EvtGen package developed jointly
	7	// for the BaBar and CLEO collaborations. If you use all or part
	8	// of it, please give an appropriate acknowledgement.
	9	//
	10	// Copyright Information:
	11	// Copyright (C) 1998 Caltech, UCSB
	12	//
	13	// Module creator:
	14	// Alexei Dvoretskii, Caltech, 2001-2002.
	15	//-----------------------------------------------------------------------
	16	#include "EvtGenBase/EvtPatches.hh"
	17
	18	// Global 3-body Dalitz decay kinematics as defined by the mass
	19	// of the mother and the daughters. Spins are not considered.
	20
	21	#include <math.h>
	22	#include <assert.h>
	23	#include <stdio.h>
	24	#include "EvtGenBase/EvtPatches.hh"
	25	#include "EvtGenBase/EvtDalitzPlot.hh"
	26	#include "EvtGenBase/EvtTwoBodyVertex.hh"
	27	#include "EvtGenBase/EvtPDL.hh"
	28	#include "EvtGenBase/EvtDecayMode.hh"
	29
	30	using namespace EvtCyclic3;
	31
	32
	33	EvtDalitzPlot::EvtDalitzPlot()
	34	: _mA(0.), _mB(0.), _mC(0.), _bigM(0.),
	35	_ldel(0.), _rdel(0.)
	36	{}
	37
	38
	39	EvtDalitzPlot::EvtDalitzPlot(double mA, double mB, double mC, double bigM,
	40	double ldel, double rdel)
	41	: _mA(mA), _mB(mB), _mC(mC), _bigM(bigM),
	42	_ldel(ldel), _rdel(rdel)
	43	{
	44	sanityCheck();
	45	}
	46
	47	EvtDalitzPlot::EvtDalitzPlot(const EvtDecayMode& mode, double ldel, double rdel )
	48	{
	49	_mA = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(A)));
	50	_mB = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(B)));
	51	_mC = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(C)));
	52	_bigM = EvtPDL::getMeanMass(EvtPDL::getId(mode.mother()));
	53
	54	_ldel = ldel;
	55	_rdel = rdel;
	56
	57	sanityCheck();
	58	}
	59
	60
	61	EvtDalitzPlot::EvtDalitzPlot(const EvtDalitzPlot& other)
	62	: _mA(other._mA), _mB(other._mB), _mC(other._mC), _bigM(other._bigM),
	63	_ldel(other._ldel), _rdel(other._rdel)
	64	{}
	65
	66
	67	EvtDalitzPlot::~EvtDalitzPlot()
68	{}
69
70
71	bool EvtDalitzPlot::operator==(const EvtDalitzPlot& other) const
72	{
73	bool ret = false;
74	if(_mA == other._mA &&
75	_mB == other._mB &&
76	_mC == other._mC &&
77	_bigM == other._bigM) ret = true;
78
79	return ret;
80	}
81
82
83	const EvtDalitzPlot* EvtDalitzPlot::clone() const
84	{
85	return new EvtDalitzPlot(*this);
86	}
87
88
89	void EvtDalitzPlot::sanityCheck() const
90	{
91	if(_mA < 0 \|\| _mB < 0 \|\| _mC < 0 \|\| _bigM <= 0 \|\| _bigM - _mA - _mB - _mC < 0.) {
92
93	printf("Invalid Dalitz plot %f %f %f %f\n",_mA,_mB,_mC,_bigM);
94	assert(0);
95	}
96	assert(_ldel <= 0.);
97	assert(_rdel >= 0.);
98	}
99
100
101	double EvtDalitzPlot::m(Index i) const {
102
103	double m = _mA;
104	if(i == B) m = _mB;
105	else
106	if(i == C) m = _mC;
107
108	return m;
109	}
110
111
112	double EvtDalitzPlot::sum() const
113	{
114	return _mA_mA + _mB_mB + _mC_mC + _bigM_bigM;
115	}
116
117
118	double EvtDalitzPlot::qAbsMin(Pair i) const
119	{
120	Index j = first(i);
121	Index k = second(i);
122
123	return (m(j) + m(k))*(m(j) + m(k));
124	}
125
126
127	double EvtDalitzPlot::qAbsMax(Pair i) const
128	{
129	Index j = other(i);
130	return (_bigM-m(j))*(_bigM-m(j));
131	}
132
133
134	double EvtDalitzPlot::qResAbsMin(EvtCyclic3::Pair i) const
135	{
136	return qAbsMin(i) - sum()/3.;
137	}
138
139	double EvtDalitzPlot::qResAbsMax(EvtCyclic3::Pair i) const
140	{
141	return qAbsMax(i) - sum()/3.;
142	}
143
144	double EvtDalitzPlot::qHelAbsMin(EvtCyclic3::Pair i) const
145	{
146	Pair j = next(i);
147	Pair k = prev(i);
148	return (qAbsMin(j) - qAbsMax(k))/2.;
149	}
150
151	double EvtDalitzPlot::qHelAbsMax(EvtCyclic3::Pair i) const
152	{
153	Pair j = next(i);
154	Pair k = prev(i);
155	return (qAbsMax(j) - qAbsMin(k))/2.;
156	}
157
158
159	double EvtDalitzPlot::mAbsMin(Pair i) const
160	{
161	return sqrt(qAbsMin(i));
162	}
163
164
165	double EvtDalitzPlot::mAbsMax(Pair i) const
166	{
167	return sqrt(qAbsMax(i));
168	}
169
170
171	// parallel
172
173	double EvtDalitzPlot::qMin(Pair i, Pair j, double q) const
174	{
175	if(i == j) return q;
176
177	else {
178
179	// Particle pair j defines the rest-frame
180	// 0 - particle common to r.f. and angle calculations
181	// 1 - particle belonging to r.f. but not angle
182	// 2 - particle not belonging to r.f.
183
184	Index k0 = common(i,j);
185	Index k2 = other(j);
186	Index k1 = other(k0,k2);
187
188	// Energy, momentum of particle common to rest-frame and angle
189	EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
190	double pk = jpair.p();
191	double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
192
193
194	// Energy and momentum of the other particle
195	EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
196	double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
197	double pj = mother.p(EvtTwoBodyKine::A);
198
199
200	// See PDG 34.4.3.1
201	return (ek+ej)(ek+ej) - (pk+pj)(pk+pj);
202	}
203	}
204
205
206	// antiparallel
207
208	double EvtDalitzPlot::qMax(Pair i, Pair j, double q) const
209	{
210
211	if(i == j) return q;
212	else {
213
214	// Particle pair j defines the rest-frame
215	// 0 - particle common to r.f. and angle calculations
216	// 1 - particle belonging to r.f. but not angle
217	// 2 - particle not belonging to r.f.
218
219	Index k0 = common(i,j);
220	Index k2 = other(j);
221	Index k1 = other(k0,k2);
222
223	// Energy, momentum of particle common to rest-frame and angle
224	EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
225	double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
226	double pk = jpair.p();
227
228	// Energy and momentum of the other particle
229	EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
230	double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
231	double pj = mother.p(EvtTwoBodyKine::A);
232
233
234	// See PDG 34.4.3.1
235	return (ek+ej)(ek+ej) - (pk-pj)(pk-pj);
236	}
237	}
238
239
240	double EvtDalitzPlot::getArea(int N, Pair i, Pair j) const
241	{
242	// Trapezoidal integral over qi. qj can be calculated.
243	// The first and the last point are zero, so they are not counted
244
245	double dh = (qAbsMax(i) - qAbsMin(i))/((double) N);
246	double sum = 0;
247
248	int ii;
249	for(ii=1;ii<N;ii++) {
250
251	double x = qAbsMin(i) + ii*dh;
252	double dy = qMax(j,i,x) - qMin(j,i,x);
253	sum += dy;
254	}
255
256	return sum * dh;
257	}
258
259
260	double EvtDalitzPlot::cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
261	{
262	if(i1 == i2) return 1.;
263
264	double qmax = qMax(i1,i2,q2);
265	double qmin = qMin(i1,i2,q2);
266
267	double cos = (qmax + qmin - 2*q1)/(qmax - qmin);
268
269	return cos;
270	}
271
272
273	double EvtDalitzPlot::e(Index i, Pair j, double q) const
274	{
275	if(i == other(j)) {
276
277	// i does not belong to pair j
278
279	return (bigM()bigM()-q-m(i)m(i))/2/sqrt(q);
280	}
281	else {
282
283	// i and k make pair j
284
285	Index k;
286	if(first(j) == i) k = second(j);
287	else k = first(j);
288
289	double e = (q + m(i)m(i) - m(k)m(k))/2/sqrt(q);
290	return e;
291	}
292	}
293
294
295	double EvtDalitzPlot::p(Index i, Pair j, double q) const
296	{
297	double en = e(i,j,q);
298	double p2 = enen - m(i)m(i);
299
300	if(p2 < 0) {
301	printf("Bad value of p2 %f %d %d %f %f\n",p2,i,j,en,m(i));
302	assert(0);
303	}
304
305	return sqrt(p2);
306	}
307
308
309	double EvtDalitzPlot::q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const
310	{
311	if(i1 == i2) return q2;
312
313	EvtCyclic3::Index f = first(i1);
314	EvtCyclic3::Index s = second(i1);
315	return m(f)m(f) + m(s)m(s) + 2e(f,i2,q2)e(s,i2,q2) - 2p(f,i2,q2)p(s,i2,q2)*cosTh;
316	}
317
318
319	double EvtDalitzPlot::jacobian(EvtCyclic3::Pair i, double q) const
320	{
321	return 2p(first(i),i,q)p(other(i),i,q); // J(BC) = 2pApB = 2pApC
322	}
323
324
325	EvtTwoBodyVertex EvtDalitzPlot::vD(Pair iRes, double m0, int L) const
326	{
327	return EvtTwoBodyVertex(m(first(iRes)),
328	m(second(iRes)),m0,L);
329	}
330
331
332	EvtTwoBodyVertex EvtDalitzPlot::vB(Pair iRes, double m0, int L) const
333	{
334	return EvtTwoBodyVertex(m0,m(other(iRes)),bigM(),L);
335	}
336
337
338	void EvtDalitzPlot::print() const
339	{
340	// masses
341	printf("Mass M %f\n",bigM());
342	printf("Mass mA %f\n",_mA);
343	printf("Mass mB %f\n",_mB);
344	printf("Mass mC %f\n",_mC);
345	// limits
346	printf("Limits qAB %f : %f\n",qAbsMin(AB),qAbsMax(AB));
347	printf("Limits qBC %f : %f\n",qAbsMin(BC),qAbsMax(BC));
348	printf("Limits qCA %f : %f\n",qAbsMin(CA),qAbsMax(CA));
349	printf("Sum q %f\n",sum());
350	printf("Limit qsum %f : %f\n",qSumMin(),qSumMax());
351	}
352
353