Hydra: dalitz_plot.C (original) (raw)

#include

#include <assert.h>

#include <time.h>

#include

#ifndef HYDRA_HOST_SYSTEM

#define HYDRA_HOST_SYSTEM CPP

#endif

#ifndef HYDRA_DEVICE_SYSTEM

#define HYDRA_DEVICE_SYSTEM TBB

#endif

#include "Minuit2/FunctionMinimum.h"

#include "Minuit2/MnUserParameterState.h"

#include "Minuit2/MnPrint.h"

#include "Minuit2/MnMigrad.h"

#include "Minuit2/MnMinimize.h"

#include <TROOT.h>

#include <TH1D.h>

#include <TH2D.h>

#include <TH3D.h>

#include <TApplication.h>

#include <TCanvas.h>

#include <TLegend.h>

#include <TLegendEntry.h>

declarg(Dplus, hydra::Vector4R)

declarg(PionA, hydra::Vector4R)

declarg(PionB, hydra::Vector4R)

declarg(Kaon , hydra::Vector4R)

template<hydra::Wave L, bool Flag=(L%2)>

template<hydra::Wave L>

struct parity<L, false>: std::integral_constant<int,1>{};

template<hydra::Wave L>

struct parity<L, true>: std::integral_constant<int,-1>{};

template<hydra::Wave L >

{

using super_type::_par;

public:

double mother_mass, double daugther1_mass, double daugther2_mass, double daugther3_mass, double radi):

super_type({c_re, c_im, mass, width}),

fLineShape(mass, width, mother_mass, daugther1_mass, daugther2_mass, daugther3_mass, radi)

{}

super_type(other),

fLineShape(other.GetLineShape())

{}

{

if(this==&other) return *this;

super_type::operator=(other);

return *this;

}

GetLineShape() const { return fLineShape; }

Evaluate(Kaon kaon, PionA pion1, PionB pion2) const {

hydra::Vector4R mother = kaon + pion1 + pion2;

hydra::Vector4R Kpi1 = kaon + pion1;

hydra::Vector4R Kpi2 = kaon + pion2;

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

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

hydra::complex contrib_12 = fLineShape((Kpi1).mass())*fAngularDist(fCosDecayAngle(mother, Kpi1, kaon));

hydra::complex contrib_13 = fLineShape((Kpi2).mass())*fAngularDist(fCosDecayAngle(mother, Kpi2, pion2));

return r;

}

private:

};

{

using super_type::_par;

public:

super_type({c_re, c_im})

{}

super_type(other)

{}

{

if(this==&other) return *this;

super_type::operator=(other);

return *this;

}

};

template<typename ...T>

{

public:

super_type(other)

{}

{

if(this==&other) return *this;

super_type::operator=(other);

return *this;

}

inline double Evaluate( T... A ) const {

auto r=add({A...});

};

private:

template

inline C add( std::initializer_list list ) const

{

C r{};

for( auto x: list)

r+=x;

return r;

}

};

template<typename ...Amplitudes>

auto make_model( Amplitudes const& ... amplitudes)

{

}

template

template<typename Amplitude, typename Model>

template<typename Backend, typename Model, typename Container >

size_t generate_dataset(Backend const& system, Model const& model, std::array<double, 3> const& masses, Container& decays, size_t nevents, size_t bunch_size);

