/////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////
////
//// 12-1-2008 Nate Baillie edits class so that it returns x from pion
////           contamination and pair contamination
////
////  ratio = exp( A(th) + B(th)*p ) "exp fit to mom"
////
////    A(th) = a + b*th  "line fit to angle"
////    B(th) = c + d*th  "line fit to angle"
////
////  ratio =  exp( a + b*th + c*p + d*th*p )
////
/////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////
////
//// USAGE OF ERROR INDEX indexFlag[] array is the errindex[]:
////
//// indexFlag[0] is unused here
////
//// indexFlag[1] == 0 means no pion background correction
//// indexflag[1] != 0 means apply pion background correction: asymmetry will be changed
////
//// indexFlag[2] == 0 means apply standard pair symmetric background correction
//// indexFlag[2] != 0 means apply standard pair symmetric backround correction with errors added
////
/////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////
////
////  USAGE: After including to your libraries (#include "Conta.h"),
////                 Declare your clas as:
////                 Conta  ContaFix;
////                 Then call with the following parameters:
////
////  ContaFix.bonus_Conta( pi_e_rat,            //// pion to electron ratio as double
////                        pi_e_ratErr,         //// pion to electron ratio error as double
////                        ep_em_rat,           //// positron to electron ratio as double
////                        ep_em_ratErr,        //// positron to electron ratio error as double
////                        iFlag,               //// index for systematic error flag as int *pointer
////                        bin_theta,           //// average polar angle of the bin in Degrees as double
////                        bin_mom,             //// average momentum of the bin in GeV as double
////                        target,              //// target string "ND3" or "NH3"
////                        beam_energy,         //// beam energy in MeV as double
////                        torus,               //// torus field as int (only the sign is used)
////                        debug(optional)      //// 1,2,3 int (optional, it prints things to your screen, don't put it for routine use)
////                      );
////
/////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////


#ifndef __Conta_h__
#define __Conta_h__

#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <math.h>
#include <fstream>
#include <string>
#include <cmath>

using namespace std;
#include "utils.h"


//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
///////////////// Define global variables and declare clas Conta /////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

//// global  parameters for standard pion contamination

double    a_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double aErr_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    b_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double bErr_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    c_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double cErr_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    d_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double dErr_mat_stand[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};

//// global parameters for total pion contamination

double    a_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double aErr_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    b_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double bErr_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    c_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double cErr_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    d_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double dErr_mat_total[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};

//// global parameters for standard ePeM contamination

double    a_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double aErr_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    b_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double bErr_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    c_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double cErr_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double    d_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};
double dErr_mat_epem[2][9][2] = {{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},},{{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},{0.,0.},}};

//// global values for positron asymmetries

double asym_posiMat[2][9][2][11][15] = {0.,0.};
double asymErr_posiMat[2][9][2][11][15] = {0.,0.};

const Int_t NEG1B_E = 8;
double eg1b_ebeam[NEG1B_E]; 

typedef struct
{
    int C1;
    string Tg;
    double E;
    int Torus;
    int iTg;
    int iE;
    int iTorus;
} indextable_T;

class Conta
{
    public:

    Conta()
    {
        Read_Index();
        Read_StdPimPar();
        Read_TotPimPar();
        Read_StdEpEmPar();
        Read_StdPosiAsym();
	eg1b_ebeam[0] = 1606.0;
	eg1b_ebeam[1] = 1723.0;
	eg1b_ebeam[2] = 2286.0;
	eg1b_ebeam[3] = 2561.0;
	eg1b_ebeam[4] = 4238.0;
	eg1b_ebeam[5] = 5615.0;
	eg1b_ebeam[6] = 5725.0;
	eg1b_ebeam[7] = 5743.0;
    };
    ~Conta(){ };

    void Find_Index( string target, double energy, int torus,
                    int &iconf_c1, int &iconf_T, int &iconf_E, int &iconf_I,
                    int &oconf_c1, int &oconf_T, int &oconf_E, int &oconf_I,
                    int db = 0 );

    void bonus_Conta( double &pi_e_rat, double &pi_e_ratErr,
                    double &ep_em_rat, double &ep_em_ratErr,
                    int *indexFlag, double theta, double mom,
                    string target, double energy, int torus,
                    int db = 0 );

    void Eg1b_Conta( double &pi_e_rat, double &pi_e_ratErr,
                    double &ep_em_rat, double &ep_em_ratErr,
                    int *indexFlag, double theta, double mom,
                    string target, double energy, int torus,
                    int db = 0 );

    void Pair_Conta( double &ep_em_rat, double &ep_em_ratErr,
                    int *indexFlag, double theta, double mom,
                    int iconf_c1, int iconf_T, int iconf_E, int iconf_I,
                    int oconf_c1, int oconf_T, int oconf_E, int oconf_I,
                    int db = 0 );

    void Pion_Conta( double &pi_e_rat, double &pi_e_ratErr,
                    int *indexFlag, double theta, double mom,
                    int conf_c1, int conf_T, int conf_E, int conf_I,
                    int db = 0 );

    private:
    void Read_Index();
    void GetOne_Index(int index,int &c1,string &target,double &energy,int &torus,int &iT,int &iE,int &iI);
    void Read_StdPimPar();
    void Read_TotPimPar();
    void Read_StdEpEmPar();
    void Read_StdPosiAsym();

    private:
    indextable_T SrIndexTable[40];


};


