#include <iostream>
#include <iomanip>
using namespace std;

int fw_binned(int *a1, int *a2, int *ad, double &f1, double &ef1, double &f2, double &ef2){
  int debug_lv=0 ;
  double range=0.5; // +/- 0.5 is physics range
  double tf1=0.5; //init value i.e. middle of physics rage
  double tf2=0.5;
  double minf1=0.5;
  double minf2=0.5;
  int tott=1000000;
  int max_test=60;
  //  int max_test=1;  
  int totd=4000;
  int tlv=0, stagelv=0, paran=0;
  double minLogL=9999999;

  for(int t=0; t<max_test; t++){
    tlv++;
    paran=(tlv/2)%2;
    stagelv=(tlv-2)/4;
    if(tlv!=1 and paran==1){
      if(tlv%2==0)tf1+=range/(pow(2,stagelv)*2);
      else if(tlv%2==1)tf1-=range/(pow(2,stagelv));
    }else if(tlv!=1 and paran==0){
      if(tlv%2==0)tf2+=range/(pow(2,stagelv)*2);
      else if(tlv%2==1)tf2-=range/(pow(2,stagelv));
    }
    double logL=0;
    //      cout<<i<<"  "<<a1[i]<<"  "<<a2[i]<<"   "<<ad[i]<<endl;


    for(int i=0; i<100; i++){ //bin loop for log likelihood

      double ob=ad[i]; //HH observed
      double ex=totd * (tf1*a1[i]/tott + tf2*a2[i]/tott); //HH expected 
      //      cout<<"ex="<<ex<<"   "<<a1[i]<<"   "<<a2[i]<<endl;
      double logP=0;
      if(ex>0) logP= -2*( -ex + ob*log(ex) - log(TMath::Factorial((int)ob)));
      //       cout<<"my logP = "<<logP<<" = "<<-ex<<" + "<<ob*log(ex)<<" - "<<log(TMath::Factorial((int)ob))<<endl;

      logL+=logP;
    }
    if(debug_lv>4) cout<<setw(4)<<tlv<<"th test, stage="<<setw(3)<<stagelv<<", para="<<paran<<", f1="<<setw(8)<<tf1<<", f2="<<setw(8)<<tf2<<",  logL="<<setw(8)<<logL;
    if(logL<minLogL){ //mini
      minLogL=logL;
      minf1=tf1;
      minf2=tf2;
    }
    if(paran==1 and tlv%2==1){
      tf1=minf1;
    }else if(paran==0 and tlv%2==1){
      tf2=minf2;
    }
    if(debug_lv>4) cout<<", minf1="<<setw(8)<<minf1<<", minf2="<<setw(8)<<minf2<<", minlogL="<<minLogL<<"   reso = "<<range/(pow(2,stagelv)*2)<<endl;
  }
  f1=minf1; f2=minf2; ef1=0; ef2=0;
  return 0;
}