#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <math.h>
#include "TFile.h"
#include "TPostScript.h"
#include <iomanip>
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
#include "TVector3.h"
#include "TMath.h"
#include "TGraphErrors.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TMarker.h"
#include "TPaveText.h"

#include "TMinuit.h"
#include "TLegend.h"
#include "TText.h"
#include "TMatrix.h"
#include "TMath.h"

using namespace std;

Double_t g_norm_raw[157][20]={0};
Double_t g_norm_data[157][20]={0};
Double_t g_norm_sim[157][20]={0};



Double_t acc[190][20][10][10];
Int_t gi;
Int_t bin_E ;
Int_t bin_T ;

struct Events
{
    double x_ad;
    double y_ad;
    double qval_event;
    double acceptance;
    Double_t weight;
    double decay_weight;
    double lambda;
};

Events event[157][20][10000];
Int_t counter_event[157][20]={0};

Events event_sim[157][20][10000];
Int_t counter_event_sim[157][20]={0};

#define SQR(a) ((a)*(a))
#define M_PI 3.14159

Int_t counter[20]={0};

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t par[], Int_t iflag);
int main() {

TLorentzVector MyProton, MyPip, MyPim, MyProton_Rest, My_Initial, MyPi0, MyPhoton, MM;
         TVector3 omega_pvec;
         Float_t m_p = 0.938;
float omega_p;
char parname[500];
char name[500],name1[500], name8[500],name9[500],name3[500],name2[500], name_sim[500];
float temp[50],temp8[1500],temp9[1500],temp2[1500];
float runnum, eventnum, top, beampol, realphi, cosrealtheta, cos_theta_p_cm, beamE, p_px, p_py, p_pz, p_E, pip_px, pip_py, pip_pz, pip_E, pim_px, pim_py, pim_pz, pim_E;
float cl_nop, cl_pi0, z, kin_beamE, kin_p_px, kin_p_py, kin_p_pz, kin_p_E, kin_pip_px, kin_pip_py, kin_pip_pz, kin_pip_E, kin_pim_px, kin_pim_py, kin_pim_pz, kin_pim_E, q_value, q_error;

Double_t norm_raw[158][20]={0};
Double_t norm_data[158][20]={0};

Double_t norm_sim[158][20]={0};

TFile *fsim = new TFile("/d/grid13/akbar/sdme/sdme_sim_2017_full.root");
TFile *fraw = new TFile("/d/grid13/akbar/sdme/sdme_raw_2017_full.root");
TFile *f = new TFile("/d/grid13/akbar/sdme/dataforsdme.root");
TFile *f11 = new TFile("/d/grid13/akbar/get_g11_sdme_final.root");

 
TFile g ("sdme_dec_new.root", "RECREATE", "Histograms from ntuples" ); 

TH2D *h2[158][20];
TH2D *hsim[158][20];
TH2D *hraw[158][20];



//
for (int i=0; i<158; i++)
for (int j=0; j<20; j++)
{
{
  float ener = 1.72+(.01*i);
  float costet = -1 + (.1*j);
  sprintf(name1,"energy_%f_costheta_%f",ener,costet); 
  
  
 hsim[i][j] = (TH2D*)fsim->Get(name1);
 hraw[i][j] = (TH2D*)fraw->Get(name1);
  
  hsim[i][j]->Divide(hraw[i][j]);
    
}
}
//

cerr<<"c"<<endl;

for (int ee=0; ee<158; ee++)
{
for (int ii=0; ii<20; ii++)
 {
   for (int jj=1; jj<=10; jj++)
    {
      for (int kk=1; kk<=10; kk++)
       {
         float xx = hsim[ee][ii]->GetBinContent(jj,kk);
         acc[ee][ii][jj-1][kk-1] = xx;
         
       }
    }
  }
}
//
for (int ee=0; ee<158; ee++)
{
for (int ii=0; ii<20; ii++)
 {
   for (int jj=1; jj<=10; jj++)
    {
      for (int kk=1; kk<=10; kk++)
       {
         float xx = hraw[ee][ii]->GetBinContent(jj,kk);
         norm_raw[ee][ii] = norm_raw[ee][ii] + xx;
         
       }
    }
  }
}

//

Int_t counter_data = 0;
cerr<<"Initializing the data :"<<endl; 
 ifstream inputfile;
    sprintf(name,"/d/grid13/akbar/crossomega/outQvalue_g12_Sit05_W_1720_3310.txt");


    inputfile.open(name,ios::in);
      while (!inputfile.eof()){
      for (int i = 0; i < 38 ; i++) inputfile >> temp[i]; 
      runnum = temp[0];
      eventnum = temp[1];
      top = temp[2];
      beampol = temp[3];
      realphi = temp[4];
      cosrealtheta = temp[5];
      cos_theta_p_cm = temp[6];
      beamE = temp[7] /1000.0;
      p_px = temp[8] /1000.0;
      p_py = temp[9] /1000.0;
      p_pz = temp[10] /1000.0;
      p_E = temp[11] /1000.0;
      pip_px = temp[12] /1000.0;
      pip_py = temp[13] /1000.0;
      pip_pz = temp[14] /1000.0;
      pip_E = temp[15] /1000.0;
      pim_px = temp[16] /1000.0;
      pim_py = temp[17] /1000.0;
      pim_pz = temp[18] /1000.0;
      pim_E = temp[19] /1000.0;
      cl_nop = temp[20];
      cl_pi0 = temp[21];
      z = temp[22];
      kin_beamE = temp[23] /1000.0;
      kin_p_px = temp[24] /1000.0;
      kin_p_py = temp[25] /1000.0;
      kin_p_pz = temp[26] /1000.0;
      kin_p_E = temp[27] /1000.0;
      kin_pip_px = temp[28] /1000.0;
      kin_pip_py = temp[29] /1000.0;
      kin_pip_pz = temp[30] /1000.0;
      kin_pip_E = temp[31] /1000.0;
      kin_pim_px = temp[32] /1000.0;
      kin_pim_py = temp[33] /1000.0;
      kin_pim_pz = temp[34] /1000.0;
      kin_pim_E = temp[35] /1000.0;
      q_value = temp[36];
      q_error = temp[37];

      
      

      if(runnum != 56520 && runnum != 56544 && runnum != 56559 && runnum != 56619 && runnum != 56637 && runnum != 56585 && runnum > 56519)
      {     
                    
	      
                     
	            

                          MyProton.SetPxPyPzE( kin_p_px,kin_p_py, kin_p_pz, kin_p_E );
                          MyPip.SetPxPyPzE( kin_pip_px,kin_pip_py, kin_pip_pz, kin_pip_E );
                          MyPim.SetPxPyPzE( kin_pim_px,kin_pim_py, kin_pim_pz, kin_pim_E);
                          MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
                          My_Initial.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE+0.938);
                         MyPhoton.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE);
                          MyPi0 = My_Initial - MyProton - MyPip - MyPim;
                          MM = MyPip + MyPim + MyPi0;
			  omega_pvec = MM.Vect();
			  omega_p = omega_pvec.Mag();



			     //Method 2 (Chris) -
			    TLorentzRotation cm_r;
			     TVector3 cm_v = My_Initial.BoostVector();
			     cm_r.Boost(-cm_v);
			     TLorentzVector omega_cm2 = cm_r * MM;
			     TLorentzVector pip_cm2 = cm_r * MyPip;
			     TLorentzVector pim_cm2 = cm_r * MyPim;
			     TLorentzVector pi0_cm2 = cm_r * MyPi0;
                            TLorentzVector pho_cm2 = cm_r * MyPhoton;

                            TVector3 zprime = (pho_cm2.Vect()).Unit();
                            TVector3 yprime = (pho_cm2.Vect().Cross(omega_cm2.Vect())).Unit();
                            TVector3 xprime = yprime.Cross(zprime);
			    double costheta_omega2 = omega_cm2.CosTheta();

			 

			     TLorentzRotation helo_r;
			     TVector3 helo_v = omega_cm2.BoostVector();
			     helo_r.Boost(-helo_v);

			     TLorentzVector pip_hel2 = helo_r * pip_cm2;
			     TLorentzVector pim_hel2 = helo_r * pim_cm2; 
			     TLorentzVector pi0_hel2 = helo_r * pi0_cm2;
                             TLorentzVector test_vector = helo_r * omega_cm2;

                             double test = (test_vector.Vect()).Mag();

                            TVector3 pion = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Unit();
                            double pionmag = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Mag();
                            double pionmag2 = pionmag*pionmag;
                            double costheta_ad = pion.Dot(zprime);
                            double theta_ad = acos(costheta_ad);

                            TVector3 zcrosspi = (zprime.Cross(pion)).Unit();

                            double cosphi_ad = yprime.Dot(zcrosspi);
			    double sinphi_ad = -1*xprime.Dot(zcrosspi);
                           double phi_ad = atan2(sinphi_ad,cosphi_ad);

//
double lambda2_uncor = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Mag2();

		  Float_t T_pip = pip_hel2.E() - pip_hel2.Mag(); // K.E. of pi+
		  Float_t T_pim = pim_hel2.E() - pim_hel2.Mag(); // K.E. of pi-
		  Float_t T_pi0 = pi0_hel2.E() - pi0_hel2.Mag(); // K.E. of pi0
	          Float_t Q = T_pip + T_pim + T_pi0; // total K.E. in omega rest frame
	     	  Float_t M = pip_hel2.Mag();
		  Float_t lambda2_max = Q*Q*((Q*Q/108.0)+(Q*M/9.0)+(M*M/3.0)); // lamdba is max when the 3 pions are
			                                                       // on the same plane and are 120 deg. apart from each other.

		  //Note: lambda max is not |pi+|^2 * |pi-|^2  (i.e., when they go perpendicular to each other, because of E-p conservation constraints with pi0)
		  Float_t lambda2 = lambda2_uncor/lambda2_max;
//
                         
                         

double W = TMath::Sqrt(2*kin_beamE*0.938+0.938*0.938);
int binenergy = ((W-1.72)/0.01);

if (binenergy < 157 && cl_pi0 > 0.01 )
{
int binteta = int ((costheta_omega2-(-1))/0.1) ;

counter_event[binenergy][binteta] = counter_event[binenergy][binteta] + 1;

norm_data[binenergy][binteta] = norm_data[binenergy][binteta] + q_value;

event[binenergy][binteta][counter_event[binenergy][binteta]].x_ad = phi_ad;

event[binenergy][binteta][counter_event[binenergy][binteta]].y_ad = theta_ad;
event[binenergy][binteta][counter_event[binenergy][binteta]].qval_event = q_value;

event[binenergy][binteta][counter_event[binenergy][binteta]].decay_weight =lambda2;

event[binenergy][binteta][counter_event[binenergy][binteta]].weight = q_value;



counter_data = counter_data+1;
cerr<<"Initialize :"<<" "<<counter_data<<" "<<q_value<<" "<<counter_event_sim[binenergy][binteta]<<" "<<norm_data[binenergy][binteta]<<endl;
//cerr<<endl;
}

else
{
  cerr<<"HIGHER ENERGY"<<endl;
}

}
}
//
cerr<<"finish initializing data"<<endl;



