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
3 changes: 2 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,9 +178,10 @@ add_subdirectory(examples/random)
add_subdirectory(examples/histograming)
add_subdirectory(examples/async)
add_subdirectory(examples/misc)
add_subdirectory(examples/phys)

if(Minuit2_FOUND)
add_subdirectory(examples/phys)
Copy link
Contributor

Choose a reason for hiding this comment

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

All examples in the directory '''phys''' depends on Minuit2.
You need to put it inside the guards ''' if(Minuit2_FOUND) ... endif(Minuit2_FOUND)'''


add_subdirectory(examples/fit)

endif(Minuit2_FOUND)
Expand Down
64 changes: 54 additions & 10 deletions examples/phys/FlatteLineShape.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -350,7 +350,7 @@ int main(int argv, char** argc)

double f0_MASS = 0.965;
double f0_MAG = 12.341;
double f0_Phase = -62.852*(M_PI / 180) ;
Copy link
Contributor

Choose a reason for hiding this comment

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

Document the numerical values, please.

double f0_Phase = -62.852*(M_PI / 180) +M_PI;
Copy link
Contributor

Choose a reason for hiding this comment

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

The primary documentation of code, is the code itself. So, please enforce double precision syntax where you are using it. For example '''180''' -> '''180.0'''

double f0_RC = f0_MAG * cos(f0_Phase);
double f0_IMC = f0_MAG * sin(f0_Phase);
double f0_rho1 = 0.165;
Expand Down Expand Up @@ -381,9 +381,8 @@ int main(int argv, char** argc)
//======================================================

//f0

auto coef_ref0 = hydra::Parameter::Create().Name("f0_RC").Value(f0_RC).Error(0.0001);
auto coef_imf0 = hydra::Parameter::Create().Name("f0_IM").Value(f0_IMC).Error(0.0001);
auto coef_ref0 = hydra::Parameter::Create().Name("f0_RC").Value(f0_RC).Error(0.0001).Limits(-100,100);
Copy link
Contributor

Choose a reason for hiding this comment

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

Define the limits in the parameters in terms of the nominal values in the parameters.

auto coef_imf0 = hydra::Parameter::Create().Name("f0_IM").Value(f0_IMC).Error(0.0001).Limits(-100,100);
auto f0Mass = hydra::Parameter::Create().Name("MASS_f0").Value(f0_MASS).Fixed();
auto f0g1 = hydra::Parameter::Create().Name("f0_g1").Value(f0_rho1).Fixed();
auto rg1og2 = hydra::Parameter::Create().Name("f0_g1xg2").Value(f0_rho2).Fixed();
Expand All @@ -403,7 +402,8 @@ int main(int argv, char** argc)

hydra::complex<double> r(0,0);

for(unsigned int i=0; i< n;i++) r += x[i];
for(unsigned int i=0; i<n; i++)
r += x[i];

return hydra::norm(r);
});
Expand Down Expand Up @@ -464,7 +464,8 @@ int main(int argv, char** argc)
f0_12_HIST,f0_13_HIST ;

double Phi_12_FF, Phi_13_FF,
f0_12_FF,f0_13_FF
f0_12_FF,f0_13_FF,
Phi_all_FF, f0_all_FF
;
#endif

Expand Down Expand Up @@ -499,6 +500,35 @@ int main(int argv, char** argc)

std::cout << std::endl <<"Toy Dataset size: "<< toy_data.size() << std::endl;

std::ofstream writer;
writer.open("toyData.txt");

if(!writer.is_open()){
std::cout << "file not open" << std::endl;
}else {

for (auto event : toy_data) {

double weight = hydra::get<0>(event);
hydra::Vector4R Kminus = hydra::get<1>(event);
hydra::Vector4R Kplus1 = hydra::get<2>(event);
hydra::Vector4R Kplus2 = hydra::get<3>(event);

double MKminusKplus1 = (Kminus + Kplus1).mass2();
double MKminusKplus2 = (Kminus + Kplus2).mass2();

writer << MKminusKplus1 << '\t' << MKminusKplus2 << std::endl;
}

std::cout << "toyfile.txt generated" <<std::endl;
}



writer.close();



}//toy data production on device

//plot toy-data
Expand Down Expand Up @@ -712,8 +742,11 @@ int main(int argv, char** argc)

auto Phi_12 = fcn.GetPDF().GetFunctor().GetFunctor(_1).GetFunctor(_0);
auto Phi_13 = fcn.GetPDF().GetFunctor().GetFunctor(_1).GetFunctor(_1);
auto Phi_all = fcn.GetPDF().GetFunctor().GetFunctor(_1);
auto f0_12 = fcn.GetPDF().GetFunctor().GetFunctor(_2).GetFunctor(_0);
auto f0_13 = fcn.GetPDF().GetFunctor().GetFunctor(_2).GetFunctor(_1);
auto f0_all = fcn.GetPDF().GetFunctor().GetFunctor(_2);


//==================================
// Draw components
Expand All @@ -729,16 +762,27 @@ int main(int argv, char** argc)
//==================================
Phi_12_FF = fit_fraction(Phi_12 , Opt_Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries);
Phi_13_FF = fit_fraction(Phi_13 , Opt_Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries);
Phi_all_FF = fit_fraction(Phi_all , Opt_Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries);
f0_12_FF = fit_fraction(f0_12 , Opt_Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries);
f0_13_FF = fit_fraction(f0_13 , Opt_Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries);
f0_all_FF = fit_fraction(f0_all , Opt_Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries);

