#include <iostream>
#include <cmath>
#include "TMinuit.h"


using namespace std;
const double B_f = 0.24;
const double alpha = 0.0072973525643;
const double alpha_f = 1.0; // for lepton case
//
const double mtop=173.;
const double s=129600;//250000;//
const double beta_val = sqrt((1 - (4 * pow(mtop,2))/s));

// Define global vectors

const double mt =173.;
const double Mw =80.15;
const double r =(pow(Mw,2))/(pow(mt,2));
const double W =pow((1-r),2) * (1+(2*r));
const double B = (1-beta_val)/(1+beta_val);


// Define the SM factors
const double SinthetaW = 0.23; // sin^2 θW
const double CosthetaW = 1.0 - SinthetaW;


const double Ay = (4.0 / 3.0) * sqrt(SinthetaW);
const double By = 0.0;
const double Az = (1.0 / (2.0 * sqrt(CosthetaW))) * (1.0 - (8.0 / 3.0) * SinthetaW);
const double Bz = (1.0 / (2.0 * sqrt(CosthetaW)));
const double ve = -1 + (4 * SinthetaW);


// Non SM Contributions ::
double Dgamma=0.05;

// Z-mass, d, d_prime
const double Mz = 91.0;
const double d = (s / (s - pow(Mz, 2))) * (1 / (16 * SinthetaW * CosthetaW));
const double d_p = s / (4 * sqrt(SinthetaW) * sqrt(CosthetaW) * (s - pow(Mz,2)));
//const double d_p = d * 4*sqrt(SinthetaW)* sqrt(CosthetaW);
const double C = 1 / (4 * SinthetaW);


std::vector<double> x_values;
std::vector<double> costheta_values;


//C(DV,Az):

double C_Dv_AZ(double C, double ve, double d_p, double Az, double Ay){
         return( - 2 * C * ((ve*d_p*Ay) - (1+pow(ve,2)) *pow(d_p,2)*Az));
}

// C(DVA,Az):
double C_DVA_AZ(double C, double ve, double d_p, double Bz){
      return (C* ((1+pow(ve,2)) *pow(d_p,2) *Bz)) ;
}


//C(EV,AZ):
double C_EV_AZ(double C, double ve, double d_p, double  Az,double Ay){
       return (2* C*(d_p *Ay - (2*ve) *pow(d_p,2) *Az));
}

//C(EVA,AZ) :
double C_EVA_AZ(double C,double ve,double Bz,double d_p) {

       return (- C * ( 2*(ve)) * pow(d_p,2)*Bz);
}


/////////////////// The depend ones ://////////////////////////////////////
double C_DA_BZ(double C,double ve,double d_p,double Bz){

        return (2* C_DVA_AZ(C,ve,d_p,Bz));
}


double C_DVA_BZ(double C,double ve,double d_p,double Az,double Ay) {

 return (0.5 * C_Dv_AZ(C,ve,d_p,Az,Ay));

}


double C_DVA_AY(double C,double ve,double d_p,double Bz) {

   return (-C* ve * d_p *Bz);
}

double C_DV_AY(double C,double ve,double d_p,double Az,double Ay){

        return ( 2*C*(Ay - (ve*d_p *Az)));
}



double C_EVA_AY(double Bz,double C,double d_p){

        return ( C* d_p* Bz);
}


double C_EV_AY(double C,double d_p, double Az){

        return (   2* C* d_p *Az);

}


double C_EA_BZ(double C,double ve,double Bz,double d_p){

return (2 * C_EVA_AZ(C,ve,Bz,d_p));

}

double C_EVA_BZ(double C,double ve,double d_p,double Bz,double Ay){
  return (0.5 *C_EV_AZ(C,ve,d_p,Az,Ay));
}


//Gammas :
double C_DA_BY(double C,double ve,double d_p,double Bz){

        return (2*C * C_DVA_AY(C,ve, d_p,Bz));
}

double C_DVA_BY(double C,double ve,double d_p,double Az,double Ay) {

 return (0.5 *C * C_DV_AY(C,ve,d_p,Az,Ay));

}

double C_EA_BY(double C,double ve,double Bz,double d_p){

return (2*C * C_EVA_AY(Bz,C,d_p));

}