//

Int_t counter_data_sim = 0;
cerr<<"Initializing the sim :"<<endl; 
 ifstream inputfile_sim;
    sprintf(name_sim,"/d/grid13/akbar/g12_analysis_new/omega_mc_allmc_allcut_sept22_newtrigger.txt");


    inputfile_sim.open(name_sim,ios::in);
      while (!inputfile_sim.eof()){
      for (int i = 0; i < 36 ; i++) inputfile_sim >> temp[i]; 
      runnum = temp[0];
      eventnum = temp[1];
      top = temp[2];
      beampol = temp[3];
      realphi = temp[4];
      cosrealtheta = temp[5];
      cos_theta_p_cm = temp[6];
      beamE = temp[7] /1000.0;
      p_px = temp[8] /1000.0;
      p_py = temp[9] /1000.0;
      p_pz = temp[10] /1000.0;
      p_E = temp[11] /1000.0;
      pip_px = temp[12] /1000.0;
      pip_py = temp[13] /1000.0;
      pip_pz = temp[14] /1000.0;
      pip_E = temp[15] /1000.0;
      pim_px = temp[16] /1000.0;
      pim_py = temp[17] /1000.0;
      pim_pz = temp[18] /1000.0;
      pim_E = temp[19] /1000.0;
      cl_nop = temp[20];
      cl_pi0 = temp[21];
      z = temp[22];
      kin_beamE = temp[23] /1000.0;
      kin_p_px = temp[24] /1000.0;
      kin_p_py = temp[25] /1000.0;
      kin_p_pz = temp[26] /1000.0;
      kin_p_E = temp[27] /1000.0;
      kin_pip_px = temp[28] /1000.0;
      kin_pip_py = temp[29] /1000.0;
      kin_pip_pz = temp[30] /1000.0;
      kin_pip_E = temp[31] /1000.0;
      kin_pim_px = temp[32] /1000.0;
      kin_pim_py = temp[33] /1000.0;
      kin_pim_pz = temp[34] /1000.0;
      kin_pim_E = temp[35] /1000.0;
      q_value = /*temp[36]*/1.0;
      q_error = /*temp[37]*/0.0;

      
      

                 
	      
                     
	            

                          MyProton.SetPxPyPzE( kin_p_px,kin_p_py, kin_p_pz, kin_p_E );
                          MyPip.SetPxPyPzE( kin_pip_px,kin_pip_py, kin_pip_pz, kin_pip_E );
                          MyPim.SetPxPyPzE( kin_pim_px,kin_pim_py, kin_pim_pz, kin_pim_E);
                          MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
                          My_Initial.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE+0.938);
                         MyPhoton.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE);
                          MyPi0 = My_Initial - MyProton - MyPip - MyPim;
                          MM = MyPip + MyPim + MyPi0;
			  omega_pvec = MM.Vect();
			  omega_p = omega_pvec.Mag();

  

			     //Method 2 (Chris) -
			    TLorentzRotation cm_r;
			     TVector3 cm_v = My_Initial.BoostVector();
			     cm_r.Boost(-cm_v);
			     TLorentzVector omega_cm2 = cm_r * MM;
			     TLorentzVector pip_cm2 = cm_r * MyPip;
			     TLorentzVector pim_cm2 = cm_r * MyPim;
			     TLorentzVector pi0_cm2 = cm_r * MyPi0;
                            TLorentzVector pho_cm2 = cm_r * MyPhoton;

                            TVector3 zprime = (pho_cm2.Vect()).Unit();
                            TVector3 yprime = (pho_cm2.Vect().Cross(omega_cm2.Vect())).Unit();
                            TVector3 xprime = yprime.Cross(zprime);
			    double costheta_omega2 = omega_cm2.CosTheta();

			 

			     TLorentzRotation helo_r;
			     TVector3 helo_v = omega_cm2.BoostVector();
			     helo_r.Boost(-helo_v);

			     TLorentzVector pip_hel2 = helo_r * pip_cm2;
			     TLorentzVector pim_hel2 = helo_r * pim_cm2; 
			     TLorentzVector pi0_hel2 = helo_r * pi0_cm2;
                             TLorentzVector test_vector = helo_r * omega_cm2;

                             double test = (test_vector.Vect()).Mag();

                            TVector3 pion = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Unit();
                            double pionmag = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Mag();
                            double pionmag2 = pionmag*pionmag;
                            double costheta_ad = pion.Dot(zprime);
                            double theta_ad = acos(costheta_ad);  

                            TVector3 zcrosspi = (zprime.Cross(pion)).Unit();

                            double cosphi_ad = yprime.Dot(zcrosspi);
			    double sinphi_ad = -1*xprime.Dot(zcrosspi);
                           double phi_ad = atan2(sinphi_ad,cosphi_ad);

