#include <iostream>
#include <fstream>

gROOT->Reset();

// Setup file including parameters
#include "plotMk_3D_setup.h"

double RN = 1;

void plotMk_3D_All() {

  //###########################################################
  // Get normalization factor
  
  if (!useOneHist) {
    double mB = 5.2790;
    double mD = 1.8646;
    TRandom *rdm1 = new TRandom();
    TRandom2 *rdm2 = new TRandom2();
    TRandom *rdm3 = new TRandom();
    TRandom3 *rdm4 = new TRandom3();
    int N = 1000000;
    int N_old = 0;
    int N_new = 0;
    double RN = 0;
    for (Int_t i=0;i<N;i++) {
      double y1 = (rdm1->Rndm(i))*2; 
      double q2 = (rdm2->Rndm(i))*11.66;
      
      // Get Fplus
      double mass=mD;
      double mtb, mbb;
      double msd, mx,mb,nf,nfp; 
      double msq,bx2,mbx,mtx;
      double zji,cji,w,gammaji,chiji,betaji_fppfm;
      double rfppfm,rfpmfm,f3fppfm,f3fpmfm,fppfm,fpmfm,al,ai,rcji,f3; 
      double mqm,msb,bb2,mup,mum,bbx2,tm,wt,r2,betaji_fpmfm;
      
      double t=q2;
      
      msb=5.2;
      msd=0.33;
      bb2=0.431*0.431;
      mbb=5.31;
      nf = 4.0;
      
      msq=1.82;
      bx2=0.45*0.45;
      mbx=0.75*2.01+0.25*1.97;
      nfp = 3.0;
      
      mtb = msb + msd;
      mtx = msq + msd;
      mb=5.279;
      mx=mass;
      
      mup=1.0/(1.0/msq+1.0/msb);
      mum=1.0/(1.0/msq-1.0/msb);
      bbx2=0.5*(bb2+bx2);
      tm=(mb-mx)*(mb-mx);
      if ( t>tm ) t=0.99*tm;
      wt=1.0+(tm-t)/(2.0*mbb*mbx);
      
      mqm = 0.1;
      
      double PI = 3.1415;
      double lqcd2 = 0.04;
      double nflav = 4;
      double As0=0.6;
      if ( mqm > 0.6 ) {
	if ( mqm < 1.85 ) {
	  nflav = 3.0;}
	
	As0 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( mqm*mqm/lqcd2);
      }
      
      double As2=0.6;
      if ( msq > 0.6 ) {
	if ( msq < 1.85 ) {
	  nflav = 3.0;}
	
	As2 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( msq*msq/lqcd2);
      }
      r2=3.0/(4.0*msb*msq)+3*msd*msd/(2*mbb*mbx*bbx2) + 
	(16.0/(mbb*mbx*(33.0-2.0*nfp)))*
	log(As0/As2);
      
      f3 = sqrt(mtx/mtb)*pow(sqrt(bx2*bb2)/bbx2,1.5) /
	((1.0+r2*(tm-t)/12.0)*(1.0+r2*(tm-t)/12.0));
      
      w = 1.0 + (( tm - t ) / ( 2.0* mb * mx ));
      rcji = ( 1/sqrt(w*w -1 ))*log( w + sqrt( w*w -1 ));
      al = (8.0 / ( 33.0 - 2.0*nfp ))*(w*rcji -1.0 );
      ai = -1.0* ( 6.0/( 33.0 - 2.0*nf));  
      
      double As1=0.6;
      if ( msb > 0.6 ) {
	if ( msb < 1.85 ) {
	  nflav = 3.0;}
	
	As1 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( msb*msb/lqcd2);
      }
      
      
      cji = pow(( As1 / As2 ),ai);
      
      zji = msq / msb;
      
      gammaji = 2+((2.0*zji)/(1-zji))*log(zji);
      gammaji = -1.0*gammaji;
      
      chiji = -1.0 - ( gammaji / ( 1- zji ));
      betaji_fppfm = gammaji - (2.0/3.0)*chiji;
      betaji_fpmfm = gammaji + (2.0/3.0)*chiji;
      
      double As3=0.6;
      double m3=sqrt(msb*msq);
      
      if ( m3 > 0.6 ) {
	if ( msq < 1.85 ) {
	  nflav = 3.0;}
	
	As3 = 12.0*PI / ( 33.0 - 2.0*nflav) /
	log( m3*m3/lqcd2);
      }
      rfppfm = cji *(1.0 + betaji_fppfm*As3/PI);
      
      rfpmfm = cji *(1.0 + betaji_fpmfm*As3/PI);
      
      f3fppfm = f3*pow(( mbb / mtb ),-0.5)*pow((mbx/mtx),0.5);
      f3fpmfm = f3*pow(( mbb / mtb ),0.5)*pow((mbx/mtx),-0.5);
      fppfm = f3fppfm* rfppfm * ( 2.0 - 
				  ( ( mtx/msq)*(1- ( (msd*msq*bb2)
						     /(2.0*mup*mtx*bbx2)))));
      fpmfm = f3fpmfm* rfpmfm * ( mtb/msq) * ( 1 - 
					       ( ( msd*msq*bb2)/
						 ( 2.0*mup*mtx*bbx2)));
      
      double fpf = (fppfm + fpmfm)/2.0;
      
      // Get R
      double mcR = 2*sqrt(mB*mD) / (mB+mD);
      
      // Get old FF
      double oldFF = mcR * fpf;
      double oldFF2 = oldFF*oldFF;
      double rto1 = (mB*mB + mD*mD - t)/(2*mB*mD);
      double rto2 = rto1*rto1;
      double rto3 = rto2 - 1;
      double rto4 = sqrt(rto3*rto3*rto3);
      double oldInt =  rto4*oldFF2;
      if (y1 <= oldInt) N_old++;
      
      // Get new FF
      double y2 = (rdm1->Rndm(i))*2; 
      double trueW = (rdm4->Rndm(i))*0.592 + 1;
      
      if (Slope = "New") {
	double trueZ 
	  = (sqrt(trueW+1) - sqrt(2))/(sqrt(trueW+1) + sqrt(2));
	double N01_new = 8*rho2*trueZ;
	double N02_new = (51*rho2 - 10)*trueZ*trueZ;
	double N03_new = (252*rho2 - 84)*trueZ*trueZ*trueZ;
	double newFF = 1 - N01_new + N02_new - N03_new;
      } else {
	double newFF = 1 - rho2*(trueW - 1);
      }
      double newFF2 = newFF*newFF;
      double rto5 = trueW*trueW - 1;
      double rto6 = sqrt(rto5*rto5*rto5);
      double newInt = rto6*newFF2;
      if (y2 <= newInt) N_new++;
    }
    
    if (N_new > 0) {
      RN = (double) N_old / N_new;
    }
    double RN_err = 0;
    if (N_new > 0 && N_old > 0) {
      double err1 = (double) N_old + N_new;
      double err2 = (double) N_old * N_new;
      double err3 = err1/err2;
      RN_err = RN * sqrt(err3);
    }
    cout << "N_old = " << N_old << endl;
    cout << "N_new = " << N_new << endl;
    cout << "RN = " << RN << endl;
    cout << "err1 = " << err1 << endl;
    cout << "err2 = " << err2 << endl;
    cout << "err3 = " << err3 << endl;
    cout << "RN_err = " << RN_err << endl;
  }

  //#######################################################################
  // Make plots

  gROOT->ProcessLine(".x plotMk_3D_ProjX_D0KPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjY_D0KPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjZ_D0KPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjX_bin_D0KPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjY_bin_D0KPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjZ_bin_D0KPi.C");

  gROOT->ProcessLine(".x plotMk_3D_ProjX_DchKPiPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjY_DchKPiPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjZ_DchKPiPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjX_bin_DchKPiPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjY_bin_DchKPiPi.C");
  gROOT->ProcessLine(".x plotMk_3D_ProjZ_bin_DchKPiPi.C");

}