//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
///////////// Functions to Read the Index File and to Find the index values //////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

void Conta::Read_Index()
{
    char sline[1010];
    ifstream indexInfo;
    //indexInfo.open("/home/claseg1/eg1b/analysis/asymFix/paramFiles/indexTable.txt",ios::in);
    indexInfo.open("./asymFix/paramFiles/indexTable.txt",ios::in);
    assert (!indexInfo.fail());
    indexInfo.getline( sline, 1000 );
    indexInfo.getline( sline, 1000 );
    indexInfo.getline( sline, 1000 );
    indexInfo.getline( sline, 1000 );

    int iLine=0;
    double xinput;
    while (indexInfo >> xinput)
    {
        SrIndexTable[iLine].C1 = int(xinput);
        indexInfo >> SrIndexTable[iLine].Tg;
        indexInfo >> SrIndexTable[iLine].E;
        indexInfo >> SrIndexTable[iLine].Torus;
        indexInfo >> SrIndexTable[iLine].iTg;
        indexInfo >> SrIndexTable[iLine].iE;
        indexInfo >> SrIndexTable[iLine].iTorus;
        iLine++;
    }
    indexInfo.close();
    //cout << "Read the index table " << endl;
    return ;

}


void Conta::GetOne_Index(int index,int &c1,string &target,double &energy,int &torus,int &iT,int &iE,int &iI)
{
    c1 = SrIndexTable[index].C1;
    target = SrIndexTable[index].Tg;
    energy = SrIndexTable[index].E;
    torus = SrIndexTable[index].Torus;
    iT = SrIndexTable[index].iTg;
    iE = SrIndexTable[index].iE;
    iI = SrIndexTable[index].iTorus;
}