double C_EVA_BY(double C,double ve,double d_p,double Bz){

  return (C*0.5 * C_EV_AY(C,d_p,Az) );
}


// fro Gamma now :
// C(EVA,Agamma) :
double C_EVA_Ay(double Bz, double C, double d_p) {

return (C* d_p * Bz);

}


// C(DVA,Ay) :
double C_DVA_Ay(double C,double d_p, double Bz){
 return ( - C * d_p * Bz);

}

//C(DV,Ay):
double C_DV_Ay(double C,double ve,double Az,double Ay)  {

 return (2*C *(Ay -(ve * d_p* Az)));

}


// C(Ev,Ay):
double C_Ev_Ay(double C,double d_p,double Az){

return (2*C *  d_p *Az);
}


// C(F1,DV) :
//C(F1:Dv) = −C(DV :Av)/2
double C_F1_DZ(double C,double ve,double d_p,double Az,double Ay){
    return (-0.5 * C_Dv_AZ(C,ve,d_p,Az,Ay));
}

//C(F4:Dv) = −C(EVA:Av)
double C_F4_DZ(double C,double ve,double Bz,double d_p){
   return (- C_EVA_AZ(C,ve,Bz,d_p));

}


double C_F1_DY(double C,double ve,double d_p,double Az,double Ay){
    return (-0.5 * C_DV_AY(C,ve,d_p,Az,Ay));
}

//C(F4:Dv) = −C(EVA:Av)
double C_F4_DY(double C,double ve,double Bz,double d_p){
   return (- C_EVA_AY(Bz,C,d_p));

}



//For FC_Z :
/// C(G1,CZ):
double C_G1_CZ(double C,double ve,double d_p,double Az,double Ay){
   return (0.5*C_Dv_AZ(C,ve,d_p,Az,Ay));
}

// C(G2,CZ):
double C_G2_CZ(double C,double ve,double d_p,double Az,double Ay ){
//C_EV_AZ(double C, double ve, double d_p, double  Az,double Ay)
        return (0.5* C_EV_AZ(C,ve,d_p,Az,Ay) );

}

// G(G3,CZ):
double C_G3_CZ(double C,double ve,double d_p,double Bz){
   return (C_DVA_AZ(C,ve,d_p,Bz));
}


//FOPR GAMMA :
double C_G1_CY(double C,double ve,double d_p,double Az,double Ay){
   return (0.5*C_DV_AY(C,ve,d_p,Az,Ay));
}

// C(G2,CZ):
double C_G2_CY(double C,double ve,double d_p,double Az ){
//C_EV_AZ(double C, double ve, double d_p, double  Az,double Ay)
        return (0.5* C_EV_AY(C,d_p,Az) );

}


// G(G3,CZ):
double C_G3_CY(double C,double ve,double d_p,double Bz){
   return (C_DVA_AY(C, ve,d_p,Bz));
}


double H1(double beta_val, double omega_plus, double omega_minus, double x_f, double W) {
    return (((1 - pow(beta_val, 2)) / x_f) *
            (1 / (2 * W * beta_val)) *
            ((pow(omega_plus, 2) * (3 - 2 * (omega_plus)) - ( pow(omega_minus, 2) * (3 - 2 * (omega_minus)) ))));


}
double H2(double beta_val, double omega_plus, double omega_minus, double x_f, double W) {
    return (1 / (4 * beta_val * W) *
            (1 + beta_val) *
            pow((1 - beta_val), 2) *
            (1 / pow(x_f, 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))));
}


double F(double W,double beta_val,double omega_plus, double omega_minus){

        return (((1+beta_val)/beta_val) * (3/(2*W)) * ( pow(omega_plus,2) - pow(omega_minus,2)  ) );
}


double G(double W, double beta_val, double omega_plus, double omega_minus, double x_f) {
    return (F(W, beta_val, omega_plus, omega_minus) +
            (3 * std::pow((1 + beta_val), 2) / (beta_val * W) *
            x_f * (omega_plus + std::log(std::abs(1 - omega_plus)) - ( omega_minus  + std::log(std::abs(1 - omega_minus))))));
}