{

double NR_MAG = 7.4;

double NR_PHI = (-18.4+180.0)*0.01745329;

double NR_CRe = NR_MAG*cos(NR_PHI);

double NR_CIm = NR_MAG*sin(NR_PHI);

double K800_MASS = 0.809 ;

double K800_WIDTH = 0.470;

double K800_MAG = 5.01;

double K800_PHI = (-163.7+180.0)*0.01745329;

double K800_CRe = K800_MAG*cos(K800_PHI);

double K800_CIm = K800_MAG*sin(K800_PHI);

double KST_892_MASS = 0.89555;

double KST_892_WIDTH = 0.0473;

double KST_892_MAG = 1.0;

double KST_892_PHI = 0.0;

double KST_892_CRe = KST_892_MAG*cos(KST_892_PHI);

double KST_892_CIm = KST_892_MAG*sin(KST_892_PHI);

double KST0_1430_MASS = 1.425;

double KST0_1430_WIDTH = 0.270;

double KST0_1430_MAG = 3.5;

double KST0_1430_PHI = (49.7-180.0)*0.01745329;

double KST0_1430_CRe = KST0_1430_MAG*cos(KST0_1430_PHI);

double KST0_1430_CIm = KST0_1430_MAG*sin(KST0_1430_PHI);

double KST2_1430_MASS = 1.4324;

double KST2_1430_WIDTH = 0.109;

double KST2_1430_MAG = 0.962;

double KST2_1430_PHI = (-29.9+180.0)*0.01745329;

double KST2_1430_CRe = KST2_1430_MAG*cos(KST2_1430_PHI);

double KST2_1430_CIm = KST2_1430_MAG*sin(KST2_1430_PHI);

double KST_1680_MASS = 1.718;

double KST_1680_WIDTH = 0.322;

double KST_1680_MAG = 6.5;

double KST_1680_PHI = (29.0)*0.01745329;

double KST_1680_CRe = KST_1680_MAG*cos(KST_1680_PHI);

double KST_1680_CIm = KST_1680_MAG*sin(KST_1680_PHI);

double D_MASS = 1.86959;

double K_MASS = 0.493677;

double PI_MASS = 0.13957061;

auto Model = make_model( K800_Resonance, KST_892_Resonance,

KST0_1430_Resonance, KST2_1430_Resonance,

KST_1680_Resonance, NR );

hydra::Vector4R B0(D_MASS, 0.0, 0.0, 0.0);

{

double M2_12 = (kaon + pion1).mass2();

double M2_13 = (kaon + pion2).mass2();

double M2_23 = (pion1+ pion2).mass2();

});

TH3D* Dalitz_Resonances = new TH3D("Dalitz_Resonances",

"Dalitz - Toy Data -;"

"M^{2}(K^{-} #pi_{1}^{+}) [GeV^{2}/c^{4}];"

"M^{2}(K^{-} #pi_{2}^{+}) [GeV^{2}/c^{4}];"

"M^{2}(#pi_{1}^{+} #pi_{2}^{+}) [GeV^{2}/c^{4}]",

100, pow(K_MASS + PI_MASS,2), pow(D_MASS - PI_MASS,2),

100, pow(PI_MASS + PI_MASS,2), pow(D_MASS - K_MASS,2));

TH3D* Dalitz_Fit = new TH3D("Dalitz_Fit",

"Dalitz - Fit -;"

"M^{2}(K^{-} #pi_{1}^{+}) [GeV^{2}/c^{4}];"

"M^{2}(K^{-} #pi_{2}^{+}) [GeV^{2}/c^{4}];"

"M^{2}(#pi_{1}^{+} #pi_{2}^{+}) [GeV^{2}/c^{4}]",

100, pow(K_MASS + PI_MASS,2), pow(D_MASS - PI_MASS,2),

100, pow(PI_MASS + PI_MASS,2), pow(D_MASS - K_MASS,2));

TH2D* Normalization = new TH2D("normalization",

"Model PDF Normalization;Norm;Error",

200, 275.0, 305.0,

200, 0.58, 0.64);

TH3D* KST800_HIST , *KST892_HIST, *KST1425_HIST, *KST1430_HIST,

*KST1680_HIST, *NR_HIST ;

double KST800_FF, KST892_FF, KST1425_FF, KST1430_FF, KST1680_FF, NR_FF;