void Conta::Find_Index( string target, double energy, int torus,
                       int &iconf_c1, int &iconf_T, int &iconf_E, int &iconf_I,
                       int &oconf_c1, int &oconf_T, int &oconf_E, int &oconf_I,
                       int db )
{


    int conf_c1;

    if( fabs( energy - 1606.0 ) < 10.0 ) conf_c1 = 1;
    if( fabs( energy - 1723.0 ) < 10.0 ) conf_c1 = 2;
    if( fabs( energy - 2286.0 ) < 10.0 ) conf_c1 = 3;
    if( fabs( energy - 2561.0 ) < 10.0 ) conf_c1 = 4;
    if( fabs( energy - 4238.0 ) < 10.0 ) conf_c1 = 5;
    if( fabs( energy - 5615.0 ) < 10.0 ) conf_c1 = 6;
    if( fabs( energy - 5725.0 ) < 10.0 ) conf_c1 = 7;
    if( fabs( energy - 5743.0 ) < 10.0 ) conf_c1 = 8;


    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    /////////////////////// read the index file //////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////


    int index_infoC1;
    string index_infoT;
    double index_infoE;
    int index_infoI;
    int index_T;
    int index_E;
    int index_I;

    for (int i=0;i<40;i++)
    {
        //// get one record from index table
        GetOne_Index(i,index_infoC1,index_infoT,index_infoE,index_infoI,index_T,index_E,index_I);

        if(
            target == index_infoT
        && fabs(energy - index_infoE) < 10.0
        && SignOf(torus) == SignOf(index_infoI)
        && conf_c1 == index_E
        && conf_c1 == index_infoC1
        )
        {

            //////////////////////////////////////////////////////////////////////////////////////////////////
            ////////////////////////////////// begin debug ///////////////////////////////////////////////////
            //////////////////////////////////////////////////////////////////////////////////////////////////
            if(db == 1 || db == 3 )
            {
                cout << "*************** NEW BIN *****************************" << endl;
                cout << endl;
                cout << "index of the configuration is determined" << endl;
                cout << "C1     : " << conf_c1 << "  =  " << index_infoC1 << endl;
                cout << "Target : " << target  << "  =  " << index_infoT << endl;
                cout << "Energy : " << energy  << "  =  " << index_infoE << endl;
                cout << "Torus  : " << torus   << "  =  " << index_infoI << endl;
            }
            //////////////////////////////////////////////////////////////////////////////////////////////////
            ////////////////////////////////// end debug /////////////////////////////////////////////////////
            //////////////////////////////////////////////////////////////////////////////////////////////////

            iconf_c1 = conf_c1;


            //// choose data, exit if data does not exist
            //// arrange data configurations according to Target and Torus
            //// find the corresponding config for opposite torus polarity
            //// oconf_c1 will be used for systematic error calculations

            if(target == "ND3" && SignOf(torus) == 1)
            {
                if      ( iconf_c1 == 1 ) { oconf_c1=1; }
                else if ( iconf_c1 == 4 ) { oconf_c1=4; }
                else if ( iconf_c1 == 5 ) { oconf_c1=5; }
                else if ( iconf_c1 == 6 ) { oconf_c1=7; }
                else if ( iconf_c1 == 7 ) { oconf_c1=7; }
                else
                {
                    cout <<"Data requested does not exist!" << endl;
                    return;
                }

            }

            if(target == "ND3" && SignOf(torus) == -1)
            {
                if      ( iconf_c1 == 1 ) { oconf_c1=1; }
                else if ( iconf_c1 == 2 ) { oconf_c1=1; }
                else if ( iconf_c1 == 4 ) { oconf_c1=4; }
                else if ( iconf_c1 == 5 ) { oconf_c1=5; }
                else if ( iconf_c1 == 7 ) { oconf_c1=7; }
                else if ( iconf_c1 == 8 ) { oconf_c1=7; }
                else
                {
                    cout <<"Data requested does not exist!" << endl;
                    return;
                }

            }

            if(target == "NH3" && SignOf(torus) == 1)
            {
                if      ( iconf_c1 == 1 ) { oconf_c1=1; }
                else if ( iconf_c1 == 3 ) { oconf_c1=4; }
                else if ( iconf_c1 == 5 ) { oconf_c1=5; }
                else if ( iconf_c1 == 6 ) { oconf_c1=7; }
                else if ( iconf_c1 == 7 ) { oconf_c1=7; }
                else
                {
                    cout <<"Data requested does not exist!" << endl;
                    return;
                }

            }

            if(target == "NH3" && SignOf(torus) == -1)
            {
                if      ( iconf_c1 == 1 ) { oconf_c1=1; }
                else if ( iconf_c1 == 2 ) { oconf_c1=1; }
                else if ( iconf_c1 == 4 ) { oconf_c1=3; }
                else if ( iconf_c1 == 5 ) { oconf_c1=5; }
                else if ( iconf_c1 == 7 ) { oconf_c1=7; }
                else if ( iconf_c1 == 8 ) { oconf_c1=7; }
                else
                {
                    cout <<"Data requested does not exist!" << endl;
                    return;
                }

            }

            iconf_T = index_T;
            iconf_E = index_E;
            iconf_I = index_I;

            oconf_T = index_T;
            oconf_E = oconf_c1;
            oconf_I = abs(index_I-1);

            //////////////////////////////////////////////////////////////////////////////////////////////////
            ////////////////////////////////// begin debug ///////////////////////////////////////////////////
            //////////////////////////////////////////////////////////////////////////////////////////////////
            if( db == 1 || db == 2 || db == 3 )
            {
                cout << "*************** NEW BIN DETAILS **********************" << endl;
                cout << endl;

                cout << "The index values are: "<< endl;
                cout << "Config" << "\t";
                cout << "index_C" << "\t";
                cout << "index_T" << "\t";
                cout << "index_E" << "\t";
                cout << "index_I" << endl;

                cout << "iConf" << "\t";
                cout << iconf_c1 << "\t";
                cout << iconf_T  << "\t";
                cout << iconf_E  << "\t";
                cout << iconf_I  << endl;

                cout << "oConf" << "\t";
                cout << oconf_c1 << "\t";
                cout << oconf_T  << "\t";
                cout << oconf_E  << "\t";
                cout << oconf_I  << endl;
            }
            //////////////////////////////////////////////////////////////////////////////////////////////////
            ////////////////////////////////// end debug /////////////////////////////////////////////////////
            //////////////////////////////////////////////////////////////////////////////////////////////////


            return;

        }

        else if(
            target == index_infoT
        && fabs(energy - index_infoE) < 10.0
        && SignOf(torus) == SignOf(index_infoI)
        && conf_c1 == index_E
        && conf_c1 != index_infoC1
        )
        {

            oconf_c1 = conf_c1;
            iconf_c1 = index_infoC1; //// use another energy for this data

            iconf_T = index_T;
            iconf_E = iconf_c1;
            iconf_I = index_I;

            oconf_T = index_T;
            oconf_E = oconf_c1;
            oconf_I = abs(index_I-1);

            //////////////////////////////////////////////////////////////////////////////////////////////////
            ////////////////////////////////// begin debug ///////////////////////////////////////////////////
            //////////////////////////////////////////////////////////////////////////////////////////////////
            if( db == 1 || db == 2 || db == 3 )
            {
                cout << "*************** NEW BIN *****************************" << endl;
                cout << endl;

                cout << "There is no data for this configuration" << endl;
                cout << "Data from another configuration will be used." << endl;
                cout << "The index values are: "<< endl;
                cout << "Config" << "\t";
                cout << "index_C" << "\t";
                cout << "index_T" << "\t";
                cout << "index_E" << "\t";
                cout << "index_I" << endl;

                cout << "iConf" << "\t";
                cout << iconf_c1 << "\t";
                cout << iconf_T  << "\t";
                cout << iconf_E  << "\t";
                cout << iconf_I  << endl;

                cout << "oConf" << "\t";
                cout << oconf_c1 << "\t";
                cout << oconf_T  << "\t";
                cout << oconf_E  << "\t";
                cout << oconf_I  << endl;
            }
            //////////////////////////////////////////////////////////////////////////////////////////////////
            ////////////////////////////////// end debug /////////////////////////////////////////////////////
            //////////////////////////////////////////////////////////////////////////////////////////////////

            return;

        }


    }


}

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
/////////// Functions to read the pion parameters and apply pion correction //////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

