#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGenericPdf.h"
#include "RooArgList.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooFormulaVar.h"
#include <iostream>
#include <cmath>
#include "TMinuit.h"
#include <vector>
#include "TFile.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TFitResult.h"
#include "RooRealVar.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooClassFactory.h"
#include "TROOT.h"
 
using namespace RooFit;
using namespace std;


void Rootfit_Test() {

//variables globally :

RooRealVar B_f("B_f", "B_f", 0.24);
RooRealVar s("s", "s", 129600);
RooRealVar MP("MP", "Pi constant", TMath::Pi());
RooRealVar alpha("alpha", "alpha constant", 0.0072973525643);
RooRealVar mtop("mtop", "mtop", 173);
RooRealVar beta_val("beta_val", "beta_val constant", 0.318433);
RooRealVar MZ("MZ", "MZ", 91);
RooRealVar d_p("d_p", "d_p", 0.633434);
RooRealVar C("C", "C constant", 1.08);
RooRealVar Bz("Bz", "Bz constant", 0.324675);
RooRealVar Az("Az", "Az", 0.229703);
RooRealVar By("By", "By constant", 0);
RooRealVar Ay("Ay", "aY constant", -0.666667);
RooRealVar B("B", "B constant", 0.516952);
RooRealVar ve("ve", "ve constant", -0.08);
RooRealVar CosthetaW("CosthetaW", "CosthetaW constant", 0.77);
RooRealVar SinthetaW("SinthetaW", "SinthetaW constant", 0.23);
RooRealVar r("r", "r constant", 0.214107);
RooRealVar W("W", "W constant", 0.882105);
RooRealVar Dgamma("Dgamma", "Dgamma constant", 0.05);
RooRealVar x("x", "x", 0.2, 1);
RooRealVar costheta("costheta", "costheta", -0.98, 0.98);
RooRealVar N("N","N",224283);
RooFormulaVar fact("fact", "gamma(N + 1)", RooArgList(N));

RooRealVar costheta_val("costheta_val","costheta_val",-0.98, 0.98);
RooRealVar x_val("x_val","x_val",0.2,1);
RooRealVar E_f("E_f","E_f",13,130);


RooFormulaVar C_Dv_AZ(
    "C_Dv_AZ",
    "-2 * C * (ve * d_p * Ay - (1 + pow(ve, 2)) * pow(d_p, 2) * Az)", 
    RooArgList(C, ve, d_p, Az, Ay)
);


/*RooFormulaVar C_DVA_BZ(
    "C_DVA_BZ",
    "0.5 * C_Dv_AZ", 
    RooArgList(C_Dv_AZ)
    //RooArgList(C, ve, d_p, Az, Ay)
);*/


RooFormulaVar C_DVA_AZ(
  "C_DVA_AZ",
  "C * (1+ pow(ve,2)) * pow(d_p,2)*Bz",
   RooArgList(C,ve,d_p,Bz)
);


RooFormulaVar C_EV_AZ(
  "C_EV_AZ",
   "2 * C *(d_p *Ay - 2 * ve * pow(d_p, 2)  *Az )",
   RooArgList(C,ve,d_p,Az,Ay)  
  		
);

RooFormulaVar C_EVA_AZ(
    "C_EVA_AZ",
    "-C * (2 * ve * pow(d_p, 2) * Bz)",
    RooArgList(C, ve, d_p, Bz)
);

RooFormulaVar C_DA_BZ(
     "C_DA_BZ",
      "2* C_DVA_AZ",
      RooArgList(C_DVA_AZ)
);

/////////////////// The depend ones ://////////////////////////////////////
RooFormulaVar C_DVA_AY(
    "C_DVA_AY",
     "-C * ve * d_p * Bz",
     RooArgList(C,ve,d_p,Bz)    
		
);


RooFormulaVar C_DV_AY(
     "C_DV_AY",	
      "2*C*(Ay - ve * d_p *Az )",     
       RooArgList(C,Ay,ve,d_p,Az)
);


/////////// Complete here :
RooFormulaVar C_EV_AY(
   "C_EV_AY",
    "2* C* d_p * Az",
    RooArgList(C,Ay,d_p,Az)   
		
);


RooFormulaVar C_EVA_AY(
    "C_EVA_AY",
    " C* d_p * Bz", 
    RooArgList(C,d_p,Bz)   
		
);


RooFormulaVar C_EA_BZ(
    "C_EA_BZ",
    "2 * C_EVA_AZ",
    RooArgList(C_EVA_AZ) 
);

RooFormulaVar C_EVA_BZ(
    "C_EVA_BZ",
     "0.5* C_EV_AZ",
     RooArgList(C_EV_AZ)    
);

RooFormulaVar C_DA_BY(
     "C_DA_BY",
     "2* C_DVA_AY",
      RooArgList(C_DVA_AY)     
		
);

//Gammas :
RooFormulaVar C_DVA_BY(
      "C_DVA_BY",
       "1/2. * C_DV_AY",
        RooArgList(C_DV_AY)      
);

RooFormulaVar C_EA_BY(
   "C_EA_BY",
    "2 *C_EVA_AY",
    RooArgList(C_EVA_AY)   
  	
);

RooFormulaVar C_EVA_BY(
   "C_EVA_BY",
   "1/2. * C_EV_AY",
   RooArgList(C_EV_AY)

);

/*
RooFormulaVar C_DVA_BZ(
     "C_DVA_BZ",
      "0.5 * C_Dv_AZ",
      RooArgList(C_Dv_AZ)     

);
*/



// C(EVA,Agamma) :
RooFormulaVar C_EVA_Ay(
   "C_EVA_Ay",
    "C* d_p * Bz",
    RooArgList(C,d_p,Bz)   
);

RooFormulaVar C_DVA_Ay(
    "C_DVA_Ay",
    "-C * d_p * Bz",
     RooArgList(C,d_p,Bz)
	
);


RooFormulaVar C_DV_Ay(
    "C_DV_Ay",
    "2 * C * (Ay - ve * d_p * Az)",
    RooArgList(C, d_p, Az, ve, Ay)
);


/// F_i 
RooFormulaVar C_F1_DZ(
    "C_F1_DZ",
    "-0.5 * C_Dv_AZ",  
      RooArgList(C_Dv_AZ)
    );


RooFormulaVar C_F4_DZ(
    "C_F4_DZ",
    "- C_EVA_AZ",  
      RooArgList(C_EVA_AZ)
    );


RooFormulaVar C_F1_DY(
    "C_F1_DY",
    "-0.5 * C_DV_AY",
    RooArgList(C_DV_AY)   
		
);

RooFormulaVar C_F4_DY(
    "C_F4_DY",
     "-C_EVA_AY",
     RooArgList(C_EVA_AY)    
);


// Z :
RooFormulaVar C_G1_CZ(
      "C_G1_CZ",
      "0.5*C_Dv_AZ",
       RooArgList(C_Dv_AZ)
);

RooFormulaVar C_G2_CZ(
    "C_G2_CZ",
     "0.5 * C_EV_AZ",
      RooArgList(C_EV_AZ)    
);

RooFormulaVar C_G3_CZ(
  "C_G3_CZ",
  "C_DVA_AZ",
   RooArgList(C_DVA_AZ)  
		
);


RooFormulaVar C_DVA_BZ(
    "C_DVA_BZ",
    "1/2. * C_DVA_AZ",
    RooArgList(C_DVA_AZ)
);


///  Gammas :
RooFormulaVar C_G1_CY(
  "C_G1_CY",
   "0.5 * C_DV_AY",
    RooArgList(C_DV_AY)
);


RooFormulaVar C_G2_CY(
	"C_G2_CY",
	 "0.5 * C_EV_AY",
	  RooArgList(C_EV_AY)
	);

RooFormulaVar C_G3_CY(
    "C_G3_CY",
    "C_DVA_AY",
     RooArgList(C_DVA_AY)
	
);




// Omega PLus -------------------------
RooFormulaVar omega_plus(
    "omega_plus",
    "((r < B) ? ((x >= B*r && x < r) ? (1 - r) : "
    "              (x >= r && x < B) ? (1 - x) : "
    "              (x >= B) ? (1 - x) : 0) : 0)",
    RooArgList(r, B) 
);

// Omega Minus ------------------------
RooFormulaVar omega_minus(
    "omega_minus",
    "((r < B) ? ((x >= B*r && x < r) ? (1 - x/B) : "
    "              (x >= r && x < B) ? (1 - x/B) : "
    "              (x >= B) ? 0 : 0) : 0)",
    RooArgList(r, B) 
);


// F -----------------------------------------------------------------
RooFormulaVar F(
    "F",
    "((1 + beta_val)/beta_val) * (3/(2*W)) * (pow( omega_plus ,2) - pow( omega_minus ,2))",
    RooArgList(W, beta_val,omega_minus,omega_plus)
);


// H1 -----------------------------------------------------------------
RooFormulaVar H1(
    "H1",
    "((1 - pow(beta_val, 2)) / x) * (1 / (2 * W * beta_val)) * (pow(omega_plus, 2) * (3 - 2 * omega_minus) - pow(omega_minus, 2) * (3 - 2 * omega_minus))",
    RooArgList(W, beta_val, omega_minus, omega_plus)
);

// H2 -----------------------------------------------------------------
RooFormulaVar H2(
    "H2",
    "(1/(4*beta_val*W)) * (1 + beta_val) * pow(1 - beta_val,2) * (1/pow(x,2)) * "
    "(pow(omega_plus,2)*(6 - 8*omega_plus + 3*pow(omega_plus,2)) - "
    "pow(omega_minus,2)*(6 - 8*omega_minus + 3*pow(omega_minus,2)))",
    RooArgList(W, beta_val,omega_minus,omega_plus)
);

// G -----------------------------------------------------------------
RooFormulaVar G(
    "G",
    "F + (3*pow(1 + beta_val,2)/(beta_val*W)) * x * ("
    "(omega_plus + std::log(std::abs(1 - omega_plus)) - "
    "(omega_minus + std::log(std::abs(1 - omega_minus)))))",
    RooArgList(F, W, beta_val,omega_minus,omega_plus)  
);


// EVA., EV, DV ...
RooFormulaVar DV(
    "DV",
    "C * (Ay * Ay - 2 * Ay * Az * ve * ve + pow(Az, 2) * (1 + pow(ve, 2)) * d_p * d_p + pow(Az, 2) * (1 + pow(ve, 2)) * pow(d_p, 2) "
    " + 2 * (Ay - Az * ve * d_p) * Dgamma - 2 * (Ay * ve * d_p - Az * (1 + pow(ve, 2)) * pow(d_p, 2)) * Dgamma)",
    RooArgList(Ay, Az, ve, d_p, C, Dgamma)
);


RooFormulaVar DA(
        "DA",
        "C * (Bz * Bz * (1 + pow(ve, 2)) * pow(d_p, 2) - 2 * Bz * ve * d_p * Dgamma + 2 * Bz * (1 + pow(ve, 2)) * pow(d_p, 2) * Dgamma)",
        RooArgList(Bz, C, ve, d_p, Dgamma)
    );


RooFormulaVar DVA(
        "DVA",
        "C * ((-Ay * Bz * d_p * ve + Az * Bz * (1 + pow(ve, 2)) * pow(d_p, 2) - Bz * ve * d_p * Dgamma) "
        " + ((Ay - ve * d_p * Az) * Dgamma + Bz * (1 + pow(ve, 2)) * pow(d_p, 2) * Dgamma) "
        " - (Ay * ve * d_p - Az * (1 + pow(ve, 2)) * pow(d_p, 2)) * Dgamma)",
        RooArgList(Ay, Az, C, Bz, d_p, ve, Dgamma)
);

RooFormulaVar EV(
        "EV",
        "2 * C * (Ay * Az * d_p - pow(Az, 2) * ve * pow(d_p, 2) + Az * d_p * Dgamma + (Ay * d_p - 2 * Az * ve * pow(d_p, 2)) * Dgamma)",
        RooArgList(Ay, Az, C, d_p, ve, Dgamma)
    );


RooFormulaVar EA(
        "EA",
        "2 * C * (-Bz * Bz * pow(d_p, 2) * ve + Bz * d_p * Dgamma - 2 * Bz * ve * pow(d_p, 2) * Dgamma)",
        RooArgList(Bz, ve, C, d_p, Dgamma)
    );

RooFormulaVar EVA(
        "EVA",
        "C * (Ay * Bz * d_p - 2 * (Az * Bz * pow(d_p, 2) * ve) + Bz * d_p * Dgamma + "
        "(Az * d_p * Dgamma) - 2 * Bz * ve * pow(d_p, 2) * Dgamma + "
        "(Ay * d_p - 2 * Az * ve * pow(d_p, 2)) * Dgamma)",
        RooArgList(Ay, ve, Az, Bz, C, d_p, Dgamma)
);



// OMEGAS for Sigma Calculation :

RooFormulaVar Omega_Zero(
    "Omega_Zero",
    "DV - (1 - 2*pow(beta_val, 2) * DA - (2*(1 - 2*pow(beta_val, 2)) * DVA + DV + DA + 2*DVA * (1 - pow(beta_val, 2)) / (2 * beta_val) * log((1 - beta_val) / (1 + beta_val))))",
    RooArgList(DV, DA, beta_val, DVA)
);

RooFormulaVar Omega_One(
    "Omega_One",
 "4 * EVA + 2*(1 - 2*pow(beta_val,2))*EA - 2*EVA + EV + EA * (1 - pow(beta_val,2))/(2*beta_val) * log((1-beta_val)/(1+beta_val))",    
   RooArgList(EVA, beta_val, EA, EV)
);

RooFormulaVar Omega_Two(
    "Omega_Two",
    "(3 - 2 * pow(beta_val, 2)) * DV + DA + 3 * (2 * (1 - 2 * pow(beta_val, 2)) * DVA - 3 * (DV + DA + 2 * DVA * (1 - pow(beta_val, 2)) / (2 * beta_val) * log((1 - beta_val) / (1 + beta_val))))",
    RooArgList(DV, DA, DVA, beta_val)
);


// Observables
RooRealVar gamma1("gamma1", "gamma1", -0.9, 2);
RooRealVar gamma2("gamma2", "gamma2", -0.5, 1);
RooRealVar gamma3("gamma3", "gamma3", -0.9, 5);
RooRealVar gamma4("gamma4", "gamma4", -0.4, 3);
RooRealVar gamma5("gamma5", "gamma5", -0.5, 2);
RooRealVar gamma6("gamma6", "gamma6", -0.4, 3);
RooRealVar gamma7("gamma7", "gamma7", -0.3, 3.);
RooRealVar gamma8("gamma8", "gamma8", -0.2, 2);


    
RooArgList params(gamma1,gamma2, gamma2, gamma3, gamma4, gamma5, gamma6,gamma7, gamma8);


RooFormulaVar sigma_total(
    "sigma_total", 
    "(3*B_f*(MP *beta_val*pow(alpha,2))/(s))* ( Omega_Zero + (1./3.)*Omega_Two)", 
    RooArgList(MP, B_f, beta_val, alpha, s,Omega_Zero,Omega_Two) 
);


RooFormulaVar SM(
    "SM",
    "(0.5 * (3 - pow(beta_val, 2)) * DV - 0.5 * (1 - (3 * pow(beta_val, 2))) * DA) + "
    "((1 - pow(beta_val, 2)) * DVA * F + 2 * (DVA * G)) + 0.5 * (DV + DA + 2 * DVA) * (2 * H1 - H2) + "
    "2 * (2 * EVA * F + 2 * (EV + EA)) * G - 2 * (2 * (EVA + EV + EA)) * H1 + "
    "0.5 * ((3 - pow(beta_val, 2)) * (DV + DA) + 6 * ((1 - pow(beta_val, 2)) * DVA)) * F + "
    "2 * DVA * G - (-3. / 2.) * (DV + DA + 2 * DVA) * (2 * H1 - H2)",
    RooArgList(EVA,EV,EA,DV, DA, DVA, F, G, H1, H2, beta_val)
);


 RooFormulaVar FDZ(
    "FDZ",
    "C_F1_DZ * "
    "(F - "
    "H1) * (1 - 3 * pow(costheta, 2)) - "
    "C_F1_DZ * "
    "G * (1 + pow(costheta, 2)) - "
    "2 * C_F4_DZ * "
    "(F + "
    "G - "
    "H1)*costheta",
    RooArgList(costheta, H1,G,F,C_F1_DZ, C_F4_DZ
    )
);


RooFormulaVar FDY(
    "FDY",
    "C_F1_DY * (F - H1) * (1 - 3 * pow(costheta, 2)) - C_F1_DY * G * (1 + pow(costheta, 2)) - 2 * C_F4_DY * (F + G - H1) * costheta",
    RooArgList(costheta, H1, G, F, C_F4_DY, C_F1_DY)
);


RooFormulaVar FCZ(
    "FCZ",
    "(-pow(beta_val, 2) * C_G1_CZ * F * (1 + pow(costheta, 2)) + "
    "2 * C_G2_CZ * (F + G - H1 * costheta) - "
    "((C_G1_CZ + (2 - pow(beta_val, 2)) * C_G3_CZ) * F + "
    "C_G3_CZ * G - "
    "((2 * C_G1_CZ + 3 * C_G3_CZ) * H1 + "
    "(C_G1_CZ + C_G3_CZ) * H2)))",
    RooArgList( costheta, beta_val, C_G1_CZ, C_G2_CZ, C_G3_CZ, F, G, H1, H2)
);


RooFormulaVar FCY(
    "FCY",
    "(-pow(beta_val, 2) * C_G1_CY * F * (1 + pow(costheta, 2)) + "
    "2 * C_G2_CY * (F + G - H1 * costheta) - "
    "((C_G1_CY + (2 - pow(beta_val, 2)) * C_G3_CY) * F + "
    "C_G3_CY * G - "
    "((2 * C_G1_CY + 3 * C_G3_CY) * H1 + "
    "(C_G1_CY + C_G3_CY) * H2)))",
    RooArgList( costheta, beta_val, C_G1_CY, C_G2_CY, C_G3_CY, F, G, H1, H2)
);


RooFormulaVar FBZ(
   "FBZ",
  "0.5 *  ((pow(beta_val,2)) * C_DA_BZ *F *(3-pow(costheta,2)) +2*(C_DVA_BZ *G)) *(1+pow(costheta,2)) -0.5 *(( C_DA_BZ + 2*(1-pow(beta_val,2)) * C_DVA_BZ) *F  -(C_DA_BZ + 2*C_DVA_BZ) *(2*H1 - H2)) * (1-(3*pow(costheta,2))) + (2*(((1-pow(beta_val,2)) * C_EA_BZ+ 2*  C_EVA_BZ) * F+ C_EA_BZ *G - (C_EA_BZ + 2*C_EVA_BZ) *H1))*costheta",
  RooArgList(costheta,beta_val,C_EA_BZ,C_EVA_BZ,C_DA_BZ,C_DVA_BZ,F, G, H1, H2)
 );


RooFormulaVar FBY(
    "FBY",
    "0.5 * ((pow(beta_val, 2)) * C_DA_BY * F * (3 - pow(costheta, 2)) + 2 * (C_DVA_BY * G)) * (1 + pow(costheta, 2)) - 0.5 * ((C_DA_BY + 2 * (1 - pow(beta_val, 2)) * C_DVA_BY) * F - (C_DA_BY + 2 * C_DVA_BY) * (2 * H1 - H2)) * (1 - (3 * pow(costheta, 2))) + (2 * (((1 - pow(beta_val, 2)) * C_EA_BY + 2 * C_EVA_BY) * F + C_EA_BY * G - (C_EA_BY + 2 * C_EVA_BY) * H1)) * costheta", 
    RooArgList(costheta, beta_val, C_EA_BY, C_EVA_BY, C_DA_BY, C_DVA_BY, F, G, H1, H2)
);


RooFormulaVar FAZ(
    "FAZ",
    "(0.5 * (3 - pow(beta_val, 2)) * C_Dv_AZ * F + 2 * C_DVA_AZ * G) * (1 + pow(costheta, 2)) - "
    "((1 - pow(beta_val, 2)) * C_DVA_AZ * F - 0.5 * (C_Dv_AZ + 2 * C_DVA_AZ * (2 * H1 - H2))) * (1 - (3 * pow(costheta, 2))) + "
    "(2 * (C_EVA_AZ * G - H1) + 2 * C_EVA_AZ * (F - H1)) * costheta",
    RooArgList(costheta, beta_val, C_Dv_AZ, C_DVA_AZ, C_EVA_AZ, F, G, H1, H2)
);


RooFormulaVar FAY(
    "FAY",
     "(0.5 * (3 - pow(beta_val, 2)) * C_DV_AY * F + 2 * C_DVA_AY * G) * (1 + pow(costheta, 2)) - "
    "((1 - pow(beta_val, 2)) * C_DVA_AY * F - 0.5 * (C_DV_AY + 2 * C_DVA_AY * (2 * H1 - H2))) * (1 - (3 * pow(costheta, 2))) + "
    "(2 * (C_EVA_AY * G - H1) + 2 * C_EVA_AY * (F - H1)) * costheta",
     RooArgList(costheta, beta_val, C_DV_AY, C_DVA_AY, C_EVA_AY, C_EV_AY, F, G, H1, H2)
);

  TFile *file = TFile::Open("data.root");
  TTree *tree = (TTree*)file->Get("Merged_tree");
   
  double mc_En[100],mc_Px[100],mc_Pz[100],mc_Py[100];
  double mc_pdgid[100];
  tree->SetBranchAddress("mc_En[100]",mc_En);
  tree->SetBranchAddress("mc_Px[100]",mc_Px);
  tree->SetBranchAddress("mc_Py[100]",mc_Py);
  tree->SetBranchAddress("mc_Pz[100]",mc_Pz);
  tree->SetBranchAddress("mc_pdgid[100]",mc_pdgid);


    double PfoEnP5,PfoPxP5,PfoPyP5,PfoPzP5;
    tree->SetBranchAddress("PfoEnP5",&PfoEnP5);
    tree->SetBranchAddress("PfoPxP5",&PfoPxP5);
    tree->SetBranchAddress("PfoPyP5",&PfoPyP5);
    tree->SetBranchAddress("PfoPzP5",&PfoPzP5);
   
    RooDataSet data2("data2", "MC Tree", RooArgSet(x, costheta));
    TLorentzVector lepton_Plus;
    for (int i = 0; i < tree->GetEntries(); i++) {
         tree->GetEntry(i);
        // TLorentzVector lepton_Plus(PfoPxP5,PfoPyP5,PfoPzP5,PfoEnP5);
        for (int j = 0; j <100; j++) {
            if( mc_pdgid[j] == 13 || mc_pdgid[j] == 11 ) {
                  lepton_Plus.SetPxPyPzE(mc_Px[j], mc_Py[j], mc_Pz[j], mc_En[j]); }
            }

         costheta_val = lepton_Plus.Z() /
                     sqrt(lepton_Plus.X() * lepton_Plus.X() +
                          lepton_Plus.Y() * lepton_Plus.Y() +
                          lepton_Plus.Z() * lepton_Plus.Z());
         E_f = lepton_Plus.E(); 
         x_val = (2*(E_f.getVal()/mtop.getVal()) * sqrt((1-beta_val.getVal())/(1+beta_val.getVal())));

	/*costheta_val = lepton_Plus.CosTheta();
        E_f = lepton_Plus.E();
        x_val = (2 * (E_f.getVal() / mtop.getVal()) * std::sqrt((1 - beta_val.getVal()) / (1 + beta_val.getVal())));
        */
	if(x_val.getVal() >0 && x_val.getVal() <=1){
        x.setVal(x_val.getVal());
        costheta.setVal(costheta_val.getVal());
       
         data2.add(RooArgSet(x, costheta));
	}
   }


   // https://root.cern.ch/doc/master/rf104__classfactory_8C.html
    RooAbsPdf *genpdf = RooClassFactory::makePdfInstance("GenPdf", 
     //224283
     "pow(1000 * pow(10, -15) * sigma_total,2) * (1/fact) * exp(-1000 * pow(10, -15) * sigma_total) *(3 * B_f * MP * beta_val * pow(alpha, 2)) / (2 * s) * "
    "(gamma1 * FDZ + gamma2 * FBY + gamma3 * FBZ + gamma4 * FDY + gamma5 * FAZ + gamma6 * FAY + gamma7 * FCZ + SM+ gamma8 * FCY)",
    RooArgSet(x, SM, costheta, fact,gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7, gamma8, 
              MP, B_f, beta_val, alpha, s, sigma_total, FDZ, FBY, FBZ, FDY, FAZ, FAY, FCZ, FCY)); 

 // Toy DATa Generation :    
//std::unique_ptr<RooDataSet> data2{genpdf->generate(RooArgSet(x, costheta), 50000)}; 
 

genpdf->fitTo(data2,RooFit::PrintLevel(1));
//genpdf->graphVizTree("model.dot");

std::cout << "gamma1: " << gamma1.getVal() << std::endl;
std::cout << "gamma2: " << gamma2.getVal() << std::endl;
std::cout << "gamma3: " << gamma3.getVal() << std::endl;
std::cout << "gamma4: " << gamma4.getVal() << std::endl;
std::cout << "gamma5: " << gamma5.getVal() << std::endl;
std::cout << "gamma6: " << gamma6.getVal() << std::endl;
std::cout << "gamma7: " << gamma7.getVal() << std::endl;
std::cout << "gamma8: " << gamma8.getVal() << std::endl;
    

// Create frame: 
RooPlot* xframe = x.frame();
RooPlot* costhetaFrame = costheta.frame();

// Plot:
data2.plotOn(xframe);
genpdf->plotOn(xframe);

data2.plotOn(costhetaFrame);
genpdf->plotOn(costhetaFrame);


TCanvas* c = new TCanvas("c", "Fit Example 2D", 800, 600);
c->Divide(2, 1); 
c->cd(1);
xframe->Draw();
c->cd(2);
costhetaFrame->Draw();

c->SaveAs("fit_plot_2D.png");

}