{

std::cout << std::endl;

std::cout << "======================================" << std::endl;

std::cout << "======= 1 - GENERATE TOY-DATA ========" << std::endl;

std::cout << "======================================" << std::endl;

auto start = std::chrono::high_resolution_clock::now();

auto end = std::chrono::high_resolution_clock::now();

std::chrono::duration<double, std::milli> elapsed = end - start;

std::cout << std::endl;

std::cout << "----------------- Device ----------------"<< std::endl;

std::cout << "| D+ -> K- pi+ pi+" << std::endl;

std::cout << "| Number of events :"<< toy_data.size() << std::endl;

std::cout << "| Time (ms) :"<< elapsed.count() << std::endl;

std::cout << "-----------------------------------------"<< std::endl;

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

}

{

std::cout << std::endl;

std::cout << "======================================" << std::endl;

std::cout << "========= 2 - PLOT TOY-DATA ==========" << std::endl;

std::cout << "======================================" << std::endl;

std::cout << std::endl << std::endl;

std::cout << "<======= [Daliz variables] { ( MSq_12, MSq_13, MSq_23) } =======>"<< std::endl;

for( size_t i=0; i<10; i++ )

std::cout << dalitz_variables[i] << std::endl;

{100,100,100},

{pow(K_MASS + PI_MASS,2), pow(K_MASS + PI_MASS,2), pow(PI_MASS + PI_MASS,2)},

{pow(D_MASS - PI_MASS,2), pow(D_MASS - PI_MASS ,2), pow(D_MASS - K_MASS,2)}

};

auto start = std::chrono::high_resolution_clock::now();

Hist_Dalitz.Fill( dalitz_variables );

auto end = std::chrono::high_resolution_clock::now();

std::chrono::duration<double, std::milli> elapsed = end - start;

std::cout << std::endl;

std::cout << "----------------- Device ----------------"<< std::endl;

std::cout << "| Sparse histogram fill" << std::endl;

std::cout << "| Number of events :"<< nentries << std::endl;

std::cout << "| Time (ms) :"<< elapsed.count() << std::endl;

std::cout << "-----------------------------------------"<< std::endl;

std::cout << std::endl;

std::cout << "Filling a ROOT Histogram... " << std::endl;

for(auto entry : Hist_Dalitz)

{

size_t bin = hydra::get<0>(entry);

double content = hydra::get<1>(entry);

unsigned int bins[3];

Hist_Dalitz.GetIndexes(bin, bins);

Dalitz_Resonances->SetBinContent(bins[0]+1, bins[1]+1, bins[2]+1, content);

}

{

std::cout << std::endl;

std::cout << "======================================" << std::endl;

std::cout << "=============== 3 - FIT ==============" << std::endl;

std::cout << "======================================" << std::endl;

std::cout << std::endl << std::endl;

std::cout << "-----------------------------------------"<<std::endl;

std::cout <<"| Initial PDF Norm: "<< Model_PDF.GetNorm() << "̣ +/- " << Model_PDF.GetNormError() << std::endl;

std::cout << "-----------------------------------------"<<std::endl;

particles.end());

ROOT::Minuit2::MnPrint::SetGlobalLevel(3);

std::cout<<fcn.GetParameters().GetMnState()<<std::endl;

auto start_d = std::chrono::high_resolution_clock::now();

FunctionMinimum minimum_d = FunctionMinimum( migrad_d(5000,250) );

auto end_d = std::chrono::high_resolution_clock::now();

std::chrono::duration<double, std::milli> elapsed_d = end_d - start_d;

std::cout << "-----------------------------------------"<<std::endl;

std::cout << "| [Migrad] Time (ms) ="<< elapsed_d.count() <<std::endl;

std::cout << "-----------------------------------------"<<std::endl;

std::cout<<"minimum: "<<minimum_d<<std::endl;

fcn.GetParameters().UpdateParameters(minimum_d);

fcn.GetPDF().GetFunctor().PrintRegisteredParameters();

auto start = std::chrono::high_resolution_clock::now();

phsp.Generate(B0, fit_data);

auto end = std::chrono::high_resolution_clock::now();

std::chrono::duration<double, std::milli> elapsed = end - start;

std::cout << std::endl;

std::cout << "----------- Device (fit data) -----------"<< std::endl;

std::cout << "| D+ -> K- pi+ pi+" << std::endl;