void Conta::Read_StdPosiAsym()
{


    //// initialize to zero

    for(int iconf_T = 0; iconf_T <= 1; iconf_T++)
    {
        for(int iconf_E = 0; iconf_E <= 8; iconf_E++)
        {
            for(int iconf_I = 0; iconf_I <= 1; iconf_I++)
            {
                for(int th = 0; th < 11; th++)
                {
                    for(int p = 0; p < 15; p++)
                    {

                        asym_posiMat[iconf_T][iconf_E][iconf_I][th][p] = 0.;
                        asymErr_posiMat[iconf_T][iconf_E][iconf_I][th][p] = 0.;

                    }
                }
            }
        }
    }



    //// read all values of positron asymmetries for all data into global arrays

    for(int iconf_T = 0; iconf_T <= 1; iconf_T++)
    {
        for(int iconf_E = 1; iconf_E <= 8; iconf_E++)
        {
            for(int iconf_I = 0; iconf_I <= 1; iconf_I++)
            {
                int fread;
                char sline[1010];


                //// read positron asymmetry
                //// temporary solution until we get parametrized asymmetries

                int th,p;
                char n_torus[1];
                char n_energy[5];
                char n_target[5];

                if     (iconf_T == 1) { sprintf(n_target,"ND3"); }
                else if(iconf_T == 0) { sprintf(n_target,"NH3"); }
                else { cout << "target config error!" << endl; return; }

                if     (iconf_E == 1) { sprintf(n_energy,"1.6"); }
                else if(iconf_E == 2) { sprintf(n_energy,"1.7"); }
                else if(iconf_E == 3) { sprintf(n_energy,"2.3"); }
                else if(iconf_E == 4) { sprintf(n_energy,"2.5"); }
                else if(iconf_E == 5) { sprintf(n_energy,"4.2"); }
                else if(iconf_E == 6) { sprintf(n_energy,"5.6"); }
                else if(iconf_E == 7) { sprintf(n_energy,"5.7"); }
                else if(iconf_E == 8) { sprintf(n_energy,"5.8"); }
                else { cout << "beamen config error!" << endl; return; }

                if     (iconf_I == 1) { sprintf(n_torus,"+"); }
                else if(iconf_I == 0) { sprintf(n_torus,"-"); }
                else { cout << "torus config error! " << endl; return; }


                char asym_posiFileName[1010];
                //sprintf(asym_posiFileName,"/home/claseg1/eg1b/analysis/asymFix/paramFiles/posiFiles/posiAsym_%s_%s%s.dat",n_target,n_energy,n_torus);
                sprintf(asym_posiFileName,"./asymFix/paramFiles/posiFiles/posiAsym_%s_%s%s.dat",n_target,n_energy,n_torus);
                //cout<<"Read positron files: "<<asym_posiFileName<<endl;
                ifstream asym_posiFile;
                asym_posiFile.open(asym_posiFileName,ios::in);
                assert ( !asym_posiFile.fail() );
                asym_posiFile.getline(sline,1000);

                while(asym_posiFile >> fread)
                {

                    th = fread;
                    asym_posiFile >> p;
                    asym_posiFile >> asym_posiMat[iconf_T][iconf_E][iconf_I][th][p]
                                  >> asymErr_posiMat[iconf_T][iconf_E][iconf_I][th][p];

                    /*
                    //// see if you read correctly
                    cout<<iconf_T<<" "<<iconf_E <<" "<<iconf_I<<" "<<th<<" "<<p<<" "
                        <<asym_posiMat[iconf_T][iconf_E][iconf_I][th][p]<< "  "
                        <<asymErr_posiMat[iconf_T][iconf_E][iconf_I][th][p]<< "  "
                        <<endl;
                    //*/

                }

                asym_posiFile.close();

            }

        }

    }

    return;

}


void Conta::Read_StdEpEmPar()
{
    int fread;
    char sline[1010];
    char ePeM_Ratio_abcdFileName[1010];
    //sprintf(ePeM_Ratio_abcdFileName,"/home/claseg1/eg1b/analysis/asymFix/paramFiles/abcd_ePeM_Ratio.txt");
    sprintf(ePeM_Ratio_abcdFileName,"./asymFix/paramFiles/abcd_ePeM_Ratio.txt");
    ifstream ePeM_Ratio_abcdFile;
    ePeM_Ratio_abcdFile.open(ePeM_Ratio_abcdFileName,ios::in);
    assert ( !ePeM_Ratio_abcdFile.fail() );
    //ePeM_Ratio_abcdFile.getline(s,1000);


    int conta_ic1;
    string conta_infoT;
    int conta_infoE;
    int conta_infoI;
    int conta_T;
    int conta_E;
    int conta_I;

    while(ePeM_Ratio_abcdFile >> fread)
    {
        conta_ic1 = fread;
        ePeM_Ratio_abcdFile
        >> conta_infoT
        >> conta_infoE
        >> conta_infoI
        >> conta_T
        >> conta_E
        >> conta_I;

        ePeM_Ratio_abcdFile
        >>    a_mat_epem[conta_T][conta_E][conta_I]
        >> aErr_mat_epem[conta_T][conta_E][conta_I]
        >>    b_mat_epem[conta_T][conta_E][conta_I]
        >> bErr_mat_epem[conta_T][conta_E][conta_I]
        >>    c_mat_epem[conta_T][conta_E][conta_I]
        >> cErr_mat_epem[conta_T][conta_E][conta_I]
        >>    d_mat_epem[conta_T][conta_E][conta_I]
        >> dErr_mat_epem[conta_T][conta_E][conta_I];

    }
    ePeM_Ratio_abcdFile.close();
    //cout << "Read the ePeM files " << endl;
    return;

}


void Conta::Pair_Conta(double &ep_em_rat, double &ep_em_ratErr,
                       int *indexFlag, double theta, double mom,
                       int iconf_c1, int iconf_T, int iconf_E, int iconf_I,
                       int oconf_c1, int oconf_T, int oconf_E, int oconf_I,
                       int db )


