Page MenuHomeHEPForge
Files F19251926

EvtDalitzResPdf.cpp
No OneTemporary
Actions

Size
3 KB
Referenced Files
None
Subscribers
None

EvtDalitzResPdf.cpp
View Options

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#include "EvtGenBase/EvtDalitzResPdf.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtDalitzCoord.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtRandom.hh"
#include <math.h>
#include <stdio.h>
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 )
{
}
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 );
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:12 AM (19 h, 7 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6552834
Default Alt Text
EvtDalitzResPdf.cpp (3 KB)

Event Timeline

HEPForge.org · .