std::cout << "| Number of events :"<< nentries << std::endl;

std::cout << "| Time (ms) :"<< elapsed.count() << std::endl;

std::cout << "-----------------------------------------"<< std::endl;

auto dalitz_weights_fit = fit_data | fit_data.GetEventWeightFunctor(fcn.GetPDF().GetFunctor());

std::cout << std::endl;

std::cout << "<======= [Daliz variables - fit] { weight : ( MSq_12, MSq_13, MSq_23) } =======>"<< std::endl;

for( size_t i=0; i<10; i++ )

std::cout << dalitz_weights_fit[i] << " : "<< dalitz_variables_fit[i] << std::endl;

{100,100,100},

{pow(K_MASS + PI_MASS,2), pow(K_MASS + PI_MASS,2), pow(PI_MASS + PI_MASS,2)},

{pow(D_MASS - PI_MASS,2), pow(D_MASS - PI_MASS ,2), pow(D_MASS - K_MASS,2)}

};

start = std::chrono::high_resolution_clock::now();

Hist_Dalitz.Fill( dalitz_variables_fit, dalitz_weights_fit);

end = std::chrono::high_resolution_clock::now();

std::cout << std::endl;

std::cout << "----------------- Device ----------------"<< std::endl;

std::cout << "| Sparse histogram fill" << std::endl;

std::cout << "| Number of events :"<< nentries << std::endl;

std::cout << "| Time (ms) :"<< elapsed.count() << std::endl;

std::cout << "-----------------------------------------"<< std::endl;

auto Opt_Model = fcn.GetPDF().GetFunctor();

auto KST800 = fcn.GetPDF().GetFunctor().GetFunctor(_1);

auto KST892 = fcn.GetPDF().GetFunctor().GetFunctor(_2);

auto KST1425 = fcn.GetPDF().GetFunctor().GetFunctor(_3);

auto KST1430 = fcn.GetPDF().GetFunctor().GetFunctor(_4);

auto KST1680 = fcn.GetPDF().GetFunctor().GetFunctor(_5);

auto NR = fcn.GetPDF().GetFunctor().GetFunctor(_6);

KST800_FF = fit_fraction(KST800 , Opt_Model, {D_MASS, K_MASS, PI_MASS}, nentries);

KST892_FF = fit_fraction(KST892 , Opt_Model, {D_MASS, K_MASS, PI_MASS}, nentries);

KST1425_FF = fit_fraction(KST1425, Opt_Model, {D_MASS, K_MASS, PI_MASS}, nentries);

KST1430_FF = fit_fraction(KST1430, Opt_Model, {D_MASS, K_MASS, PI_MASS}, nentries);

KST1680_FF = fit_fraction(KST1680, Opt_Model, {D_MASS, K_MASS, PI_MASS}, nentries);

std::cout << "KST800_FF :" << KST800_FF << std::endl;

std::cout << "KST892_FF :" << KST892_FF << std::endl;

std::cout << "KST1425_FF :" << KST1425_FF << std::endl;

std::cout << "KST1430_FF :" << KST1430_FF << std::endl;

std::cout << "KST1680_FF :" << KST1680_FF << std::endl;

std::cout << "NR_FF :" << NR_FF << std::endl;

std::cout << "Sum :"

<< KST800_FF +

KST892_FF +

KST1425_FF +

KST1430_FF +

KST1680_FF +

NR_FF << std::endl;