{
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    /////////////////////// find bin numbers /////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    makeBins();

    for ( int bin=0; bin<15; bin++ )
    {
        if ( (mom >= pBin_min[bin]) && (mom < pBin_max[bin]) && (mom != 0.0) )
        {
            pBin=bin;
            break;
        } else { pBin = 0; }
    }

    for ( int bin=0; bin<11; bin++ )
    {
        if ( (theta >= thetaBin_min[bin]) && (theta < thetaBin_max[bin]) && (theta != 0.0) )
        {
            thetaBin=bin;
            break;
        } else { thetaBin = 0; }
    }



    //////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////// begin debug ///////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////
    if( db == 1 || db == 2 )
    {
        cout << endl;
        cout << "******* PAIR CONTAMINATION CORRECTION FOR THE FOLLOWING KINEMATIC : "<< endl;
        cout <<"theta = " << theta << endl;
        cout <<"mom   = " << mom << endl;
        cout << endl;
    }
    //////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////// end debug /////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////


    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////// make corrections //////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    //// the element of the matrices that we read above that corresponds to
    //// our current configuration is element [iconf_T][iconf_E][iconf_I]

    double a1;
    double a1Err;
    double b1;
    double b1Err ;
    double c1;
    double c1Err;
    double d1;
    double d1Err;
    double X1;
    double X1_min;
    double X1_max;
    double X1_err;

    double a2;
    double a2Err;
    double b2;
    double b2Err ;
    double c2;
    double c2Err;
    double d2;
    double d2Err;
    double X2;
    double X2_min;
    double X2_max;
    double X2_err;

    double x;
    double xErr;


    a1       = a_mat_epem[iconf_T][iconf_E][iconf_I];
    a1Err = aErr_mat_epem[iconf_T][iconf_E][iconf_I];
    b1       = b_mat_epem[iconf_T][iconf_E][iconf_I];
    b1Err = bErr_mat_epem[iconf_T][iconf_E][iconf_I];
    c1       = c_mat_epem[iconf_T][iconf_E][iconf_I];
    c1Err = cErr_mat_epem[iconf_T][iconf_E][iconf_I];
    d1       = d_mat_epem[iconf_T][iconf_E][iconf_I];
    d1Err = dErr_mat_epem[iconf_T][iconf_E][iconf_I];

    a2       = a_mat_epem[oconf_T][oconf_E][oconf_I];
    a2Err = aErr_mat_epem[oconf_T][oconf_E][oconf_I];
    b2       = b_mat_epem[oconf_T][oconf_E][oconf_I];
    b2Err = bErr_mat_epem[oconf_T][oconf_E][oconf_I];
    c2       = c_mat_epem[oconf_T][oconf_E][oconf_I];
    c2Err = cErr_mat_epem[oconf_T][oconf_E][oconf_I];
    d2       = d_mat_epem[oconf_T][oconf_E][oconf_I];
    d2Err = dErr_mat_epem[oconf_T][oconf_E][oconf_I];


    //// find average systematic error for all theta-mom values

    double X1_sum_num = 0.0;
    double X1_sum_den = 0.0;
    double X2_sum_num = 0.0;
    double X2_sum_den = 0.0;
    double t_binave, p_binave;
    for ( int bint=2; bint<=7; bint++ )
    {
        for ( int binp=3; binp<=10; binp++ )
        {

            t_binave = thetaBin_ave[bint];
            p_binave = pBin_ave[binp];

            X1 = exp(a1 + b1*t_binave + c1*p_binave + d1*t_binave*p_binave );
            X1_max = exp((a1+a1Err) + (b1+b1Err)*t_binave + (c1+c1Err)*p_binave + (d1+d1Err)*t_binave*p_binave );
            X1_min = exp((a1-a1Err) + (b1-b1Err)*t_binave + (c1-c1Err)*p_binave + (d1-d1Err)*t_binave*p_binave );
            X1_err = fabs(X1_max-X1_min)/2.0;
            if(X1_err > 0.0)
            {
                X1_sum_num +=  X1/(X1_err*X1_err);
                X1_sum_den +=  1.0/(X1_err*X1_err);
            }
            else
            {
                X1_sum_num += X1;
                X1_sum_den += 1.0;
            }

            X2 = exp(a2 + b2*t_binave + c2*p_binave + d2*t_binave*p_binave );
            X2_max = exp((a2+a2Err) + (b2+b2Err)*t_binave + (c2+c2Err)*p_binave + (d2+d2Err)*t_binave*p_binave );
            X2_min = exp((a2-a2Err) + (b2-b2Err)*t_binave + (c2-c2Err)*p_binave + (d2-d2Err)*t_binave*p_binave );
            X2_err = fabs(X2_max-X2_min)/2.0;
            if(X2_err > 0.0)
            {
                X2_sum_num +=  X2/(X2_err*X2_err);
                X2_sum_den +=  1.0/(X2_err*X2_err);
            }
            else
            {
                X2_sum_num += X2;
                X2_sum_den += 1.0;
            }

        }

    }

    X1 = X1_sum_num/X1_sum_den;
    X1_err = sqrt(1.0/X1_sum_den);

    X2 = X2_sum_num/X2_sum_den;
    X2_err = sqrt(1.0/X2_sum_den);


    //// now get the final values

    x =  exp( a1 + b1*theta + c1*mom + d1*theta*mom ); //// local value for this [theta,p] point
    xErr = fabs(X1 - X2)/2.0;                          //// average syst error for all values

    ep_em_rat = x;
    ep_em_ratErr = xErr;

    if(indexFlag[2] != 0)
    {
        x = x + xErr;
    }


    //////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////// begin debug ///////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////

    if(db == 1)
    {
        cout << "x    = " << x <<endl;
        cout << "xErr = " << xErr <<endl;
        cout << endl;
    }
    if(db == 2)
    {
        cout << "epem abcd parameters used for configuration: "<< endl;
        cout << "a1 = " << a1 << endl;
        cout << "b1 = " << b1 << endl;
        cout << "c1 = " << c1 << endl;
        cout << "d1 = " << d1 << endl;
        cout << "average e+/e- ratio of this configuration : " << X1 << endl;
        cout << endl;

        cout << "epem abcd parameters used for error calculaton: "<< endl;
        cout << "a2 = " << a2 << endl;
        cout << "b2 = " << b2 << endl;
        cout << "c2 = " << c2 << endl;
        cout << "d2 = " << d2 << endl;
        cout << "average  e+/e- ratio of counter configuration : " << X2 << endl;
        cout << endl;

        cout << "final e+/e- ratio = " << x << "    err = " << xErr << endl;
        cout << endl;

        char asym_posiFileName[1010];
        cout << "file used for positron asymmetry is: " << endl;
        {
            int th,p;
            char n_torus[1];
            char n_energy[5];
            char n_target[5];

            if     (iconf_T == 1) { sprintf(n_target,"ND3"); }
            else if(iconf_T == 0) { sprintf(n_target,"NH3"); }
            else { cout << "target config error!" << endl; return; }

            if     (iconf_E == 1) { sprintf(n_energy,"1.6"); }
            else if(iconf_E == 2) { sprintf(n_energy,"1.7"); }
            else if(iconf_E == 3) { sprintf(n_energy,"2.3"); }
            else if(iconf_E == 4) { sprintf(n_energy,"2.5"); }
            else if(iconf_E == 5) { sprintf(n_energy,"4.2"); }
            else if(iconf_E == 6) { sprintf(n_energy,"5.6"); }
            else if(iconf_E == 7) { sprintf(n_energy,"5.7"); }
            else if(iconf_E == 8) { sprintf(n_energy,"5.8"); }
            else { cout << "beamen config error!" << endl; return; }

            if     (iconf_I == 1) { sprintf(n_torus,"+"); }
            else if(iconf_I == 0) { sprintf(n_torus,"-"); }
            else { cout << "torus config error! " << endl; return; }

            //sprintf(asym_posiFileName,"/home/claseg1/eg1b/analysis/asymFix/paramFiles/posiFiles/posiAsym_%s_%s%s.dat",n_target,n_energy,n_torus);
            sprintf(asym_posiFileName,"./asymFix/paramFiles/posiFiles/posiAsym_%s_%s%s.dat",n_target,n_energy,n_torus);
        }
        cout << asym_posiFileName << endl;
    }

    //////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////// end debug /////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////



    return;


}