double DV(double Ay, double Az, double d_p, double ve,double C) {
    double term1 =( Ay * Ay - 2 * Ay * Az * ve * d_p + Az * Az * (1 + pow(ve,2)) * d_p * d_p);
    double term2 =2*(Ay- Az * ve* d_p)*Dgamma;
    double term3 =-2*( Ay * ve* d_p -Az *(1+pow(ve,2)) *pow(d_p,2))*Dgamma;
    return (C*(term1+term2+term3));
}


double DA(double Bz,double C,double ve, double d_p) {
    double term1 = (Bz * Bz * (1 + pow(ve,2)) * pow(d_p,2));
    double term2 = -2 *Bz * ve * d_p *Dgamma;
    double term3 = 2 * Bz *(1+pow(ve,2))*pow(d_p,2)*Dgamma;
    return (C*(term1+term2+term3));
}


double DVA(double Ay, double Az,double C,double Bz, double d_p, double ve) {
    double term1 = (-Ay * Bz * d_p *ve  + Az * Bz * (1 + pow(ve,2)) * d_p * d_p - Bz*ve*d_p*Dgamma);
    double term2 = (Ay -ve * d_p *Az)* Dgamma + Bz*(1+ pow(ve,2)) *pow(d_p,2)*Dgamma;
    double term3 = - (Ay * ve * d_p -Az *(1+pow(ve,2)) *pow(d_p,2))*Dgamma;

    return (C*(  term1+term2+term3));
  //  return term1;
}

double EV(double Ay, double Az, double C,double d_p,double ve) {
    double term1 = 2*C * (Ay * Az * d_p - pow(Az,2)*ve * d_p * d_p  + Az *d_p *Dgamma + (Ay *d_p -2*Az *ve * pow(d_p,2))*Dgamma );
    return term1;
}


double EA(double Bz, double ve, double C,double d_p) {
    double term1 = 2*C *(-Bz * Bz * pow(d_p,2)*ve  + Bz *d_p *Dgamma - 2*Bz * ve *pow(d_p,2)*Dgamma) ;
    return term1;
}

double EVA(double Ay, double ve,double Az, double Bz, double C,double d_p) {

    double term1 = 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) ;
    return term1;
}




//sm :
double Function_Theta(double x_f, double thetaCosine) {

double omega_minus=0,omega_plus=0;
double omega=0;

if (r < B) {

     if(x_f >= B*r && x_f <  r) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- x_f;
    }
    else if ( x_f >= B && x_f <= 1) {
         omega_minus = 0;
         omega_plus = 1- x_f;
    }


}

    // Terms for C0
    double C0_term1 = 0.5 * (3 - pow(beta_val, 2)) * DV(Ay, Az, d_p, ve, C) -
                      0.5 * (1 - (3 * pow(beta_val, 2))) * DA(Bz, C, ve, d_p);

    double C0_term2 = alpha_f * ((1 - pow(beta_val, 2)) * DVA(Ay, Az, C, Bz, d_p, ve));
    double C0_term3 = 2 * alpha_f * (DVA(Ay, Az, C, Bz, d_p, ve) * G(W,beta_val,omega_plus,omega_minus,x_f));
    double C0_term4 = 0.5*( (DV(Ay, Az, d_p, ve, C) + DA(Bz, C, ve, d_p) +
                       2 * alpha_f * (DVA(Ay, Az, C, Bz, d_p, ve)))) * (2*H1(beta_val,omega_plus,omega_minus, x_f, W)-H2(beta_val,omega_plus,omega_minus, x_f, W));

    // First Coefficient C0
    double C0 = ( (C0_term1 + C0_term2) * F(beta_val, W, omega_plus,omega_minus)) +  C0_term3 + C0_term4;

    // C1 Decomposition
    double C1_term1 = 2 * (2 * EVA(Ay, ve, Az, Bz, C, d_p) +
                            alpha_f * ((1 - pow(beta_val, 2)) * EA(Bz, ve, C, d_p))) * F(beta_val,W, omega_plus,omega_minus);
    double C1_term2 =  2* alpha_f * (EV(Ay, Az, d_p, ve, C) + EA(Bz, ve, C, d_p)) * G(W,beta_val,omega_plus,omega_minus, x_f);
