Skip to content

Develop #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 16 commits into
base: develop
Choose a base branch
from
Open

Develop #25

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
2bcf807
first commit - Flatte Lineshape
JuanBSLeite Mar 7, 2018
a13c230
testing flatté
JuanBSLeite Mar 7, 2018
7eaaf92
flatte lineshape
JuanBSLeite Mar 9, 2018
301b2c2
clear unused code
JuanBSLeite Mar 9, 2018
1cffbcd
cleaning some code
JuanBSLeite Mar 10, 2018
6d59a21
flatte
JuanBSLeite Mar 17, 2018
6413653
changes for PR
JuanBSLeite Mar 19, 2018
311ffa2
changes for PR
JuanBSLeite Mar 19, 2018
771d054
last change to PR
JuanBSLeite Mar 19, 2018
1299a65
if/else to boolean operations
JuanBSLeite Mar 22, 2018
d7550a5
How to use the Bolean algebra...
AAAlvesJr Mar 22, 2018
d5de8fa
boolean operation
JuanBSLeite Mar 23, 2018
123edd0
negative sqrt correction/.//
AAAlvesJr Mar 23, 2018
3a2fe66
FlatteLineshape based in EvtGen
JuanBSLeite Apr 4, 2018
49d88c9
Flatte running in cuda backend
JuanBSLeite Apr 12, 2018
8bba5bd
flatte unused paramaters removed
JuanBSLeite Apr 12, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 16 additions & 15 deletions examples/phys/FlatteLineShape.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,27 +238,28 @@ public:

};