//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////// Functions to read the ePeM files and to apply epem correction ///////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

void Conta::Read_StdPimPar()
{
    int fread;
    char sline[1010];
    char stand_bumpConta_abcdFileName[1010];
    //sprintf(stand_bumpConta_abcdFileName,"/home/claseg1/eg1b/analysis/asymFix/paramFiles/abcd_stand_bumpConta.txt");
    sprintf(stand_bumpConta_abcdFileName,"./asymFix/paramFiles/abcd_stand_bumpConta.txt");
    ifstream stand_bumpConta_abcdFile;
    stand_bumpConta_abcdFile.open(stand_bumpConta_abcdFileName,ios::in);
    assert ( !stand_bumpConta_abcdFile.fail() );
    //stand_bumpConta_abcdFile.getline(sline,1000);

    int conta_ic1;
    string conta_infoT;
    int conta_infoE;
    int conta_infoI;
    int conta_T;
    int conta_E;
    int conta_I;

    while(stand_bumpConta_abcdFile >> fread)
    {
        conta_ic1 = fread;
        stand_bumpConta_abcdFile
        >> conta_infoT
        >> conta_infoE
        >> conta_infoI
        >> conta_T
        >> conta_E
        >> conta_I;

        stand_bumpConta_abcdFile
        >>    a_mat_stand[conta_T][conta_E][conta_I]
        >> aErr_mat_stand[conta_T][conta_E][conta_I]
        >>    b_mat_stand[conta_T][conta_E][conta_I]
        >> bErr_mat_stand[conta_T][conta_E][conta_I]
        >>    c_mat_stand[conta_T][conta_E][conta_I]
        >> cErr_mat_stand[conta_T][conta_E][conta_I]
        >>    d_mat_stand[conta_T][conta_E][conta_I]
        >> dErr_mat_stand[conta_T][conta_E][conta_I];

    }
    stand_bumpConta_abcdFile.close();
    //cout << "Read the StdPim tables " << endl;
    return ;

}