//
double lambda2_uncor = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Mag2();

		  Float_t T_pip = pip_hel2.E() - pip_hel2.Mag(); // K.E. of pi+
		  Float_t T_pim = pim_hel2.E() - pim_hel2.Mag(); // K.E. of pi-
		  Float_t T_pi0 = pi0_hel2.E() - pi0_hel2.Mag(); // K.E. of pi0
	          Float_t Q = T_pip + T_pim + T_pi0; // total K.E. in omega rest frame
	     	  Float_t M = pip_hel2.Mag();
		  Float_t lambda2_max = Q*Q*((Q*Q/108.0)+(Q*M/9.0)+(M*M/3.0)); // lamdba is max when the 3 pions are
			                                                       // on the same plane and are 120 deg. apart from each other.

		  //Note: lambda max is not |pi+|^2 * |pi-|^2  (i.e., when they go perpendicular to each other, because of E-p conservation constraints with pi0)
		  Float_t lambda2 = lambda2_uncor/lambda2_max;
//                     
                         

double W = TMath::Sqrt(2*kin_beamE*0.938+0.938*0.938);
int binenergy = ((W-1.72)/0.01);
//int binenergy = ((W-1.72)/0.04);

if (binenergy < 157 && cl_pi0 > 0.01 /*&& sect_status != 4 && sect_status == 2*/)
{
int binteta = int ((costheta_omega2-(-1))/0.1) ;

counter_event_sim[binenergy][binteta] = counter_event_sim[binenergy][binteta] + 1;

norm_sim[binenergy][binteta] = norm_sim[binenergy][binteta] + 1;

event_sim[binenergy][binteta][counter_event_sim[binenergy][binteta]].x_ad = phi_ad;

event_sim[binenergy][binteta][counter_event_sim[binenergy][binteta]].y_ad = theta_ad;
event_sim[binenergy][binteta][counter_event_sim[binenergy][binteta]].qval_event = q_value;

event_sim[binenergy][binteta][counter_event_sim[binenergy][binteta]].decay_weight = lambda2;

event_sim[binenergy][binteta][counter_event_sim[binenergy][binteta]].weight = q_value;




counter_data_sim = counter_data_sim+1;
cerr<<"Initialize SIM :"<<counter_data_sim<<" "<<counter_event_sim[binenergy][binteta]<<" "<<binenergy<<endl;
//cerr<<endl;
}