//}
    double C1_term3 = -2 * (2 * EVA(Ay, ve, Az, Bz, C, d_p) +
                             alpha_f * (EV(Ay, Az, d_p, ve, C) + EA(Bz, ve, C, d_p))) * H1(beta_val,omega_plus,omega_minus, x_f, W);

    // Second Coefficient C1
    double C1 = C1_term1 + C1_term2 + C1_term3;

    // C2 Decomposition
    double C2_term1 = 0.5*( (3 - pow(beta_val, 2)) *
                       (DV(Ay, Az, d_p, ve, C) + DA(Bz, C, ve, d_p)) +
                       6* alpha_f * ((1 - pow(beta_val, 2)) * DVA(Ay, Az, C, Bz, d_p, ve)));
    double C2_term2 = alpha_f * 2 * DVA(Ay, Az, C, Bz, d_p, ve) * G(W,beta_val,omega_plus,omega_minus, x_f);

    double C2_term3 = (-3 / 2) *( (DV(Ay, Az, d_p, ve, C) + DA(Bz, C, ve, d_p) + 2*alpha_f* DVA(Ay, Az, C, Bz, d_p, ve)))*(2*H1(beta_val,omega_plus,omega_minus, x_f, W) -H2(beta_val,omega_plus,omega_minus, x_f, W) ) ;

    // Third Coefficient C2
    double C2 = (C2_term1 * F(beta_val,W, omega_plus,omega_minus)) +
                C2_term2 + C2_term3;


    // print the co
    return ( ( C0 +(C1*thetaCosine)+C2*pow(thetaCosine,2)));

}



double F_AZ(double x_f, double thetacosine){

double omega_minus=0,omega_plus=0;
double omega=0;

if (r < B) {

    if(x_f >= B*r && x_f <  r) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- x_f;
    }
     else if ( x_f >= B && x_f <= 1) {
         omega_minus = 0;
         omega_plus = 1- x_f;
    }

}

double term1 = (0.5 *(3 - pow(beta_val,2)) * C_Dv_AZ(C,ve,d_p,Az,Ay) * F(W,beta_val,omega_plus,omega_minus) + 2* C_DVA_AZ(C,ve,d_p,Bz) * G(W,beta_val,omega_plus,omega_minus,x_f));

double term2 = -( (1 - pow(beta_val,2)) * C_DVA_AZ(C,ve,d_p,Bz) * F(W,beta_val,omega_plus,omega_minus) - 0.5 * ( C_Dv_AZ(C,ve,d_p,Az,Ay)  + 2*C_DVA_AZ(C,ve,d_p,Bz) * (2*H1(beta_val,omega_plus,omega_minus,x_f,W)- H2(beta_val,omega_plus,omega_minus,x_f,W))));

double term3 = (2 *( C_EV_AZ(C,ve,d_p,Az,Ay) *( G(W,beta_val,omega_plus,omega_minus,x_f) - H1(beta_val,omega_plus,omega_minus,x_f, W)) + 2*C_EVA_AZ(C,ve,Bz,d_p) *(F(W,beta_val,omega_plus,omega_minus) -H1(beta_val,omega_plus,omega_minus,x_f, W))));


      return ( term1 *(1 + pow(thetacosine,2))  + term2*(1 - (3 *pow(thetacosine,2)) ) + term3*thetacosine);
 }



double F_Ay(double x_f, double thetacosine){

double omega_minus=0.0,omega_plus=0.0;
double omega=0.0;

if (r < B) {
    if(x_f >= B*r && x_f <  r) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- x_f;
    }
    else if ( x_f >= B && x_f <= 1) {
         omega_minus = 0;
         omega_plus = 1- x_f;
    }
}





double term1 =( 0.5 *(3-pow(beta_val,2)) * C_DV_AY(C,ve,d_p,Az,Ay) * F(W,beta_val,omega_plus,omega_minus) + 2* C_DVA_AY(C,ve,d_p,Bz) *G(W,beta_val,omega_plus,omega_minus,x_f));