void Conta::Read_TotPimPar()
{
    int fread;
    char sline[1010];
    char total_pimConta_abcdFileName[1010];
    //sprintf(total_pimConta_abcdFileName,"/home/claseg1/eg1b/analysis/asymFix/paramFiles/abcd_total_pimConta.txt");
    sprintf(total_pimConta_abcdFileName,"./asymFix/paramFiles/abcd_total_pimConta.txt");
    ifstream total_pimConta_abcdFile;
    total_pimConta_abcdFile.open(total_pimConta_abcdFileName,ios::in);
    assert ( !total_pimConta_abcdFile.fail() );
    //total_pimConta_abcdFile.getline(sline,1000);

    int conta_ic1;
    string conta_infoT;
    int conta_infoE;
    int conta_infoI;
    int conta_T;
    int conta_E;
    int conta_I;

    while(total_pimConta_abcdFile >> fread)
    {
        conta_ic1 = fread;
        total_pimConta_abcdFile
        >> conta_infoT
        >> conta_infoE
        >> conta_infoI
        >> conta_T
        >> conta_E
        >> conta_I;

        total_pimConta_abcdFile
        >>    a_mat_total[conta_T][conta_E][conta_I]
        >> aErr_mat_total[conta_T][conta_E][conta_I]
        >>    b_mat_total[conta_T][conta_E][conta_I]
        >> bErr_mat_total[conta_T][conta_E][conta_I]
        >>    c_mat_total[conta_T][conta_E][conta_I]
        >> cErr_mat_total[conta_T][conta_E][conta_I]
        >>    d_mat_total[conta_T][conta_E][conta_I]
        >> dErr_mat_total[conta_T][conta_E][conta_I];

    }
    total_pimConta_abcdFile.close();
    //cout << "Read the TotPim tables " << endl;
    return ;

}



void Conta::Pion_Conta(double &pi_e_rat, double &pi_e_ratErr,
                       int *indexFlag, double theta, double mom,
                       int conf_c1, int conf_T, int conf_E, int conf_I,
                       int db )
{

    //// determine the momentum range for appropriate correction
    int range1 = 0;
    int range2 = 0;
    char mrange[100];

    if(mom < 2.7)
    {
        range1 = 1; range2 = 0;
    }
    if(mom >= 2.7 && mom <= 3.0)
    {
        range1 = 1; range2 = 1;
    }
    if(mom > 3.0)
    {
        range1 = 0; range2 = 1;
    }


    /*

    //// Put some other tool here to catch situations where no parameters can be found!!

    if(found1 == 0 && range1 == 1)
    {
        cout << "didn't find pion 1 parameters" << endl;
        return;
    }
    if(found2 == 0 && range2 == 1)
    {
        cout << "didn't find pion 2 parameters" << endl;
        return;
    }

    //*/

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////// make corrections //////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    //// the element of the matrices that we read above that corresponds to
    //// our current configuration is element [conf_T][conf_E][conf_I]

    double a;
    double aErr;
    double b;
    double bErr ;
    double c;
    double cErr;
    double d;
    double dErr;

    double x;

    if(range1 == 1 && range2 == 0)
    {
        a = a_mat_stand[conf_T][conf_E][conf_I];
        b = b_mat_stand[conf_T][conf_E][conf_I];
        c = c_mat_stand[conf_T][conf_E][conf_I];
        d = d_mat_stand[conf_T][conf_E][conf_I];
        sprintf(mrange,"low momentum range p < 2.7 GeV");
    }

    else if(range1 == 0 && range2 == 1)
    {
        a = a_mat_total[conf_T][conf_E][conf_I];
        b = b_mat_total[conf_T][conf_E][conf_I];
        c = c_mat_total[conf_T][conf_E][conf_I];
        d = d_mat_total[conf_T][conf_E][conf_I];
        sprintf(mrange,"high momentum range p > 3.0 GeV");
    }

    else if(range1 == 1 && range2 == 1)
    {
        a = (a_mat_stand[conf_T][conf_E][conf_I]
             +a_mat_total[conf_T][conf_E][conf_I])/2.0;
        b = (b_mat_stand[conf_T][conf_E][conf_I]
             +b_mat_total[conf_T][conf_E][conf_I])/2.0;
        c = (c_mat_stand[conf_T][conf_E][conf_I]
             +c_mat_total[conf_T][conf_E][conf_I])/2.0;
        d = (d_mat_stand[conf_T][conf_E][conf_I]
             +d_mat_total[conf_T][conf_E][conf_I])/2.0;
        sprintf(mrange,"medium momentum range p > 2.7 && p < 3.0 GeV");
    }
    else
    {
        return;
    }

    x =  exp(a + b*theta + c*mom + d*theta*mom );
pi_e_rat = x;
pi_e_ratErr = x;
    
    if( x >= 1 )
    {
        cout << "dude, invalid pion contamination value!!"<< endl;
        cout << "******* ERROR IN PION CONTAMINATION CORRECTIONS : "<< endl;
        cout << mrange << endl;
        cout << "abcd used for pion contamination: "<< endl;
        cout << "a = " << a << endl;
        cout << "b = " << b << endl;
        cout << "c = " << c << endl;
        cout << "d = " << d << endl;
	cout << "Theta (deg) = " << theta << ", Momentum(GeV) = " << mom << endl;
    }
    //////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////// begin debug ///////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////
    if(db == 1)
    {
        cout << "******* PION CONTAMINATION CORRECTIONS : "<< endl;
        cout << mrange << endl;
        cout << "Pion to Electron ratio    : " << x << endl;
	cout << "Theta (deg) = " << theta << ", Momentum(GeV) = " << mom << endl;
        cout << endl;
    }
    if(db == 2)
    {
        cout << endl;
        cout << endl;
        cout << endl;
        cout << "******* PION CONTAMINATION CORRECTIONS : "<< endl;
        cout << mrange << endl;
        cout << "abcd used for pion contamination: "<< endl;
        cout << "a = " << a << endl;
        cout << "b = " << b << endl;
        cout << "c = " << c << endl;
        cout << "d = " << d << endl;
	cout << "Theta (deg) = " << theta << ", Momentum(GeV) = " << mom << endl;
        cout << endl;
    }
    //////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////// end debug /////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////



    return;


}

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
////////////////////// Function that Applies all corrections /////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