else
{
  cerr<<"simm HIGHER"<<endl;
}

}

//
cerr<<"finish initializing sim"<<endl;

double par0_init[157][20]={0.0025};
double par1_init[157][20]={0.0025};
double par2_init[157][20]={0.0025};

TH1D *g11_p0[157];
TH1D *g11_p1[157];
TH1D *g11_p2[157];



for (int ii=0; ii<157; ii++)
for (int jj=0; jj<20; jj++)
{
{
  g_norm_data[ii][jj] = norm_data[ii][jj];
  g_norm_sim[ii][jj] = norm_sim[ii][jj];
  //g_norm_raw[ii][jj] = norm_raw[ii][jj];
}
}

//
TH1D *sdme0[157];

cerr<<"A"<<endl;
for (int i=0; i<157; i++)
{
 float elow = 1.72+(.01*i);

sprintf(name1,"par0_sdme_E%f",elow);
sprintf(name2,"SDME_E%f_to_E%f",elow,elow+0.04);
sdme0[i] = new TH1D(name1,name2, 20, -1, 1);
/*sdme0[i]->SetMaximum(0.8);
sdme0[i]->SetMinimum(0.0);*/

if (i<110)
{
sprintf(parname,"par0_w_cross_W%f",elow);
g11_p0[i] = (TH1D*)f11->Get(parname);

  for (int j=0; j<20; j++)
 {
   par0_init[i][j] = g11_p0[i]->GetBinContent(j+1);
 }
}

}