{

std::vector integrals;

std::vector integrals_error;

for(auto x: fcn.GetPDF().GetNormCache() ){

integrals.push_back(x.second.first);

integrals_error.push_back(x.second.second);

}

auto integral_bounds = std::minmax_element(integrals.begin(),

integrals.end());

auto integral_error_bounds = std::minmax_element(integrals_error.begin(),

integrals_error.end());

Normalization->GetXaxis()->SetLimits(*integral_bounds.first, *integral_bounds.second);

Normalization->GetYaxis()->SetLimits(*integral_error_bounds.first, *integral_error_bounds.second);

for(auto x: fcn.GetPDF().GetNormCache() ){

Normalization->Fill(x.second.first, x.second.second );

}

std::cout << "Filling a ROOT Histogram... " << std::endl;

for(auto entry : Hist_Dalitz)

{

size_t bin = hydra::get<0>(entry);

double content = hydra::get<1>(entry);

unsigned int bins[3];

Hist_Dalitz.GetIndexes(bin, bins);

Dalitz_Fit->SetBinContent(bins[0]+1, bins[1]+1, bins[2]+1, content);

}

Dalitz_Fit->Scale(Dalitz_Resonances->Integral()/Dalitz_Fit->Integral() );

KST800_HIST->Scale( KST800_FF*Dalitz_Fit->Integral()/KST800_HIST->Integral() );

KST892_HIST->Scale( KST892_FF*Dalitz_Fit->Integral()/KST892_HIST->Integral() );

KST1425_HIST->Scale( KST1425_FF*Dalitz_Fit->Integral()/KST1425_HIST->Integral() );

KST1430_HIST->Scale( KST1430_FF*Dalitz_Fit->Integral()/KST1430_HIST->Integral() );

KST1680_HIST->Scale( KST1680_FF*Dalitz_Fit->Integral()/KST1680_HIST->Integral() );

NR_HIST->Scale( NR_FF*Dalitz_Fit->Integral()/NR_HIST->Integral() );

TH1* hist2D=0;

TCanvas* canvas_3= new TCanvas("canvas_3", "Dataset", 500, 500);

hist2D = Dalitz_Resonances->Project3D("yz");

hist2D->SetTitle("");

hist2D->Draw("colz");

canvas_3->SaveAs("plots/dalitz/Dataset1.pdf");

TCanvas* canvas_4= new TCanvas("canvas_4", "Dataset", 500, 500);

hist2D = Dalitz_Resonances->Project3D("xy");

hist2D->SetTitle("");

hist2D->Draw("colz");

canvas_4->SaveAs("plots/dalitz/Dataset2.pdf");

TCanvas* canvas_5= new TCanvas("canvas_5", "Fit", 500, 500);

hist2D = Dalitz_Fit->Project3D("yz");

hist2D->SetTitle("");

hist2D->SetStats(0);

hist2D->Draw("colz");

canvas_5->SaveAs("plots/dalitz/FitResult1.pdf");

TCanvas* canvas_6= new TCanvas("canvas_4", "Phase-space FLAT", 500, 500);

hist2D = Dalitz_Fit->Project3D("xy");

hist2D->SetTitle("");

hist2D->SetStats(0);

hist2D->Draw("colz");

canvas_6->SaveAs("plots/dalitz/FitResult2.pdf");

TH1* hist=0;

const char* axis =0;

auto KST800_Color = kViolet-5;

auto KST892_Color = kBlue;

auto KST1425_Color = kGreen-2;

auto KST1430_Color = kOrange-3;

auto KST1680_Color = kYellow-2;

auto NR_Color = kBlack;

double X1NDC = 0.262458;

double Y1NDC = 0.127544;

double X2NDC = 0.687708;

double Y2NDC = 0.35;

axis = "x";

TCanvas* canvas_x= new TCanvas("canvas_x", "", 600, 750);

auto legend_x = new TLegend( X1NDC, Y1NDC, X2NDC, Y2NDC);

legend_x->SetEntrySeparation(0.3);

legend_x->SetNColumns(2);

legend_x->SetBorderSize(0);

hist= Dalitz_Fit->Project3D(axis)->DrawCopy("hist");

hist->SetLineColor(2);

hist->SetLineWidth(2);

hist->SetMinimum(0.001);

hist->SetStats(0);

hist->SetTitle("");

legend_x->AddEntry(hist,"Fit","l");

hist= Dalitz_Resonances->Project3D(axis)->DrawCopy("e0same");

hist->SetMarkerStyle(8);

hist->SetMarkerSize(0.6);

hist->SetStats(0);

legend_x->AddEntry(hist,"Data","lep");

hist = KST800_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST800_Color);