double term2 = - ( (1-pow(beta_val,2))  *C_DVA_AY(C,ve,d_p,Bz) * F(W,beta_val,omega_plus,omega_minus) - 0.5 * (C_DV_AY(C,ve,d_p,Az,Ay) + 2*C_DVA_AY(C,ve,d_p,Bz) * (2*H1(beta_val,omega_plus,omega_minus,x_f,W) - H2(beta_val,omega_plus,omega_minus,x_f,W))));

double term3 = 2 * (C_EV_AY(C,d_p,Az) *( G(W,beta_val,omega_plus,omega_minus,x_f) - H1(beta_val,omega_plus,omega_minus,x_f, W)) + 2*C_EVA_AY(Bz,C,d_p)   *(F(W,beta_val,omega_plus,omega_minus) -H1(beta_val,omega_plus,omega_minus,x_f, W)));



    return ( term1*(1+pow(thetacosine,2)) + (term3*thetacosine) + term2*(1 - (3 * pow(thetacosine,2))) );

  }



// FB Z:
double F_BZ(double x_f, double thetacosine){

double omega_minus=0,omega_plus=0;
double omega=0;
if (r < B) {

    if(x_f >= B*r && x_f <  r) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- x_f;
    }
    else if ( x_f >= B && x_f <= 1) {
         omega_minus = 0;
         omega_plus = 1- x_f;
    }
}


    double ter1 = 0.5 *  (pow(beta_val,2)) * C_DA_BZ(C,ve,d_p,Bz) *F(W,beta_val,omega_plus,omega_minus) *(3-pow(thetacosine,2)) +2*(C_DVA_BZ(C,ve,d_p,Az,Ay) *G(W,beta_val,omega_plus,omega_minus,x_f));


   //C_DVA_AZ(double C, double ve, double d_p, double Bz)

    double ter2 = -0.5 *(( C_DA_BZ(C,ve,d_p,Bz) + 2*(1-pow(beta_val,2)) * C_DVA_AZ(C,ve,d_p, Bz) ) *F(W,beta_val,omega_plus,omega_minus)  -(C_DA_BZ(C,ve,d_p,Bz) + 2*C_DVA_BZ(C,ve,d_p,Az,Ay)) *(2*H1(beta_val,omega_plus,omega_minus,x_f, W) - H2(beta_val,omega_plus,omega_minus,x_f,W)));

    double ter3 = (  2*(((1-pow(beta_val,2)) * C_EA_BZ(C,ve,Bz,d_p)+ 2*  C_EVA_BZ(C,ve,d_p,Bz,Ay)) * F(W,beta_val,omega_plus,omega_minus) + C_EA_BZ(C,ve,Bz,d_p) *G(W,beta_val,omega_plus,omega_minus,x_f)  - (C_EA_BZ(C,ve,Bz,d_p) + 2*C_EVA_BZ(C,ve,d_p,Bz,Ay)) *H1(beta_val,omega_plus,omega_minus,x_f, W)));

    return ( ter1*(1+pow(thetacosine,2)) +ter2*(1-(3*pow(thetacosine,2)))+(ter3*thetacosine));



}