TH1D *sdme1[157];


for (int i=0; i<157; i++)
{
 float elow = 1.72+(.01*i);

sprintf(name1,"par1_sdme_E%f",elow);
sprintf(name2,"SDME_E%f_to_E%f",elow,elow+0.04);
sdme1[i] = new TH1D(name1,name2, 20, -1, 1);
/*sdme1[i]->SetMaximum(0.2);
sdme1[i]->SetMinimum(-0.2);*/

if (i<110)
{
sprintf(parname,"par1_w_cross_W%f",elow);
g11_p1[i] = (TH1D*)f11->Get(parname);

  for (int j=0; j<20; j++)
 {
   par1_init[i][j] = g11_p1[i]->GetBinContent(j+1);
 }
}

}

TH1D *sdme2[157];
cerr<<"b"<<endl;

for (int i=0; i<157; i++)
{
 float elow = 1.72+(.01*i);

sprintf(name1,"par2_sdme_E%f",elow);
sprintf(name2,"SDME_E%f_to_E%f",elow,elow+0.04);
sdme2[i] = new TH1D(name1,name2, 20, -1, 1);
/*sdme2[i]->SetMaximum(0.2);
sdme2[i]->SetMinimum(-0.2);*/

if (i<110)
{
sprintf(parname,"par2_w_cross_W%f",elow);
g11_p2[i] = (TH1D*)f11->Get(parname);

  for (int j=0; j<20; j++)
 {
   par2_init[i][j] = g11_p2[i]->GetBinContent(j+1);
 }
}

}