hist->SetLineWidth(2);

legend_x->AddEntry(hist,"#kappa","l");

hist = KST892_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST892_Color);

hist->SetLineWidth(2);

legend_x->AddEntry(hist,"K*(892)","l");

hist = KST1680_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1680_Color);

hist->SetLineWidth(2);

legend_x->AddEntry(hist,"K_{1}(1680)","l");

hist = KST1425_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1425_Color);

hist->SetLineWidth(2);

legend_x->AddEntry(hist,"K*_{0}(1425)","l");

hist = KST1430_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1430_Color);

hist->SetLineWidth(2);

legend_x->AddEntry(hist,"K*_{2}(1430)","l");

hist = NR_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(NR_Color);

hist->SetLineStyle(5);

hist->SetLineWidth(2);

legend_x->AddEntry(hist,"NR","l");

canvas_x->SaveAs("plots/dalitz/Proj_X.pdf");

canvas_x->SetLogy(1);

legend_x->Draw();

canvas_x->SaveAs("plots/dalitz/Proj_LogX.pdf");

axis = "y";

TCanvas* canvas_y=new TCanvas("canvas_y", "", 600, 750);

auto legend_y = new TLegend( X1NDC, Y1NDC, X2NDC, Y2NDC);

legend_y->SetEntrySeparation(0.3);

legend_y->SetNColumns(2);

legend_y->SetBorderSize(0);

hist= Dalitz_Fit->Project3D(axis)->DrawCopy("hist");

hist->SetLineColor(2);

hist->SetLineWidth(2);

hist->SetMinimum(0.001);

hist->SetStats(0);

hist->SetTitle("");

legend_y->AddEntry(hist,"Fit","l");

hist= Dalitz_Resonances->Project3D(axis)->DrawCopy("e0same");

hist->SetMarkerStyle(8);

hist->SetMarkerSize(0.6);

hist->SetStats(0);

legend_y->AddEntry(hist,"Data","lep");

hist = KST800_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST800_Color);

hist->SetLineWidth(2);

legend_y->AddEntry(hist,"#kappa","l");

hist = KST892_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST892_Color);

hist->SetLineWidth(2);

legend_y->AddEntry(hist,"K*(892)","l");

hist = KST1680_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1680_Color);

hist->SetLineWidth(2);

legend_y->AddEntry(hist,"K_{1}(1680)","l");

hist = KST1425_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1425_Color);

hist->SetLineWidth(2);

legend_y->AddEntry(hist,"K*_{0}(1425)","l");

hist = KST1430_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1430_Color);

hist->SetLineWidth(2);

legend_y->AddEntry(hist,"K*_{2}(1430)","l");

hist = NR_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(NR_Color);

hist->SetLineStyle(5);

hist->SetLineWidth(2);

legend_y->AddEntry(hist,"NR","l");

canvas_y->SaveAs("plots/dalitz/Proj_Y.pdf");

canvas_y->SetLogy(1);

legend_y->Draw();

canvas_y->SaveAs("plots/dalitz/Proj_LogY.pdf");

axis = "z";

TCanvas* canvas_z= new TCanvas("canvas_z", "", 600, 750);

auto legend_z = new TLegend( X1NDC, Y1NDC, X2NDC, Y2NDC);

legend_z->SetEntrySeparation(0.3);

legend_z->SetNColumns(2);

legend_z->SetBorderSize(0);

hist= Dalitz_Fit->Project3D(axis)->DrawCopy("hist");

hist->SetLineColor(2);

hist->SetLineWidth(2);

hist->SetMinimum(0.001);

hist->SetStats(0);

hist->SetTitle("");

legend_z->AddEntry(hist,"Fit","l");

hist= Dalitz_Resonances->Project3D(axis)->DrawCopy("e0same");

hist->SetMarkerStyle(8);

hist->SetMarkerSize(0.6);

hist->SetStats(0);