template<unsigned int CHANNEL,hydra::Wave L> class Flatte: public hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 5>
template<unsigned int CHANNEL,hydra::Wave L> class Flatte: public hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 3>
{

constexpr static unsigned int _I1 = CHANNEL-1;
constexpr static unsigned int _I2 = (CHANNEL!=3)*CHANNEL;
constexpr static unsigned int _I3 = 3-( (CHANNEL-1) + (CHANNEL!=3)*CHANNEL );

using hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 5>::_par;
using hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 3>::_par;

public:

Flatte() = delete;

Flatte(hydra::Parameter const& c_re, hydra::Parameter const& c_im, hydra::Parameter const& mean, hydra::Parameter const& rho1 , hydra::Parameter const& rho2,
double mother_mass, double daugther1_mass, double daugther2_mass, double daugther3_mass, double radi):
hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 5>{c_re, c_im, mean, rho1, rho2},
fLineShape(mean,rho1,rho2,mother_mass,daugther1_mass,daugther2_mass,daugther3_mass,radi)
Flatte(hydra::Parameter const& c_re, hydra::Parameter const& c_im, hydra::Parameter const& mean,
std::vector<std::vector<double>> const& params,double mother_mass, double daugther1_mass, double daugther2_mass,
double daugther3_mass, double radi):
hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 3>{c_re, c_im, mean},
fLineShape(mean,params,mother_mass,daugther1_mass,daugther2_mass,daugther3_mass,radi)
{}

__hydra_dual__ Flatte( Flatte<CHANNEL,L> const& other):
hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 5>(other),
hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 3>(other),
fLineShape(other.GetLineShape())
{}

Expand All @@ -268,7 +269,7 @@ public:
{
if(this==&other) return *this;

hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 5>::operator=(other);
hydra::BaseFunctor<Flatte<CHANNEL,L>, hydra::complex<double>, 3>::operator=(other);
fLineShape=other.GetLineShape();

return *this;
Expand All @@ -285,15 +286,12 @@ public:
hydra::Vector4R p3 = p[_I3];

fLineShape.SetParameter(0,_par[2]);
fLineShape.SetParameter(1,_par[3]);
fLineShape.SetParameter(2,_par[4]);

double theta = fCosDecayAngle( (p1+p2+p3), (p1+p2), p1 );
double angular = fAngularDist(theta);

auto r = hydra::complex<double>(_par[0],_par[1])*fLineShape((p1+p2).mass())*angular;

return r;
return r;
}

private:
Expand Down Expand Up @@ -360,8 +358,11 @@ int main(int argv, char** argc)
double Kplus_MASS = 0.493677; // K+ mass
double Kminus_MASS = Kplus_MASS;

//Flatté
double pi_MASS = 0.13957018;
std::vector<std::vector<double>> params = {{pi_MASS,pi_MASS,f0_rho1},{Kplus_MASS,Kplus_MASS,f0_rho2}};

//======================================================
//======================================================
//Phi
auto mass = hydra::Parameter::Create().Name("MASS_Phi" ).Value(Phi_MASS ).Error(0.01).Fixed();
auto width = hydra::Parameter::Create().Name("WIDTH_Phi").Value(Phi_Width).Error(0.001).Fixed();
Expand All @@ -388,8 +389,8 @@ int main(int argv, char** argc)
auto rg1og2 = hydra::Parameter::Create().Name("f0_g1xg2").Value(f0_rho2).Fixed();


Flatte<1,hydra::SWave> f0_Resonance_12(coef_ref0,coef_imf0,f0Mass,f0g1,rg1og2,D_MASS,Kminus_MASS,Kplus_MASS,Kplus_MASS,1.5);
Flatte<3,hydra::SWave> f0_Resonance_13(coef_ref0,coef_imf0,f0Mass,f0g1,rg1og2,D_MASS,Kminus_MASS,Kplus_MASS,Kplus_MASS,1.5);
Flatte<1,hydra::SWave> f0_Resonance_12(coef_ref0,coef_imf0,f0Mass,params,D_MASS,Kminus_MASS,Kplus_MASS,Kplus_MASS,1.5);
Flatte<3,hydra::SWave> f0_Resonance_13(coef_ref0,coef_imf0,f0Mass,params,D_MASS,Kminus_MASS,Kplus_MASS,Kplus_MASS,1.5);

auto f0_Resonance = (f0_Resonance_12 + f0_Resonance_13);

Expand Down
105 changes: 44 additions & 61 deletions hydra/functions/FlatteLineShape.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
#include <hydra/Complex.h>
#include <hydra/functions/Utils.h>
#include <hydra/functions/BlattWeisskopfFunctions.h>
#include <hydra/Vector4R.h>
#include <tuple>
#include <limits>
#include <stdexcept>
Expand All @@ -61,8 +62,8 @@ namespace hydra {


template<unsigned int CHANNEL, unsigned int ArgIndex = 0>
class FlatteLineShape : public BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 3> {
using BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 3>::_par;
class FlatteLineShape : public BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 1> {
using BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 1>::_par;

public:

Expand All @@ -78,11 +79,13 @@ namespace hydra {
* @param daugther3_mass daughter particle 2 mass
* @param radi decay vertex radio.
*/
FlatteLineShape(hydra::Parameter const& mean, hydra::Parameter const& rho1 , hydra::Parameter const& rho2,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the length of the vector std::vector<std::vector<double>> params and what does it holds?
I think can be a good idea to document the meaning of this parameter.
Also, are these parameters supposed to be optimized in a fit? Case yes, consider to use the correct interface to store and make them "trackable" by Minuit2 and copy constructors on device and host sides.
If the number of such parameters is not fixed, but still constant at compile time, please, add a template parameter to the class to implement then, and set the default value to the most common case (my guess is most of the time it is 2...) . After set them using this->SetParameter(param_index, param_value) in the body of the constructor.

FlatteLineShape(hydra::Parameter const& mean,
std::vector<std::vector<double>> params,
double mother_mass,
double daugther1_mass, double daugther2_mass, double daugther3_mass,
double radi) :
BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 3>{mean, rho1, rho2},
BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>,1>{mean},
fParams(params),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unfortunately, this implementation does not fit the Hydra's requirements, since this class is not copyable on device side for CUDA backends. This happens because you are using a std::vector as a member. This will works on host and device CPU-based backends, but not in CUDA, because std::vector are not implemented there...

fDaughter1Mass(daugther1_mass),
fDaughter2Mass(daugther2_mass),
fBachelorMass(daugther3_mass),
Expand All @@ -91,7 +94,8 @@ namespace hydra {

__hydra_host__ __hydra_device__
FlatteLineShape(FlatteLineShape<CHANNEL, ArgIndex> const &other) :
BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 3>(other),
BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex<double>, 1>(other),
fParams(other.GetParams()),
fDaughter1Mass(other.GetDaughter1Mass()),
fDaughter2Mass(other.GetDaughter2Mass()),
fBachelorMass(other.GetDaughter3Mass()),
Expand All @@ -104,8 +108,9 @@ namespace hydra {
if (this == &other) return *this;

BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>,
hydra::complex<double>, 3>::operator=(other);
hydra::complex<double>, 7>::operator=(other);

fParams = other.GetParams();
fDaughter1Mass = other.GetDaughter1Mass();
fDaughter2Mass = other.GetDaughter2Mass();
fBachelorMass = other.GetDaughter3Mass();
Expand All @@ -115,6 +120,18 @@ namespace hydra {
return *this;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

L.110 : BaseFunctor<FlatteLineShape<CHANNEL, ArgIndex>, hydra::complex, 7>::operator=(other);**
This should be 1 -------------------------------------------------------------------------------------------^^^

This assignment operator can't compile either.


__hydra_host__ __hydra_device__ inline
std::vector<std::vector<double>>GetParams() const {
return fParams;
}

__hydra_host__ __hydra_device__ inline
void SetDaughter1Mass(std::vector<std::vector<double>> _params) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it seems to me this member is setting the wrong class member... Should be fDaughterMass1 instead of fParams, right?

fParams = _params;
}



__hydra_host__ __hydra_device__ inline
double GetDaughter1Mass() const {
return fDaughter1Mass;
Expand Down Expand Up @@ -170,84 +187,49 @@ namespace hydra {
hydra::complex<double> Evaluate(unsigned int, T *x) const {

const double m = x[ArgIndex];

const double resonance_mass = _par[0];
const double rho1 = _par[1];
const double rho2 = _par[2];

return m > (fDaughter1Mass + fDaughter2Mass) && m < (fMotherMass - fBachelorMass) ?
LineShape(m, resonance_mass, rho1, rho2) : hydra::complex<double>(0.0, 0.0);
return LineShape( m , resonance_mass) ;

}

template<typename T>
__hydra_host__ __hydra_device__ inline
hydra::complex<double> Evaluate(T x) const {

double m = get<ArgIndex>(x);

const double m = get<ArgIndex>(x);
const double resonance_mass = _par[0];
const double rho1 = _par[1];
const double rho2 = _par[2];

return m > (fDaughter1Mass + fDaughter2Mass) && m < (fMotherMass - fBachelorMass) ?
LineShape(m, resonance_mass, rho1, rho2) : hydra::complex<double>(0.0, 0.0);
return LineShape( m , resonance_mass) ;
}

private:
__hydra_host__ __hydra_device__ inline
hydra::complex<double>
sqrtCplx(const double in) const { return (in > 0) ? hydra::complex<double>(::sqrt(in), 0.0) : hydra::complex<double>(0.0, ::sqrt(-in)); }

__hydra_dual__ inline
hydra::complex<double> LineShape(const double s, const double resonance_mass, const double g1 , const double g2) const {

double pipmass = 0.13957018;
double pi0mass = 0.1349766;
double kpmass = 0.493677;
double k0mass = 0.497614;

double twopimasssq = 4 * pipmass * pipmass;
double twopi0masssq = 4 * pi0mass * pi0mass;
double twokmasssq = 4 * kpmass * kpmass;
double twok0masssq = 4 * k0mass * k0mass;

double rhokk_real = 0;
double rhokk_imag = 0;
double rhopipi_real = 0;
double rhopipi_imag = 0;


constexpr static double _1p3 = 0.333333333;
constexpr static double _2p3 = 0.666666667;
static const double two_pi0_factor = ::sqrt(1 - twopi0masssq / s );
static const double two_pi_factor = ::sqrt(1 - twopimasssq / s );
static const double two_pi0_factor_img = ::sqrt((-1 + twopi0masssq / s)*(twopi0masssq >= s) );
static const double two_pi_factor_img = ::sqrt((-1 + twopimasssq / s)*(twopimasssq >= s)) );

rhopipi_real = (s >= twopi0masssq)*_1p3*twopi0masssq;

rhopipi_imag = (s >= twopi0masssq && s < twokmasssq)* (_2p3 * two_pi_factor_img )+
(s < twopi0masssq)*_1p3*two_pi0_factor_img;

constexpr static double _1p2 = 0.5;
static const double two_k0_factor = ::sqrt(1 - twok0masssq / s );
static const double two_k_factor = ::sqrt(1 - twokmasssq / s );
static const double two_k0_factor_img = ::sqrt((-1 + twok0masssq / s)*(twok0masssq>=s) );
static const double two_k_factor_img = ::sqrt((-1 + twokmasssq / s)*(twokmasssq>=s) );
__hydra_host__ __hydra_device__ inline
hydra::complex<double>
LineShape(const double s, const double resonance_mass) const {

hydra::complex<double> w;

rhokk_real = (s >= twokmasssq && s < twok0masssq)*_1p2*two_k_factor +
(s >= twok0masssq)*_1p2*two_k0_factor;
for(size_t i = 0; i < fParams.size() ; i++) {

rhokk_imag = (s < twok0masssq)*_1p2*two_k0_factor_img;
double m1a = fParams[i][0];
double m1b = fParams[i][1];
double g = fParams[i][2];

w += g * g * sqrtCplx((1 - (m1a - m1b) * (m1a - m1b) / s*s ) *
(1 - (m1a + m1b) * (m1a + m1b) / s*s ));

double A = (resonance_mass*resonance_mass - s) + resonance_mass*(rhopipi_imag*g1 + rhokk_imag*g2);
double B = resonance_mass*(rhopipi_real*g1 + rhokk_real*g2);
double C = 1.0 / (A*A + B*B);
}

hydra::complex<double> retur(A*C, B*C);
hydra::complex<double> denom = resonance_mass - s*s - hydra::complex<double>(0.0,1.0)*w;

return retur;
hydra::complex<double> ampl = hydra::complex<double>(1.0,0.0)/denom;

return ampl;

}

Expand All @@ -256,6 +238,7 @@ namespace hydra {
double fBachelorMass;
double fMotherMass;
double fRadi;
std::vector<std::vector<double>> fParams;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this member is not copyable on CUDA...


};

Expand Down