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

#include "EvtGenBase/EvtPatches.hh"
/*******************************************************************************
 * Project: BaBar detector at the SLAC PEP-II B-factory
 * Package: EvtGenBase
 *    File: $Id: EvtDalitzResPdf.cpp,v 1.3 2009-03-16 15:54:07 robbep Exp $
 *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
 *
 * Copyright (C) 2002 Caltech
 *******************************************************************************/

#include <stdio.h>
#include <math.h>
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtDalitzResPdf.hh"
#include "EvtGenBase/EvtDalitzCoord.hh"
#include "EvtGenBase/EvtRandom.hh"
using namespace EvtCyclic3;

EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzPlot& dp, 
				 double _m0, double _g0, EvtCyclic3::Pair pair)
  : EvtPdf<EvtDalitzPoint>(), 
  _dp(dp), _m0(_m0), _g0(_g0), _pair(pair)
{}


EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzResPdf& other)
  : EvtPdf<EvtDalitzPoint>(other), 
  _dp(other._dp),_m0(other._m0), _g0(other._g0), _pair(other._pair)
{}


EvtDalitzResPdf::~EvtDalitzResPdf()
{}

EvtValError EvtDalitzResPdf::compute_integral(int N) const
{
  assert(N != 0);
  
  EvtCyclic3::Pair i = _pair;
  EvtCyclic3::Pair j = EvtCyclic3::next(i);

  // Trapezoidal integral

  double dh = (_dp.qAbsMax(j) - _dp.qAbsMin(j))/((double) N);
  double sum = 0;
  
  int ii;
  for(ii=1;ii<N;ii++) {
    
    double x = _dp.qAbsMin(j) + ii*dh;
    double min = (_dp.qMin(i,j,x) - _m0*_m0)/_m0/_g0;
    double max = (_dp.qMax(i,j,x) - _m0*_m0)/_m0/_g0;
    double itg = 1/EvtConst::pi*(atan(max) - atan(min));
    sum += itg;
  }
  EvtValError ret(sum*dh,0.); 
  
  return ret;
}


EvtDalitzPoint EvtDalitzResPdf::randomPoint()
{
  // Random point generation must be done in a box encompassing the 
  // Dalitz plot


  EvtCyclic3::Pair i = _pair;
  EvtCyclic3::Pair j = EvtCyclic3::next(i);  
  double min = 1/EvtConst::pi*atan((_dp.qAbsMin(i) - _m0*_m0)/_m0/_g0);
  double max = 1/EvtConst::pi*atan((_dp.qAbsMax(i) - _m0*_m0)/_m0/_g0);

  int n = 0;
  while(n++ < 1000) {

    double qj = EvtRandom::Flat(_dp.qAbsMin(j),_dp.qAbsMax(j));
    double r = EvtRandom::Flat(min,max);
    double qi = tan(EvtConst::pi*r)*_g0*_m0 + _m0*_m0;
    EvtDalitzCoord x(i,qi,j,qj);
    EvtDalitzPoint ret(_dp,x);
    if(ret.isValid()) return ret;
  }
  
  // All generated points turned out to be outside of the Dalitz plot
  // (in the outer box)
  
  printf("No point generated for dalitz plot after 1000 tries\n");
  return EvtDalitzPoint(0.,0.,0.,0.,0.,0.);
}


double EvtDalitzResPdf::pdf(const EvtDalitzPoint& x) const
{
  EvtCyclic3::Pair i = _pair;
  double dq = x.q(i) - _m0*_m0;
  return 1/EvtConst::pi*_g0*_m0/(dq*dq + _g0*_g0*_m0*_m0);
}