// FB Y:
double F_BY(double x_f, double thetacosine){

double w_minus,w_plus;
double omega;

if (r < B) {
    if(x_f >= B*r && x_f <  r) {
         w_minus = (1-x_f)/B;
         w_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         w_minus = (1-x_f)/B;
         w_plus = 1- x_f;
    }
    else if ( x_f >= B && x_f <= 1) {
         w_minus = 0;
         w_plus = 1- x_f;
    }
}

 omega = w_plus - w_minus;


//C_DVA_BY(double C,double ve,double d_p,double Az,double Ay)

//C_EA_BY(double C,double ve,double Bz,double d_p)
//C_EVA_BY(double C,double ve,double d_p,double Bz)
//C_DA_BY(double C,double ve,double d_p,double Bz)
//
 double ter1 = 0.5 *  (pow(beta_val,2) * C_DA_BY(C,ve,d_p,Bz) *F(B_f,W,beta_val,omega) *(3-pow(thetacosine,2))) +2*(C_DVA_BY(C, ve,d_p, Az,Ay)*G(W,beta_val,omega,x_f,B_f)) * (1+pow(thetacosine,2));

  double ter2 = -0.5 *(( C_DA_BY(C,ve,d_p,Bz) +  2*(1-pow(beta_val,2) * C_DVA_BY(C, ve,d_p, Az,Ay) ))*F(B_f,W,beta_val,omega)  - (C_DA_BY(C,ve, d_p, Bz)+2*C_DVA_BY(C, ve,d_p, Az,Ay) *(2*H1(B_f,beta_val,omega,x_f, W) - H2( B_f,beta_val,omega,x_f,W)))) *(1-3*pow(thetacosine,2));

  double ter3 = (  2*( (1-pow(beta_val,2) *C_EA_BY(C,ve,Bz, d_p) + 2*C_EVA_BY(C,ve,d_p, Bz))* F(B_f,W,beta_val,omega) + C_EA_BY(C,ve,Bz, d_p)*G(W,beta_val,omega,x_f,B_f)  - (C_EA_BY(C,ve,Bz, d_p) + 2*C_EVA_BY(C,ve,d_p,Bz) *H1(B_f,beta_val,omega,x_f, W))) *thetacosine);

   return (ter1 + ter2 +ter3);
}




double F_DZ(double x_f, double thetacosine){

double omega_minus=0,omega_plus=0;
double omega=0;
if (r < B) {
    if(x_f >= B*r && x_f <  r) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- x_f;
    }
    else if ( x_f >= B && x_f <= 1) {
         omega_minus = 0;
         omega_plus = 1- x_f;
    }
}



   double term1 = C_F1_DZ(C,ve,d_p,Az,Ay) *(F(W,beta_val,omega_plus,omega_minus) - H1(beta_val,omega_plus,omega_minus,x_f, W)) * (1-3*pow(thetacosine,2)) - C_F1_DZ(C,ve,d_p,Az, Ay)*G(W,beta_val,omega_plus,omega_minus,x_f) *(1+pow(thetacosine,2));
   double term2 = -2*C_F4_DZ(C,ve,Bz,d_p) *(F(W,beta_val,omega_plus,omega_minus)  + G(W,beta_val,omega_plus,omega_minus,x_f) - H1(beta_val,omega_plus,omega_minus,x_f, W));

      return ( term1+(term2*thetacosine));


   }

double F_DY(double x_f, double thetacosine){

double omega_minus=0,omega_plus=0;
double omega=0;


if (r < B) {
    if(x_f >= B*r && x_f <  r) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- r;
     }
     else if ( x_f >= r && x_f < B) {
         omega_minus = 1-x_f/B;
         omega_plus = 1- x_f;
    }
    else if ( x_f >= B && x_f <= 1) {
         omega_minus = 0;
         omega_plus = 1- x_f;
    }
}

   double term1 = C_F1_DY(C,ve,d_p,Az,Ay) * ( F(W,beta_val,omega_plus,omega_minus) - H1(beta_val,omega_plus,omega_minus,x_f, W)) * (1-(3*pow(thetacosine,2))) - C_F1_DY(C,ve,d_p,Az,Ay)*G(W,beta_val,omega_plus,omega_minus,x_f) *(1+pow(thetacosine,2));
   double term2 = -2* C_F4_DY(C,ve,Bz,d_p)*(F(W,beta_val,omega_plus,omega_minus)  + G(W,beta_val,omega_plus,omega_minus,x_f) - H1(beta_val,omega_plus,omega_minus,x_f, W));

     return ( term1+(term2*thetacosine));



 }



void PreprocessData(TTree * tree) {    
       
       x_values.clear();
       costheta_values.clear();
    
    
    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);

    TLorentzVector lepton_Plus,Top_Plus;
    double w_minus,w_plus;

     double thetacosine;

     double PfoEnP5,PfoPxP5,PfoPyP5,PfoPzP5;

    tree->SetBranchAddress("PfoEnP5",&PfoEnP5);
    tree->SetBranchAddress("PfoPxP5",&PfoPxP5);
    tree->SetBranchAddress("PfoPyP5",&PfoPyP5);
    tree->SetBranchAddress("PfoPzP5",&PfoPzP5);


    for (int i = 0; i < tree->GetEntries(); i++) {
        tree->GetEntry(i);
        TLorentzVector lepton_Plus;

        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]);
            }
        }
                double costheta = lepton_Plus.CosTheta();
                double E_f = lepton_Plus.E();
                double x = (2 * (E_f/mtop) * std::sqrt((1 - beta_val)/(1 + beta_val)));
                x_values.push_back(x);
                costheta_values.push_back(costheta);  
    }

}