std::cout << "Phi_12_FF :" << Phi_12_FF << std::endl;
std::cout << "Phi_13_FF :" << Phi_13_FF << std::endl;
std::cout << "Phi_all_FF :" << Phi_all_FF << std::endl;
std::cout << "f0_12_FF :" << f0_12_FF << std::endl;
std::cout << "f0_13_FF :" << f0_13_FF << std::endl;
std::cout << "f0_all_FF :" << f0_all_FF << std::endl;
std::cout << "Sum :"
<< Phi_12_FF + Phi_13_FF + f0_12_FF + f0_13_FF << std::endl;

std::cout << "Phi_12_FF :" << fit_fraction(Phi_Resonance_12, Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries) << std::endl;
std::cout << "Phi_13_FF :" << fit_fraction(Phi_Resonance_13, Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries) << std::endl;
std::cout << "Phi_FF :" << fit_fraction(Phi_Resonance, Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries) << std::endl;
std::cout << "f0_12_FF: " << fit_fraction(f0_Resonance_12, Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries) << std::endl;
std::cout << "f0_13_FF: " << fit_fraction(f0_Resonance_13, Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries) << std::endl;
std::cout << "f0_FF: " << fit_fraction(f0_Resonance, Model, {D_MASS, Kminus_MASS, Kplus_MASS}, nentries) << std::endl;

#ifdef _ROOT_AVAILABLE_

{
Expand Down Expand Up @@ -815,12 +859,12 @@ int main(int argv, char** argc)

Dalitz_Fit.Scale(Dalitz_Resonances.Integral()/Dalitz_Fit.Integral() );

Phi_12_HIST.Scale( Phi_12_FF*Dalitz_Fit.Integral()/Phi_12_HIST.Integral() );
Phi_13_HIST.Scale( Phi_13_FF*Dalitz_Fit.Integral()/Phi_13_HIST.Integral() );
Phi_12_HIST.Scale(Phi_12_FF*Dalitz_Fit.Integral()/Phi_12_HIST.Integral() );
Phi_13_HIST.Scale(Phi_13_FF*Dalitz_Fit.Integral()/Phi_13_HIST.Integral() );
f0_12_HIST.Scale( f0_12_FF*Dalitz_Fit.Integral()/f0_12_HIST.Integral() );
f0_13_HIST.Scale( f0_13_FF*Dalitz_Fit.Integral()/f0_13_HIST.Integral() );

//=============================================================
//=============================================================
//projections
TH1* hist2D=0;

Expand Down Expand Up @@ -1133,6 +1177,7 @@ double fit_fraction( Amplitude const& amp, Model const& model, std::array<double
// Create PhaseSpace object for D-> K K K
hydra::PhaseSpace<3> phsp{Kminus_MASS, Kplus_MASS, Kplus_MASS};


//norm lambda
auto Norm = hydra::wrap_lambda( []__host__ __device__ (unsigned int n, hydra::complex<double>* x){

Expand All @@ -1146,7 +1191,6 @@ double fit_fraction( Amplitude const& amp, Model const& model, std::array<double
auto amp_int = phsp.AverageOn(hydra::device::sys, D, functor, nentries);
auto model_int = phsp.AverageOn(hydra::device::sys, D, model, nentries);


return amp_int.first/model_int.first;

}
Expand Down
4 changes: 2 additions & 2 deletions hydra/functions/FlatteLineShape.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -218,9 +218,9 @@ namespace hydra {


if(s >= twopimasssq)
rhopipi_real = (2. / 3) * TMath::Sqrt(1 - twopimasssq / s ); // Above pi+pi- threshold
rhopipi_real = (2. / 3) * ::sqrt(1 - twopimasssq / s ); // Above pi+pi- threshold
Copy link
Contributor

@AAAlvesJr AAAlvesJr Mar 19, 2018 •
edited
Loading

Choose a reason for hiding this comment

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

This whole piece of code could generate branching and does unnecessary variable redefinition, which makes execution slow down on gpus. I consider a defect of C++ that variables are not default const, like in Rust, but this is another discussion. The code is also confusing. So deploy Boolean arithmetic here. Example:

instead to write

 if(condition) x = a+b;
else x=a-b;

write (1 line !)

x = a + (condition==true) *(b) - (condition!=true) *(b)

and so on...
Please, reflect about this.

else
rhopipi_imag = (2. / 3) * TMath::Sqrt(-1 + twopimasssq / s);
rhopipi_imag = (2. / 3) * ::sqrt(-1 + twopimasssq / s);
Copy link
Contributor

Choose a reason for hiding this comment

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

idem...

Copy link
Contributor

Choose a reason for hiding this comment

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

Please, try to eliminate the if/else instructions using boolean arithmetic. Avoid value., enforce value.0.


if(s >= twopi0masssq)
rhopipi_real = (1. / 3) * TMath::Sqrt(1 - twopi0masssq / s ); // Above pi0pi0 threshold
Expand Down