TH1D *sdme3[157];
cerr<<"b"<<endl;

for (int i=0; i<157; i++)
{
 float elow = 1.72+(.01*i);

sprintf(name1,"par3_sdme_E%f",elow);
sprintf(name2,"SDME_E%f_to_E%f",elow,elow+0.01);
sdme3[i] = new TH1D(name1,name2, 20, -1, 1);


}





cerr<<"d"<<endl;
cerr<<"Start minuit fitting"<<endl;
cerr<<endl;
for (int ii=10; ii<157; ii++)
for (int jj=0; jj<20; jj++)
 {
 {
  bin_E = ii;
  bin_T = jj;
  if (counter_event[ii][jj] >0)
  {
  cerr<<"start bin :  "<<ii<<endl;
  cerr<<endl;
  int flag;
  gMinuit = new TMinuit(3); 
  //gMinuit = new TMinuit(4); 
  gMinuit -> SetPrintLevel(0); 
  gMinuit->SetFCN(fcn);
  
 //gMinuit->mnparm(0,"par_0",0.0025,/*0.0025*/0.001,0.0,0.0,flag);
  //gMinuit->mnparm(1,"par_1",0.0025,/*0.0025*/0.001,0.0,0.0,flag);
  //gMinuit->mnparm(2,"par_2",0.0025,/*0.0025*/0.001,0.0,0.0,flag);
 gMinuit->mnparm(0,"par_0",par0_init[ii][jj]-0.05,0.0025,0.0,0.0,flag); //initial, step, limit, limit
  gMinuit->mnparm(1,"par_1",par1_init[ii][jj]+0.05,0.0025,0.0,0.0,flag);
  gMinuit->mnparm(2,"par_2",par2_init[ii][jj]-0.05,0.0025,0.0,0.0,flag);
  //gMinuit->Migrad();
 Double_t arglist[10];
   Int_t ierflg = 0;
 arglist[0] = 5000;
   arglist[1] = 1.;
   //gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
   gMinuit->Migrad();

  /*gMinuit->mnparm(0,"par_0",0.0025,0.0025,0.0,0.0,flag);
  gMinuit->mnparm(1,"par_1",0.0025,0.0025,0.0,0.0,flag);
  gMinuit->mnparm(2,"par_2",0.0025,0.0025,0.0,0.0,flag);
  //gMinuit->mnparm(3,"par_3",0.0025,0.0025,0.0001,0.0,flag); //additional parameter
  //gMinuit->mnparm(3,"par_3",0.0025,0.0025,0.0,0.0,flag); //additional parameter
  gMinuit->Migrad();*/
  //gMinuit->mnexcm("HESSE");

  double parreturn[3], parreturnError[3];
  //double parreturn[4], parreturnError[4];
  
  for(int parn = 0; parn < 3/*4*/; parn++)
   {
    gMinuit->GetParameter(parn,parreturn[parn],parreturnError[parn]);
   }

  sdme0[bin_E]->SetBinContent(bin_T+1,parreturn[0]);
  sdme0[bin_E]->SetBinError(bin_T+1,parreturnError[0]);

  sdme1[bin_E]->SetBinContent(bin_T+1,parreturn[1]);
  sdme1[bin_E]->SetBinError(bin_T+1,parreturnError[1]);

  sdme2[bin_E]->SetBinContent(bin_T+1,parreturn[2]);
  sdme2[bin_E]->SetBinError(bin_T+1,parreturnError[2]);

  /*sdme3[bin_E]->SetBinContent(bin_T+1,parreturn[3]);
  sdme3[bin_E]->SetBinError(bin_T+1,parreturnError[3]);*/

  cerr<<endl;
  cerr<<parreturn[0]<<"   "<<parreturnError[0]<<endl;
  cerr<<parreturn[1]<<"   "<<parreturnError[1]<<endl;
  cerr<<parreturn[2]<<"   "<<parreturnError[2]<<endl;
  //cerr<<"NORM :"<<parreturn[3]<<"   "<<norm_data[bin_E][bin_T]<<" "<<norm_sim[bin_E][bin_T]<<" "<<norm_raw[bin_E][bin_T]<<endl;
  cerr<<endl;
  cerr<<"finish bin :  "<<gi<<endl;
  cerr<<endl;
  }

 }
 }