void negativeLogLikelihood(double &nll,double *params) {
  
   double gamma1 = params[0];
   double gamma2 = params[1];
   double gamma3 = params[2];
   double gamma4 = params[3];
   double gamma5 = params[4];
   double gamma6 = params[5];

   nll = 0.0;

    // nbr of events == N:
       double N=42000;
       double fact=0;
       for(int i=1;i <= N;i++){
          fact += (log(i));
      }

   for (size_t i = 0; i < x_values.size(); i++) {
	double p1 = F_DZ(x_values[i], costheta_values[i]);
        double p2 = F_BY(x_values[i], costheta_values[i]);
        double p3 = F_BZ(x_values[i], costheta_values[i]);
        double p4 = F_AZ(x_values[i], costheta_values[i]);
        double p5 = F_Ay(x_values[i], costheta_values[i]);
        double p6 = F_DY(x_values[i], costheta_values[i]);  
        // SM one :
	double p7 = Function_Theta(x_values[i], costheta_values[i]);

       	// Likelihood contribution 
        double likelihood_event = (1/8.7) * 3*B_f*(M_PI *beta_val*pow(alpha,2))/(2*s) *(gamma1 * p1 + gamma2 * p2 + gamma3 * p3 + gamma4 * p4 + gamma5 * p5 + gamma6 * p6 +p7);

      // Sum of negative logs
     if(likelihood_event > 0 ) {
        // Accumulate the NLL
        nll += -std::log(likelihood_event);
    }
  
 }

       // multiply by a constant :
       nll *= - N*std::log(8.7) * (fact) * 8.7;
     
}

double Comp_Matrix[6][6];

void RunMinimization() {

    const int nParams = 6;  
   
    TMinuit *minuit = new TMinuit(nParams);
    minuit->SetPrintLevel(1);

    // minimization function
      minuit->SetFCN([](int &npar, double *gin, double &f, double *par, int iflag) {
        double nll;
	negativeLogLikelihood(nll, par);
	f = nll;
    });

 /*
   */
    
    double stepSize = 0.0001;
    int ierflg;

    //starting value of gammas:
    double startGamma = 0.08; 

    // gamma1 to gamma6
    minuit->mnparm(0, "gamma1", startGamma, stepSize, 0.05, 0.9, ierflg);
    minuit->mnparm(1, "gamma2", startGamma, stepSize, 0.04, 0.7, ierflg);
    minuit->mnparm(2, "gamma3", startGamma, stepSize, 0.04, 0.5, ierflg);
    minuit->mnparm(3, "gamma4", startGamma, stepSize, 0.02, 0.6, ierflg);
    minuit->mnparm(4, "gamma5", startGamma, stepSize, 0.02, 0.4, ierflg);
    minuit->mnparm(5, "gamma6", startGamma, stepSize, 0.02, 0.9, ierflg);

    // Minimization
    double arglist[6];
    arglist[0] = 5000; 
    arglist[1] = 0.001;
    int ierflag;
    minuit->mnexcm("MINIMIZE", arglist, 2, ierflag);

  // Retrieve results
    double gamma_opt[6], gamma_err[6];
    for (int i = 0; i < 6; i++) {
        minuit->GetParameter(i, gamma_opt[i], gamma_err[i]);
    }

    std::cout << "Fit Results:\n";
    for (int i = 0; i < 6; i++) {
        std::cout << "gamma" << i << ": " << gamma_opt[i] << " ± " << gamma_err[i] << "\n";
    }

    delete minuit;
}

void TMinuit() {

   TFile *file = TFile::Open("data.root", "READ");
   TTree *tree = (TTree*)(file->Get("tree_data"));

   // run this on all events :
   PreprocessData(tree);

   RunMinimization();
}