double EvtDalitzResPdf::pdfMaxValue() const
{
  return 1/(EvtConst::pi*_g0*_m0);
}
Commit	Line	Data
da0e9ce3	1	#include "EvtGenBase/EvtPatches.hh"
	2	/*******************************************************************************
	3	* Project: BaBar detector at the SLAC PEP-II B-factory
	4	* Package: EvtGenBase
0ca57c2f	5	* File: $Id: EvtDalitzResPdf.cpp,v 1.3 2009-03-16 15:54:07 robbep Exp $
da0e9ce3	6	* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
	7	*
	8	* Copyright (C) 2002 Caltech
	9	*******************************************************************************/
	10
	11	#include <stdio.h>
	12	#include <math.h>
	13	#include "EvtGenBase/EvtPatches.hh"
	14	#include "EvtGenBase/EvtConst.hh"
	15	#include "EvtGenBase/EvtDalitzResPdf.hh"
	16	#include "EvtGenBase/EvtDalitzCoord.hh"
	17	#include "EvtGenBase/EvtRandom.hh"
	18	using namespace EvtCyclic3;
	19
	20	EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzPlot& dp,
	21	double _m0, double _g0, EvtCyclic3::Pair pair)
	22	: EvtPdf<EvtDalitzPoint>(),
	23	_dp(dp), _m0(_m0), _g0(_g0), _pair(pair)
	24	{}
	25
	26
	27	EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzResPdf& other)
	28	: EvtPdf<EvtDalitzPoint>(other),
	29	_dp(other._dp),_m0(other._m0), _g0(other._g0), _pair(other._pair)
	30	{}
	31
	32
	33	EvtDalitzResPdf::~EvtDalitzResPdf()
	34	{}
	35
	36	EvtValError EvtDalitzResPdf::compute_integral(int N) const
	37	{
	38	assert(N != 0);
	39
	40	EvtCyclic3::Pair i = _pair;
	41	EvtCyclic3::Pair j = EvtCyclic3::next(i);
	42
	43	// Trapezoidal integral
	44
	45	double dh = (_dp.qAbsMax(j) - _dp.qAbsMin(j))/((double) N);
	46	double sum = 0;
	47
	48	int ii;
	49	for(ii=1;ii<N;ii++) {
	50
	51	double x = _dp.qAbsMin(j) + ii*dh;
	52	double min = (_dp.qMin(i,j,x) - _m0*_m0)/_m0/_g0;
	53	double max = (_dp.qMax(i,j,x) - _m0*_m0)/_m0/_g0;
	54	double itg = 1/EvtConst::pi*(atan(max) - atan(min));
	55	sum += itg;
	56	}
	57	EvtValError ret(sum*dh,0.);
	58
	59	return ret;
	60	}
	61
	62
	63	EvtDalitzPoint EvtDalitzResPdf::randomPoint()
	64	{
	65	// Random point generation must be done in a box encompassing the
	66	// Dalitz plot
	67
	68
	69	EvtCyclic3::Pair i = _pair;
70	EvtCyclic3::Pair j = EvtCyclic3::next(i);
71	double min = 1/EvtConst::piatan((_dp.qAbsMin(i) - _m0_m0)/_m0/_g0);
72	double max = 1/EvtConst::piatan((_dp.qAbsMax(i) - _m0_m0)/_m0/_g0);
73
74	int n = 0;
75	while(n++ < 1000) {
76
77	double qj = EvtRandom::Flat(_dp.qAbsMin(j),_dp.qAbsMax(j));
78	double r = EvtRandom::Flat(min,max);
79	double qi = tan(EvtConst::pir)_g0_m0 + _m0_m0;
80	EvtDalitzCoord x(i,qi,j,qj);
81	EvtDalitzPoint ret(_dp,x);
82	if(ret.isValid()) return ret;
83	}
84
85	// All generated points turned out to be outside of the Dalitz plot
86	// (in the outer box)
87
88	printf("No point generated for dalitz plot after 1000 tries\n");
0ca57c2f	89	return EvtDalitzPoint(0.,0.,0.,0.,0.,0.);
da0e9ce3	90	}
	91
	92
	93	double EvtDalitzResPdf::pdf(const EvtDalitzPoint& x) const
	94	{
	95	EvtCyclic3::Pair i = _pair;
	96	double dq = x.q(i) - _m0*_m0;
	97	return 1/EvtConst::pi_g0_m0/(dqdq + _g0_g0_m0_m0);
	98	}
	99
	100
	101	double EvtDalitzResPdf::pdfMaxValue() const
	102	{
	103	return 1/(EvtConst::pi_g0_m0);
	104	}
	105