cerr<<"WRITING"<<endl;
for (int i=10; i<157; i++)

{
sdme0[i]->Write();
sdme1[i]->Write();
sdme2[i]->Write();
//sdme3[i]->Write();
}

}



void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t par[], Int_t iflag)
{

  Double_t nll = 0.0;
  Double_t nll_sim = 0.0;
  Double_t nll_sim_prob = 0.0;
  Double_t prob = 0.0;
  Double_t prob_sim = 0.0;
  Double_t delta_N = 0.0;
  Double_t weight_data_total = 0.0;

  Double_t ecm = 1.72 + 0.01*bin_E + 0.005;
  Double_t ws = ecm*ecm;
  Double_t fs = TMath::Sqrt((ws-1.72*1.72)*(ws-0.16*0.16))/(4*ws*TMath::Power(2*M_PI,4));

  Double_t N_data = g_norm_data[bin_E][bin_T];
  Double_t N_raw = g_norm_raw[bin_E][bin_T];
  Double_t N_sim = g_norm_sim[bin_E][bin_T];
  Double_t N_event = counter_event[bin_E][bin_T];
  
  Double_t norm_total = (4*M_PI*N_data/(N_sim*0.89));

  

  for (Int_t i =1; i<=counter_event[bin_E][bin_T]; i++)
  {
    double x = event[bin_E][bin_T][i].x_ad;
    double y = event[bin_E][bin_T][i].y_ad;
    Double_t w = event[bin_E][bin_T][i].weight;
    Double_t w_decay = event[bin_E][bin_T][i].decay_weight;
    weight_data_total = weight_data_total + w*w_decay;

    
        
    prob = 0.239*(0.5*(1-par[0])+0.5*(3*par[0]-1)*(cos(y)*cos(y))-par[1]*(sin(y)*sin(y))*cos(2*x)-1.414*par[2]*sin(2*y)*cos(x));
    

    if(prob>0) nll = nll + w*TMath::Log(prob);
  }

  for (Int_t i =1; i<=counter_event_sim[bin_E][bin_T]; i++)
  {
    
    double x_sim = event_sim[bin_E][bin_T][i].x_ad;
    double y_sim = event_sim[bin_E][bin_T][i].y_ad;
    Double_t w_sim = event_sim[bin_E][bin_T][i].weight;
    Double_t w_sim_decay = event_sim[bin_E][bin_T][i].decay_weight;
    

   
      prob_sim = 0.239*(0.5*(1-par[0])+0.5*(3*par[0]-1)*(cos(y_sim)*cos(y_sim))-par[1]*(sin(y_sim)*sin(y_sim))*cos(2*x_sim)-1.414*par[2]*sin(2*y_sim)*cos(x_sim));
    
    

    if(prob_sim>0) nll_sim = nll_sim + w_sim_decay*(prob_sim);
    
  }


f = -2*(nll-norm_total*(nll_sim));

}