legend_z->AddEntry(hist,"Data","lep");

hist = KST800_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST800_Color);

hist->SetLineWidth(2);

legend_z->AddEntry(hist,"#kappa","l");

hist = KST892_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST892_Color);

hist->SetLineWidth(2);

legend_z->AddEntry(hist,"K*(892)","l");

hist = KST1680_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1680_Color);

hist->SetLineWidth(2);

legend_z->AddEntry(hist,"K_{1}(1680)","l");

hist = KST1425_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1425_Color);

hist->SetLineWidth(2);

legend_z->AddEntry(hist,"K*_{0}(1425)","l");

hist = KST1430_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(KST1430_Color);

hist->SetLineWidth(2);

legend_z->AddEntry(hist,"K*_{2}(1430)","l");

hist = NR_HIST->Project3D(axis)->DrawCopy("histCsame");

hist->SetLineColor(NR_Color);

hist->SetLineStyle(5);

hist->SetLineWidth(2);

legend_z->AddEntry(hist,"NR","l");

canvas_z->SaveAs("plots/dalitz/Proj_Z.pdf");

canvas_z->SetLogy(1);

legend_z->Draw();

canvas_z->SaveAs("plots/dalitz/Proj_LogZ.pdf");

TCanvas* canvas_7=new TCanvas("canvas_7", "Normalization", 500, 500);

Normalization->Draw("colz");

}

template<typename Backend, typename Model, typename Container >

size_t generate_dataset(Backend const& system, Model const& model, std::array<double, 3> const& masses, Container& decays, size_t nevents, size_t bunch_size)

{

const double D_MASS = masses[0];

const double K_MASS = masses[1];

const double PI_MASS = masses[2];

hydra::Vector4R D(D_MASS, 0.0, 0.0, 0.0);

do {

phsp.SetSeed(S());

phsp.Generate(D, _data );

auto sample = _data.Unweight(model, -1, S());

decays.insert(decays.end(), sample.begin(), sample.end());

} while(decays.size()<nevents );

decays.erase(decays.begin()+nevents, decays.end());

return decays.size();

}

template<typename Amplitude, typename Model>

{

const double D_MASS = masses[0];

const double K_MASS = masses[1];

const double PI_MASS = masses[2];

hydra::Vector4R D(D_MASS, 0.0, 0.0, 0.0);

});

return amp_int.first/model_int.first;

}

template

{

const double D_MASS = masses[0];

const double K_MASS = masses[1];

const double PI_MASS = masses[2];

TH3D* Component = new TH3D(name,

";"

"M^{2}(K^{-} #pi^{+}) [GeV^{2}/c^{4}];"

"M^{2}(#pi^{+} #pi^{+}) [GeV^{2}/c^{4}]",

100, pow(K_MASS + PI_MASS,2), pow(D_MASS - PI_MASS,2),

100, pow(PI_MASS + PI_MASS,2), pow(D_MASS - K_MASS,2));

hydra::Vector4R D(D_MASS, 0.0, 0.0, 0.0);

{

double M2_12 = (kaon + pion1).mass2();

double M2_13 = (kaon + pion2).mass2();

double M2_23 = (pion1+ pion2).mass2();

});

phsp.Generate(D, events);

auto dalitz_weights = events | events.GetEventWeightFunctor(functor);

{100,100,100},

{pow(K_MASS + PI_MASS,2), pow(K_MASS + PI_MASS,2), pow(PI_MASS + PI_MASS,2)},

{pow(D_MASS - PI_MASS,2), pow(D_MASS - PI_MASS ,2), pow(D_MASS - K_MASS,2)}

};

Hist_Component.Fill( dalitz_variables, dalitz_weights );

for(auto entry : Hist_Component){

size_t bin = hydra::get<0>(entry);

double content = hydra::get<1>(entry);

unsigned int bins[3];

Hist_Component.GetIndexes(bin, bins);

Component->SetBinContent(bins[0]+1, bins[1]+1, bins[2]+1, content);

}

return Component;

}