void Conta::Eg1b_Conta( double &pi_e_rat, double &pi_e_ratErr,
                        double &ep_em_rat, double &ep_em_ratErr,
                        int *indexFlag, double theta, double mom,
                        string target, double energy, int torus,
                        int db )

{

    //// energy should be in MeV
    //// target string "ND3" or "NH3"
    //// theta in degrees
    //// momentum in GeV

    int iconf_c1, iconf_T, iconf_E, iconf_I;
    int oconf_c1, oconf_T, oconf_E, oconf_I;

    Find_Index( target, energy, torus,
               iconf_c1, iconf_T, iconf_E, iconf_I,
               oconf_c1, oconf_T, oconf_E, oconf_I,db);

    Pair_Conta( ep_em_rat, ep_em_ratErr,
               indexFlag, theta, mom,
               iconf_c1, iconf_T, iconf_E, iconf_I,
               oconf_c1, oconf_T, oconf_E, oconf_I,db);

    Pion_Conta( pi_e_rat,pi_e_ratErr,
               indexFlag, theta, mom,
               iconf_c1, iconf_T, iconf_E, iconf_I,db);

    return;


}

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
/////////// Function that determines which eg1b energies the bonus energy is between 
////////////// then returns a weighted average of those ratios for bonus///////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

void Conta::bonus_Conta( double &pi_e_rat, double &pi_e_ratErr,
                        double &ep_em_rat, double &ep_em_ratErr,
                        int *indexFlag, double theta, double mom,
                        string target, double energy, int torus,
                        int db )

{

  double pie1, epie1, pie2, epie2; 
  double epem1, eepem1,epem2, eepem2; 
  double weight1, weight2, e1, e2;

  e1 = e2 = weight1 = weight2 = -1000.0;

  if(energy < eg1b_ebeam[0]) e1 = e2 = eg1b_ebeam[0];
  if(energy > eg1b_ebeam[NEG1B_E-1]) e1 = e2 = eg1b_ebeam[NEG1B_E-1];
  
  for(int ii = 0; ii < NEG1B_E; ii++)
    {
      if(energy == eg1b_ebeam[ii]) 
	{e1 = e2 = eg1b_ebeam[ii]; break;}

      if(energy > eg1b_ebeam[ii] && energy < eg1b_ebeam[ii+1])
	{
	  e1 = eg1b_ebeam[ii];
	  e2 = eg1b_ebeam[ii+1];
	  weight1 = 1/(energy - eg1b_ebeam[ii]);
	  weight2 = 1/(eg1b_ebeam[ii+1] - energy);
	  break;
	}
      
    }

  //check to see if something went wrong in determining the energies and weights
  if(e1 == -1000.0 || e2 == -1000.0 || weight1 == -1000.0 || weight2 == -1000.0)
    {
      cout << "ERROR!!!!! Cannot find eg1b beam energies to use to determine contaminations" << endl;
      return;
    }

  if(e1 == e2)
    {

      Eg1b_Conta( pi_e_rat, pi_e_ratErr,
		  ep_em_rat, ep_em_ratErr,
		  indexFlag, theta, mom,
		  target, e1, torus, db );
    }

  if(e1 != e2)
    {
      Eg1b_Conta( pie1, epie1,
		  epem1, eepem1,
		  indexFlag, theta, mom,
		  target, e1, torus, db );

      Eg1b_Conta( pie2, epie2,
		  epem2, eepem2,
		  indexFlag, theta, mom,
		  target, e2, torus, db );
      
      pi_e_rat    = ( pie1 * weight1 + pie2 * weight2 )   / (weight1 + weight2);
      pi_e_ratErr = ( epie1 * weight1 + epie2 * weight2 ) / (weight1 + weight2);
      ep_em_rat   = ( epem1 * weight1 + epem2 * weight2 ) / (weight1 + weight2);
      ep_em_ratErr= ( eepem1 * weight1 + eepem2 * weight2)/ (weight1 + weight2);
    }

  if(db == 1 || db == 2)
    {
      printf("BoNuS contamination ratios for energy = %.0f MeV:\n", energy);
      printf("Lower EG1B Energy = %.0f, weight = %f\n",e1,weight1);
      printf("pi/e = %f +- %f\n", pie1,epie1);
      printf("e+/e- = %f +- %f\n\n", epem1,eepem1);
      printf("Upper EG1B Energy = %.0f, weight = %f\n",e2,weight2);
      printf("pi/e = %f +- %f\n", pie2,epie2);
      printf("e+/e- = %f +- %f\n\n", epem2,eepem2);
      printf("Final Value for pi/e = %f +- %f\n", pi_e_rat,pi_e_ratErr);
      printf("Final Value for e+/e- = %f +- %f\n\n", ep_em_rat,ep_em_ratErr);      
    }

}



#endif //// __Conta_h__