#include<TVector3.h>

Double_t Ecorr_func (Double_t *x, Double_t *par){
    
    // parameterization of energy correction function: E_corr = E/Ecoor;
    
    // constants
    Double_t q0 = par[0];
    Double_t q1 = par[1];
    Double_t q2 = par[2];
    Double_t q3 = par[3];
    Double_t q4 = par[4];
    Double_t q5 = par[5];
    Double_t q6 = par[6];
    Double_t q7 = par[7];
    Double_t q8 = par[8];
    
    Double_t Eg = x[0];
    Double_t f = 0;
    
    f = Eg>0? q0 - q1*exp(q2*Eg+q3) - q4/(q5 + q6*exp(-q7*Eg+q8)) : 0;
     
   // cout << " Ecorr parameters=" << q0 << " " << q1 << " " << q2 << " " << q3 << " " <<  q4 << " " << q5 << " " << q6 << " " << q7 << " " << q8 << endl;
   // cout << "Ecorr_func: Eg=" << Eg << " 1/f=" << 1/f << endl;
    
    return f;
}
Double_t leakage_func (Double_t *x, Double_t *par){
    
    // leakage parameterization based on BCAL NIM paper Fig. 10
    // leakage as a fraction
    
    // constants
    Double_t k1 = par[0];
    Double_t k2 = par[1];
    Double_t k3 = par[2];
    // Double_t ktheta = par[3];
    Double_t ktheta = 90.;          // eliminate theta dependence.
    Double_t k4 = par[4];
    Double_t k5 = par[5];
    
    Double_t Eg = x[0];
    Double_t f1=0;
    Double_t f2=0;
    
    Double_t factor=1;         // decrease angle dependence; nominal factor=1
    Double_t lambda=factor*25;
    Double_t k1m = k1*exp(-(ktheta-12)/lambda);
    lambda=50*factor;
    Double_t factor2=1;       // decrease angle dependence; nominal factor=1
    Double_t k2m = factor2*k2*exp(-(ktheta-12)/lambda);
    lambda = 30*factor;
    Double_t k3m = k3*exp(-(ktheta-12)/lambda);
    
    
    // Front leakage
    f1 = Eg>0? k1m + k2m*exp(-Eg/k3m): 0;   // return zero for negative Eg
    // f1 = Eg>0? 0.2/Eg: 0;                   // fixed leakage
    
    lambda = 30;
    Double_t k4m = k4*exp(-(90-ktheta)/lambda);
    Double_t k5m = k5;
    
    // Back leakage
    // f2 = Eg>0? k4m*(1-exp(-Eg/k5m)) : 0;     // constant at high energy
    f2 = Eg>0? k4m*sqrt(Eg/k5m) : 0;
    // f2 = 0;                                 // no leakage
    
    return f1+f2;
}

Double_t eff_func (Double_t *x, Double_t *par){
    
    // Dummy efficiency function taken from Qiao's talk GlueX-doc-3844 for ease of comparison
    // f = 0.4 + sqrt(Eg);
    
    // constants
    // Double_t const PI = 3.14159;
    Double_t a = par[0];
    Double_t b = par[1];
    
    Double_t MPI =0.139570;
    Double_t PI=3.14159;
    
    Double_t Eg = x[0];
    Double_t f=0;
    
    f = Eg>0? a + b*sqrt(Eg): 0;   // return zero for negative Eg
    if (f > 1) f = 1;             // cap efficiency at 1
    
    f=1;                          // 100%
    return f;
}
Double_t tanLabCM_func (Double_t *x, Double_t *par){
    
    // Tan (angle_Lab) given the angle_Lab (degrees) in the CM system
    // betastar: Velocity (units of c) of CM relative to Lab (i.e. negative)
    // betaCM: Velocity of particle (units of c) in CM
    
    // constants
    // Double_t const PI = 3.14159;
    Double_t betastar= par[0];
    Double_t betaCM= par[1];
    
    Double_t MPI =0.139570;
    Double_t PI=3.14159;
    
    Double_t thetaCM = x[0]*PI/180.;
    
    Double_t tanLab=0;
    
    if (betastar < -1 || betastar >1) {
        cout << " ***tanLabCM_func: betastar=" << betastar << endl;
        return tanLab;
    }
    Double_t gammastar =  1/sqrt(1-betastar*betastar);
    
    tanLab = (1/gammastar)*(sin(thetaCM)/(cos(thetaCM) + betastar/betaCM));
    Double_t thetaLab = atan(tanLab)*180./PI;
    
    // cout << " betastar=" << betastar << " betaCM=" << betaCM << " thetaCM=" << thetaCM << " tanLab=" << tanLab << " thetalab=" << thetaLab << endl;
    
    return	thetaLab>0? thetaLab : thetaLab + 180.;
}
void NonLinCheck(void)
{
// Consider transformation of angles from pi0->gg between CM and lab for the study of efficiency calculations in the calorimeters
    // Dec 17, 2018.
//

  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit(111111);
  gStyle->SetPadRightMargin(0.15);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadBottomMargin(0.15);
    
    
    char string[256];
    TString tag="Test";
    TString filename = "NonLinCheck_"+tag;
    // define function
    
    // Define historgrams
    Int_t nbins=20;
    
    // initialize random number seeds
    
    UInt_t  iseed = 2342411151;
    TRandom1 *r1 = new TRandom1();
    r1->SetSeed(iseed);
    iseed = 2342411155;
    TRandom1 *r2 = new TRandom1();
    r2->SetSeed(iseed);
    iseed = 2342411157;
    TRandom1 *r3 = new TRandom1();
    r3->SetSeed(iseed);
    iseed = 2342411153;
    TRandom1 *r4 = new TRandom1();
    r4->SetSeed(iseed);
    iseed = 2342411353;
    TRandom1 *r5 = new TRandom1();
    r5->SetSeed(iseed);
    iseed = 2342411333;
    TRandom1 *r6 = new TRandom1();
    r6->SetSeed(iseed);
    iseed = 2342411555;
    TRandom1 *r7 = new TRandom1();
    r7->SetSeed(iseed);
    
    
    // kinematics for Ma -> Mb + X (for pi0 Mpi0 -> g g)

    Double_t PI=3.14159;
    Double_t MPI =0.139570;
    // Double_t MPI =0.8;
    Double_t paLab0 =8.;
    Double_t Mg = 0;
    Double_t Ma = MPI;
    Double_t Mb = Mg;
    Double_t EaLab0 = sqrt(paLab0*paLab0 + Ma*Ma);
    Double_t betastar = -paLab0/EaLab0;
    
    Double_t EbCM = Ma/2;
    Double_t pbCM = sqrt(EbCM*EbCM - Mb*Mb);
    Double_t betabCM = pbCM/EbCM;
    Double_t Ethresh = 0.5;      // 50 Mev Threshold
    Double_t BCALa=0.052;
    Double_t BCALb=0.036;
    /*Double_t BCALa=0.01;
    Double_t BCALb=0.01;
    Double_t Ethresh = 0.001;   // 1 Mev Threshold*/
    
    const Int_t npts=1000000;
    // const Int_t npts=10000000;
    
    Double_t xmin=0;
    Double_t xmax=2;
    Double_t aff=0.4;
    Double_t beff=1;
    Int_t npar = 2;
    
    TF1 *eff = new TF1("eff",eff_func,xmin,xmax,npar);
    eff->SetParameters(aff,beff);
    eff->SetParNames("a","b");
    
    npar = 9;
    Double_t q0=0.97;
    Double_t q1=0;
    Double_t q2=0;
    Double_t q3=0;
    Double_t q4=0;
    Double_t q5=1;
    Double_t q6=1;
    Double_t q7=0;
    Double_t q8=0;
    
    xmin=0;
    xmax=paLab0;
    TF1 *Ecorr = new TF1("Ecorr",Ecorr_func,xmin,xmax,npar);
    Ecorr->SetParameters(q0,q1,q2,q3,q4,q5,q6,q7,q8);
    Ecorr->SetParNames("q0","q1","q2","q3","q4","q5","q6","q7","q8");
    
    Double_t k1=0.06;
    Double_t k2=0.1;
    // Double_t k3=0.5;    // nominal value
    Double_t k3=4;     // increase energy dependence
    Double_t ktheta=90;
    /* Double_t k4=0.03;       // constants for exponential fit
    Double_t k5=0.3;*/
    
    Double_t k4=0.03;       // constants for exponential fit
    Double_t k5=1;

    npar = 6;
    xmin=0;
    xmax=10;
    TF1 *leakage = new TF1("leakage",leakage_func,xmin,xmax,npar);
    leakage->SetParameters(k1,k2,k3,ktheta,k4,k5);
    leakage->SetParNames("k1","k2","k3","ktheta","k4","k5");
    
    TF1 *inverse = new TF1("inverse","1 + 0.5*[0]*(x+[1])/(x*[1])",0,paLab0);
    inverse->SetParameters(0.2,1);
    inverse->SetParNames("leak","E2");
    
    TF1 *inverse2 = new TF1("inverse2","[1] + 0.5*[0]*(x+[2])",0,paLab0);
    inverse2->SetParameters(0.2,1,1);
    inverse2->FixParameter(2,1);   // fix E2=1
    inverse2->SetParNames("leak","gain","E2");
    
    
    // fits to leakage fraction ~ 0.2/E. These should match at f=4.
    /*root [0] TF1 *e1 = new TF1("e1","log( ([1]-x)/[0]) + ([1] - x - [0])/x",1,4);
    root [1] TF1 *e2 = new TF1("e2","1 + log(x/[0]) - 2*[0]/[1]",4,7);
    root [2] e1->SetParameter(0,0.5);
    root [3] e1->SetParameter(1,8);
    root [4] e2->SetParameter(0,0.5);
    root [5] e2->SetParameter(1,8);*/

    TCanvas *c0 = new TCanvas("c0", "c0",200,10,1000,1000);
    c0->Divide(2,3);
    
    c0->cd(1);
    ktheta = 12.;
    TF1 *leakage_draw = new TF1("leakage_draw",leakage_func,xmin,xmax,npar);
    leakage_draw->SetParameters(k1,k2,k3,ktheta,k4,k5);
    leakage_draw->SetParNames("k1","k2","k3","ktheta","k4","k5");
    leakage_draw->Draw();
    

    // Define histograms
    
    nbins=150;
    TH1D *h1_costhestar = new TH1D ("h1_costhestar","Cos(theta_CM)",nbins,-1.5,1.5);
    TH1D *h1_costhestar_eff = new TH1D ("h1_costhestar_eff","Efficiency vs Cos(theta_CM)",nbins,-1.5,1.5);   // must have same binning as previous
    
    TH1D *h1_costhe = new TH1D ("h1_costhe","Cos(theta_Lab)",nbins,-1.5,1.5);
    TH1D *h1_DeltaCos = new TH1D ("h1_DeltaCos","CostheCM - CostheCMMeas",nbins,-1.5,1.5);
    TH1D *h1_costhestar_meas = new TH1D ("h1_costhestar_meas","Cos(theta_CM) Measured",nbins,-1.5,1.5);
    TH2D *h2_Cos_measVSCos = new TH2D("h2_Cos_measVSCos","Cos Meas vs Cos",nbins,-1.5,1.5,nbins,-1.5,1.5);
    TH1D *h1_E1 = new TH1D("h1_E1","E1 (GeV)",100,0,paLab0);
    TH1D *h1_E1meas = new TH1D("h1_E1meas","E1 meas (GeV)",100,0,paLab0);
    TH1D *h1_E2meas = new TH1D("h1_E2meas","E2 meas (GeV)",100,0,paLab0);
    TH1D *h1_Emeas = new TH1D("h1_Emeas","E1+E2 meas (GeV)",200,0,paLab0);
    TH1D *h1_the1Lab = new TH1D("h1_the1Lab","the1Lab (deg)",100,0,100);
    TH1D *h1_the2Lab = new TH1D("h1_the2Lab","the2Lab (deg)",100,0,100);
    TH1D *h1_patheta = new TH1D("h1_patheta","patheta (deg)",100,0,100);
    TH1D *h1_E1diff = new TH1D("h1_E1diff","(E1 - E1meas)/E1 (GeV)",100,-0.3,0.3);
    TH1D *h1_E2diff = new TH1D("h1_E2diff","(E2 - E2meas)/E2 (GeV)",100,-0.3,0.3);
    Double_t ratio_min=0.7;
    Double_t ratio_max=1.3;
    TH1D *h1_pi0ratio = new TH1D("h1_pi0ratio","pi0meas/pi0",100,ratio_min,ratio_max);
    TH1D *h1_pi0ratio_corr1 = new TH1D("h1_pi0ratio_corr1","pi0meas/pi0 (E2_corr)",100,ratio_min,ratio_max);
    TH1D *h1_pi0ratio_corr2 = new TH1D("h1_pi0ratio_corr2","pi0meas/pi0 (E2_corr and E1_corr)",100,ratio_min,ratio_max);
    TH1D *h1_pi0ratio_DE2 = new TH1D("h1_pi0ratio_DE2","pi0meas/pi0 (DE2<0.5 GeV)",100,ratio_min,ratio_max);
    TH2D *h2_pi0ratioVSE1 = new TH2D ("h2_pi0ratioVSE1","pi0mratio vs E1 (GeV)",100,0,paLab0+2,100,ratio_min,ratio_max);
    TH2D *h2_pi0ratioVSE1_corr1 = new TH2D ("h2_pi0ratioVSE1_corr1","pi0mratio (E2_corr) vs E1 (GeV)",100,0,paLab0+2,100,ratio_min,ratio_max);
    TH2D *h2_pi0ratioVSE1_corr2 = new TH2D ("h2_pi0ratioVSE1_corr2","pi0mratio (E2_corr and E1_corr) vs E1 (GeV)",100,0,paLab0+2,100,ratio_min,ratio_max);
    TH2D *h2_pi0ratioVSE1_DE2 = new TH2D ("h2_pi0ratioVSE1_DE2","pi0mratio (DE2<0.5 GeV) vs E1 (GeV)",100,0,paLab0+2,100,ratio_min,ratio_max);
    
    TH2D *h2_costheVScosthestar = new TH2D ("h2_costheVScosthestar","Cos(theta_Lab) vs Cos(theta_CM)",100,-1,1,100,-1,1);
    TH2D *h2_E2VSE1 = new TH2D ("h2_E2VSE1","E2 (GeV) vs E1 (GeV)",100,0,paLab0,100,0,paLab0);
    TH2D *h2_E2measVSE1meas = new TH2D ("h2_E2measVSE1meas","E2meas (GeV) vs E1meas (GeV)",100,0,paLab0,100,0,paLab0);
    // generate photon angular distribution in the CM, then boost to the lab.
    TH2D *h2_the1LabVSpatheta = new TH2D ("h2_the1LabVSpatheta","the1Lab VS patheta (deg)",100,0,100,100,0,100);
    TH2D *h2_the2LabVSpatheta = new TH2D ("h2_the2LabVSpatheta","the2Lab VS patheta (deg)",100,0,100,100,0,100);
    TH2D *h2_the2LabVSthe1Lab = new TH2D ("h2_the2LabVSthe2Lab","the2Lab VS the1Lab (deg)",100,0,100,100,0,100);
    
    
    // Open input file
    ifstream parameters;
    Int_t iter=1;
    TString infile;
    infile.Form(filename+"_iter%d.in",iter);
    parameters.open (infile.Data());
    if (!parameters) {
        cout << "ERROR: Failed to open data file= " << infile.Data() << endl;
        return 1;   // failed to open data file
    }
    else {
        cout << "Openned input file=" << infile << endl;
    }
    
    TString line;
    while (line.ReadLine(parameters)){
        
        TObjArray *tokens = line.Tokenize(" ");
        Int_t ntokens = tokens->GetEntries();
        q0 = (((TObjString*)tokens->At(0))->GetString()).Atof();
        q1 = (((TObjString*)tokens->At(1))->GetString()).Atof();
        q2 = (((TObjString*)tokens->At(2))->GetString()).Atof();
        q3 = (((TObjString*)tokens->At(3))->GetString()).Atof();
        q4 = (((TObjString*)tokens->At(4))->GetString()).Atof();
        q5 = (((TObjString*)tokens->At(5))->GetString()).Atof();
        q6 = (((TObjString*)tokens->At(6))->GetString()).Atof();
        q7 = (((TObjString*)tokens->At(7))->GetString()).Atof();
        q8 = (((TObjString*)tokens->At(8))->GetString()).Atof();
        
        
    }

    
    cout << " Input parameters=" << q0 << " " << q1 << " " << q2 << " " << q3 << " " <<  q4 << " " << q5 << " " << q6 << " " << q7 << " " << q8 << endl;

    
    Ecorr->SetParameters(q0,q1,q2,q3,q4,q5,q6,q7,q8);  // set parameters to input file

    // open list output file

    ofstream outlist;
    TString ofilename;
    ofilename.Form("%s_iter%d.out",filename.Data(),iter+1);
    outlist.open (ofilename.Data());
    
    Double_t Estar = Ma/2;
    Double_t thetaCMmin=0;
    Double_t thetaCMmax=180;
    Double_t cosCMmax=cos(thetaCMmin*PI/180.);
    Double_t cosCMmin=cos(thetaCMmax*PI/180.);
    cout << " thetaCMmin=" << thetaCMmin << " thetaCMmax=" << thetaCMmax << " cosCMmin=" << cosCMmin << " cosCMmax=" << cosCMmax << endl;
    Double_t EaLabmin=0;
    
    Double_t thetamin=10;
    Double_t thetamax=90;
    Double_t cosmax=cos(thetamin*PI/180.);
    Double_t cosmin=cos(thetamax*PI/180.);
    cout << " thetamin=" << thetamin << " thetamax=" << thetamax << " cosmin=" << cosmin << " cosmax=" << cosmax << endl;

    Double_t DeltaE1_min=2;
    Double_t DeltaE1_max=2.5;
    Double_t DeltaE2_min=0.95;
    Double_t DeltaE2_max=1.05;
    // Double_t DeltaE2_min=0.5;
    // Double_t DeltaE2_max=4;

    for (Int_t jj=0; jj<npts; jj++) {
        Double_t EaLab = EaLabmin + r4->Uniform()*(EaLab0-EaLabmin);     // pick random Energy between 0 and EaLab0
        // Double_t sigmaE = EaLab*sqrt(BCALa*BCALa/EaLab+BCALb*BCALb);   // compute
        Double_t paLab = sqrt(EaLab*EaLab - Ma*Ma);
        // Double_t pacosthe = cosmin + (cosmax-cosmin)*r6->Uniform();
        // Double_t patheta = acos(pacosthe)*180./PI;
        Double_t patheta = thetamin + (thetamax-thetamin)*r6->Uniform();
        Double_t pacosthe = cos(patheta*PI/180.);
        // cout << " patheta=" << patheta << endl;;
        // Double_t paphi = 360*r3->Uniform();
        Double_t paphi = 360*r3->Uniform();
        Double_t pasinthe = sin(patheta*PI/180.);
        Double_t papx = paLab*pasinthe*cos(paphi*PI/180.);
        Double_t papy = paLab*pasinthe*sin(paphi*PI/180.);
        Double_t papz = paLab*pacosthe;
        TLorentzVector pi0lab (papx,papy,papz,EaLab);
        TLorentzVector pi0boost (0,0,paLab,EaLab);    // want to boost along the pi0 direction
        
        Double_t costhestar = cosCMmin + (cosCMmax-cosCMmin)*r1->Uniform();
        Double_t thetastar = acos(costhestar)*180./PI;
        // Double_t phistar = 360*r5->Uniform();
        Double_t phistar = 360*r5->Uniform();
        Double_t sinthestar = sin(thetastar*PI/180.);
        Double_t pxstar = Estar*sinthestar*cos(phistar*PI/180.);
        Double_t pystar = Estar*sinthestar*sin(phistar*PI/180.);
        Double_t pzstar = Estar*costhestar;
        TLorentzVector photonstar1(pxstar,pystar,pzstar,Estar);
        TLorentzVector photonstar2(-pxstar,-pystar,-pzstar,Estar);
        TLorentzRotation Boost(pi0boost.BoostVector() );
        TLorentzVector photon1 = Boost * photonstar1;
        TLorentzVector photon2 = Boost * photonstar2;
        TLorentzVector photonsum = photon1 + photon2;
        Double_t thetasum = photonsum.Theta()*180./PI;
        
        /*cout << " photonstar1="; photonstar1.Print();
        cout << " photonstar2="; photonstar2.Print();
        cout << " photon1="; photon1.Print();
        cout << " photon2="; photon2.Print();
        cout << endl;*/
        
             TVector3 zhatp = pi0lab.Vect().Unit();
             TVector3 zhat(0,0,1);
             TVector3 xhat(1,0,0);
             TVector3 yhat(0,1,0);
             TVector3 yhatp = zhat.Cross(zhatp).Unit();
             TVector3 xhatp = yhatp.Cross(zhatp);
             
             /*cout << " xhat="; xhat.Print(); cout << " yhat="; yhat.Print(); cout << " zhat="; zhat.Print(); cout << endl;
             cout << " xhapt="; xhatp.Print(); cout << " yhatp="; yhatp.Print(); cout << " zhatp="; zhatp.Print(); cout << endl;
             // cout << " xhat=" << xhat.Print() << endl;
             // << " yhat=" << yhat.Print() << " zhat=" << zhat.Print() << endl;*/
             
             TLorentzVector photon1Lab(photon1.Vect().Dot(xhatp),photon1.Vect().Dot(yhatp),photon1.Vect().Dot(zhatp),photon1.E());
             TLorentzVector photon2Lab(photon2.Vect().Dot(xhatp),photon2.Vect().Dot(yhatp),photon2.Vect().Dot(zhatp),photon2.E());
             TLorentzVector photonsumLab = photon1Lab + photon2Lab;
             Double_t thetasumLab = photonsumLab.Theta()*180./PI;
             Double_t the1Lab = photon1Lab.Theta()*180./PI;
             Double_t costhe1Lab = cos(photon1Lab.Theta());
             Double_t the2Lab = photon2Lab.Theta()*180./PI;
             Double_t costhe2Lab = cos(photon2Lab.Theta());
             // cout << " thetasumLab=" << thetasumLab << " the1Lab=" << the1Lab << " the2Lab=" << the2Lab << endl << endl;
        
        Double_t E1 = photon1.E();
        Double_t E2 = photon2.E();
        
        // Here we add leakage based on a simplified model
        
        leakage->SetParameters(k1,k2,k3,the1Lab,k4,k5);
        Double_t genleak;
        genleak=r7->Exp(leakage->Eval(E1));            // randomize leakage
        // genleak = leakage->Eval(E1);
        Double_t E1_miss = genleak<0.9? E1*genleak : E1*0.9;    // limit leakage to 90% of value to avoid negative energies
        
        leakage->SetParameters(k1,k2,k3,the2Lab,k4,k5);
        // genleak = r7->Exp(leakage->Eval(E2));         // randomize leakage
        genleak = leakage->Eval(E2);
        Double_t E2_miss = genleak<0.9? E2*genleak : E2*0.9;
        
        Double_t sigma = E1*sqrt(BCALa*BCALa/E1+BCALb*BCALb);
        Double_t genres;
        genres = r2->Gaus(E1-E1_miss,sigma);
        Double_t E1meas = genres>0? genres: 0.;
        sigma = E2*sqrt(BCALa*BCALa/E2+BCALb*BCALb);
        genres = r2->Gaus(E2-E2_miss,sigma);
        Double_t E2meas = genres>0? genres: 0.;
        
        // cout << " E1=" << E1 << " ktheta=" << the1Lab << " E_miss1=" << E1_miss << " E2_miss=" << E2_miss << endl;
        
        if (E1meas < E2meas) {    // always take E1meas to have the higher energy. Switch both gen and meas energies
                   TLorentzVector photontemp = photon1;
                   photon1 = photon2;
                   photon2 = photontemp;
                   Double_t Etemp = E1;
                   E1 = E2;
                   E2 = Etemp;
                   Etemp = E1meas;
                   E1meas = E2meas;
                   E2meas = Etemp;
                   costhestar = -costhestar;
               }
        // Note that variables may be overwritten depending of selection of E1 > E2
        
        Double_t the1 = photon1.Theta()*180./PI;
        Double_t costhe1 = cos(photon1.Theta());
        Double_t the2 = photon2.Theta()*180./PI;
        Double_t costhe2 = cos(photon2.Theta());
        the1Lab = photon1Lab.Theta()*180./PI;
        costhe1Lab = cos(photon1Lab.Theta());
        the2Lab = photon2Lab.Theta()*180./PI;
        costhe2Lab = cos(photon2Lab.Theta());
        
        /*cout << " pi0lab:"; pi0lab.Print();
        cout << " photonsum:"; photonsum.Print();
        cout << " photonstar1:"; photonstar1.Print();
        cout << " photonstar2:"; photonstar2.Print();
        cout << " photon1:";  photon1.Print();
        cout << " photon2:"; photon2.Print();
        cout << " EaLab=" << EaLab << " E1+E2=" << E1+E2 << " patheta=" << patheta << " thetasum=" << thetasum << " the1=" << the1 << " the2=" << the2 << endl;
        */
        
        /*cout << " photon1Lab:";  photon1Lab.Print();
        cout << " photon2Lab:"; photon2Lab.Print();
        cout << endl;*/
        Double_t ratio = E1meas/E1;
        TLorentzVector photon1meas(photon1Lab.Vect().X()*ratio,photon1Lab.Vect().Y()*ratio,photon1Lab.Vect().Z()*ratio,photon1Lab.E()*ratio);
        ratio = E2meas/E2;
        TLorentzVector photon2meas(photon2Lab.Vect().X()*ratio,photon2Lab.Vect().Y()*ratio,photon2Lab.Vect().Z()*ratio,photon2Lab.E()*ratio);
        TLorentzVector pi0meas = photon1meas + photon2meas;
        Double_t pi0mass_ratio = pi0meas.M()/photonsumLab.M();
        Double_t nonlin_corr1 = 1./Ecorr->Eval(E1meas);
        Double_t nonlin_corr2 = 1./Ecorr->Eval(E2meas);
        
        TLorentzVector photon1meas_corr = nonlin_corr1*photon1meas;
        TLorentzVector photon2meas_corr = nonlin_corr2*photon2meas;
        TLorentzVector pi0meas_corr1 = photon1meas_corr + photon2meas;
        Double_t pi0mass_ratio_corr1 = pi0meas_corr1.M()/photonsumLab.M();
        TLorentzVector pi0meas_corr2 = photon1meas_corr + photon2meas_corr;
        Double_t pi0mass_ratio_corr2 = pi0meas_corr2.M()/photonsumLab.M();

        /*if (E1meas>DeltaE1_min && E1meas<DeltaE1_max) {
            cout << " E1=" << E1 << " E2=" << E2 << " E1meas=" << E1meas << " E2meas=" << E2meas << " E1_miss=" << E1_miss << " E2_miss=" << E2_miss << " Leakage1=" << E1_miss/E1 << " Leakage2=" << E2_miss/E2 << " pi0mass_ratio=" << pi0mass_ratio << endl;
        }*/
        
        // cout << " photonsumLab.M=" << photonsumLab.M() << " pi0meas.M=" << pi0meas.M() << " ratio=" << pi0mass_ratio << endl;
        
        if (E2meas < Ethresh) continue;     // impose threshold on smallest energy
        if (E1meas < Ethresh) continue;     // impose threshold on largest energy (redundant except for testing)
        if (r4->Uniform() > eff->Eval(E1meas)) continue;   // efficiency for detecting E1.
        if (r4->Uniform() > eff->Eval(E2meas)) continue;   // efficiency for detecting E2.
        if (E2meas>DeltaE2_max || E2meas<DeltaE2_min) continue;   // impose condition on E2 for all analyzes
        
        // cout << " E1meas=" << E1meas << " nonlin_corr1=" << nonlin_corr1 << " E2meas=" << E2meas << " nonlin_corr2=" << nonlin_corr2 << " mass_corr1=" << pi0mass_ratio_corr1 << " mass_corr2=" << pi0mass_ratio_corr2 << endl;

        h1_E1->Fill(E1);
        h1_E1meas->Fill(E1meas);
        h1_Emeas->Fill(E1meas+E2meas);
        h2_E2VSE1->Fill(E1,E2);
        h1_costhestar->Fill(costhestar);
        h1_costhe->Fill(costhe1);
        h2_costheVScosthestar->Fill(costhestar,costhe1);
        if (E2meas<DeltaE2_max && E2meas>DeltaE2_min) h2_E2measVSE1meas->Fill(E1meas,E2meas);
  
        Double_t costhestar2 = (E1meas-E2meas)/sqrt( (E1meas+E2meas)*(E1meas+E2meas) - Ma*Ma);
        
        // if (costhestar2 > 1) costhestar2 = 1;     // impose physicsl limits
        // if (costhestar2 < -1) costhestar2 = -1;
        h1_costhestar_meas->Fill(costhestar2);
        h1_DeltaCos->Fill(costhestar2 - costhestar);
        h2_Cos_measVSCos->Fill(costhestar,costhestar2);
        // cout << " costhestar=" << costhestar << " costhestar2=" << costhestar2 << endl << endl;
        
        h1_the1Lab->Fill(the1Lab);
        h1_the2Lab->Fill(the2Lab);
        h1_patheta->Fill(patheta);
        h2_the1LabVSpatheta->Fill(patheta,the1Lab);
        h2_the2LabVSpatheta->Fill(patheta,the2Lab);
        h2_the2LabVSthe1Lab->Fill(the1Lab,the2Lab);
        
        h1_E1diff->Fill((E1-E1meas)/E1);
        h1_E2diff->Fill((E2-E2meas)/E2);
        h1_pi0ratio->Fill(pi0mass_ratio);
        if (E2meas<DeltaE2_max && E2meas>DeltaE2_min) h1_pi0ratio_DE2->Fill(pi0mass_ratio);
        h2_pi0ratioVSE1->Fill(E1meas,pi0mass_ratio);
        if (E2meas<DeltaE2_max && E2meas>DeltaE2_min) h2_pi0ratioVSE1_DE2->Fill(E1meas,pi0mass_ratio);
        h1_pi0ratio_corr1->Fill(pi0mass_ratio_corr1);
        h2_pi0ratioVSE1_corr1->Fill(E1meas,pi0mass_ratio_corr1);
        h1_pi0ratio_corr2->Fill(pi0mass_ratio_corr2);
        h2_pi0ratioVSE1_corr2->Fill(E1meas,pi0mass_ratio_corr2);
        
        
        // cout << " jj=" << jj << " costhestar=" << costhestar << endl;
    }
    
    TString text;
    
    TCanvas *c2 = new TCanvas("c2", "c2",200,10,1000,1000);
    c2->Divide(2,3);
    
    c2->cd(1);
    Double_t ymin=0.1;
    Double_t ymax=h1_costhestar->GetMaximum()*1.5;
    gPad->SetLogy();
    text.Form("p_{#pi}=%.2f, E_{thresh}=%.2f, a=%.3f, b=%.3f, iter=%d",paLab0,Ethresh,BCALa,BCALb,0);
    h1_costhestar->SetTitle(text);
    // h1_costhestar->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhestar->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhestar->GetXaxis()->SetTitleSize(0.05);
    h1_costhestar->GetYaxis()->SetTitleSize(0.05);
    h1_costhestar->GetXaxis()->SetTitle("Cos(theta_CM)");
    h1_costhestar->GetYaxis()->SetTitle("Events");
    h1_costhestar->SetLineColor(1);
    h1_costhestar->Draw();
    
    ymin=0.1;
    ymax=h1_costhe->GetMaximum()*1.5;
    c2->cd(2);
    h1_costhe->SetTitle(text);
    gPad->SetLogy();
    // h1_costhe->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhe->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhe->GetXaxis()->SetTitleSize(0.05);
    h1_costhe->GetYaxis()->SetTitleSize(0.05);
    h1_costhe->GetXaxis()->SetTitle("Cos(theta_Lab)");
    h1_costhe->GetYaxis()->SetTitle("Events");
    h1_costhe->SetLineColor(1);
    h1_costhe->Draw();
    
    c2->cd(3);
    h2_costheVScosthestar->SetTitle(text);
    // h2_costheVScosthestar->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_costheVScosthestar->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_costheVScosthestar->GetXaxis()->SetTitleSize(0.05);
    h2_costheVScosthestar->GetYaxis()->SetTitleSize(0.05);
    h2_costheVScosthestar->GetXaxis()->SetTitle("Cos(theta_CM)");
    h2_costheVScosthestar->GetYaxis()->SetTitle("Cos(theta_Lab)");
    h2_costheVScosthestar->SetLineColor(1);
    h2_costheVScosthestar->Draw("colz");
    
    c2->cd(4);
    h2_E2VSE1->SetTitle(text);
    // h2_E2VSE1->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_E2VSE1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_E2VSE1->GetXaxis()->SetTitleSize(0.05);
    h2_E2VSE1->GetYaxis()->SetTitleSize(0.05);
    h2_E2VSE1->GetXaxis()->SetTitle("E1 (GeV)");
    h2_E2VSE1->GetYaxis()->SetTitle("E2 (GeV)");
    h2_E2VSE1->SetLineColor(1);
    h2_E2VSE1->Draw("colz");
    
    c2->cd(5);
    h2_E2measVSE1meas->SetTitle(text);
    // h2_E2measVSE1meas->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_E2measVSE1meas->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_E2measVSE1meas->GetXaxis()->SetTitleSize(0.05);
    h2_E2measVSE1meas->GetYaxis()->SetTitleSize(0.05);
    h2_E2measVSE1meas->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_E2measVSE1meas->GetYaxis()->SetTitle("E2meas (GeV)");
    h2_E2measVSE1meas->SetLineColor(1);
    h2_E2measVSE1meas->Draw("colz");
                          
    c2->cd(6);
    h1_DeltaCos->SetTitle(text);
	// h1_E2measVSE1meas->GetXaxis()->SetRangeUser(xmin,xmax);
    // h1_DeltaCos->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_DeltaCos->GetXaxis()->SetTitleSize(0.05);
    h1_DeltaCos->GetYaxis()->SetTitleSize(0.05);
    h1_DeltaCos->GetXaxis()->SetTitle("CosTheCM_meas - CosTheCM");
    h1_DeltaCos->GetYaxis()->SetTitle("Events");
    h1_DeltaCos->SetLineColor(1);
    h1_DeltaCos->Fit("gaus");
    h1_DeltaCos->Draw();
    
    TCanvas *c3 = new TCanvas("c3", "c3",200,10,1000,1000);
    c3->Divide(2,3);
    
    ymin=0;
    ymax=h1_costhestar_meas->GetMaximum()*1.5;
    c3->cd(1);
    h1_costhestar_meas->SetTitle(text);
    // h1_costhestar_meas->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhestar_meas->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhestar_meas->GetXaxis()->SetTitleSize(0.05);
    h1_costhestar_meas->GetYaxis()->SetTitleSize(0.05);
    h1_costhestar_meas->GetXaxis()->SetTitle("Cos(theta_CM) meas");
    h1_costhestar_meas->GetYaxis()->SetTitle("Events");
    h1_costhestar_meas->SetLineColor(1);
    h1_costhestar_meas->Draw();
    
    c3->cd(2);
    h2_Cos_measVSCos->SetTitle(text);
    // h2_Cos_measVSCos->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_Cos_measVSCos->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_Cos_measVSCos->GetXaxis()->SetTitleSize(0.05);
    h2_Cos_measVSCos->GetYaxis()->SetTitleSize(0.05);
    h2_Cos_measVSCos->GetXaxis()->SetTitle("Cos(theta_CM)");
    h2_Cos_measVSCos->GetYaxis()->SetTitle("Cos(theta_CM) meas");
    h2_Cos_measVSCos->SetLineColor(1);
    h2_Cos_measVSCos->Draw("colz");
    
    ymin=0;
    ymax=h1_E1->GetMaximum()*1.5;
    c3->cd(3);
    h1_E1->SetTitle(text);
    // h1_E1->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_E1->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_E1->GetXaxis()->SetTitleSize(0.05);
    h1_E1->GetYaxis()->SetTitleSize(0.05);
    h1_E1->GetXaxis()->SetTitle("E1 (GeV)");
    h1_E1->GetYaxis()->SetTitle("Events");
    h1_E1->SetLineColor(1);
    h1_E1->Draw();
    
    ymin=0;
    ymax=h1_E1meas->GetMaximum()*1.5;
    c3->cd(4);
    h1_E1meas->SetTitle(text);
    // h1_E1meas->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_E1meas->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_E1meas->GetXaxis()->SetTitleSize(0.05);
    h1_E1meas->GetYaxis()->SetTitleSize(0.05);
    h1_E1meas->GetXaxis()->SetTitle("E1meas (GeV)");
    h1_E1meas->GetYaxis()->SetTitle("Events");
    h1_E1meas->SetLineColor(1);
    h1_E1meas->Draw();
    
    ymin=0;
    ymax=h1_the1Lab->GetMaximum()*1.5;
    c3->cd(5);
    h1_the1Lab->SetTitle(text);
    // h1_the1Lab->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_the1Lab->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_the1Lab->GetXaxis()->SetTitleSize(0.05);
    h1_the1Lab->GetYaxis()->SetTitleSize(0.05);
    h1_the1Lab->GetXaxis()->SetTitle("the1Lab (GeV)");
    h1_the1Lab->GetYaxis()->SetTitle("Events");
    h1_the1Lab->SetLineColor(1);
    h1_the1Lab->Draw();
    
    ymin=0;
    ymax=h1_patheta->GetMaximum()*1.5;
    c3->cd(6);
    h1_patheta->SetTitle(text);
    // h1_patheta->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_patheta->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_patheta->GetXaxis()->SetTitleSize(0.05);
    h1_patheta->GetYaxis()->SetTitleSize(0.05);
    h1_patheta->GetXaxis()->SetTitle("patheta (GeV)");
    h1_patheta->GetYaxis()->SetTitle("Events");
    h1_patheta->SetLineColor(1);
    h1_patheta->Draw();
    

    TCanvas *c4 = new TCanvas("c4", "c4",200,10,1000,1000);
    c4->Divide(2,3);
    
    ymin=0;
    ymax=h1_the2Lab->GetMaximum()*1.5;
    c4->cd(1);
    h1_the2Lab->SetTitle(text);
    // h1_the2Lab->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_the2Lab->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_the2Lab->GetXaxis()->SetTitleSize(0.05);
    h1_the2Lab->GetYaxis()->SetTitleSize(0.05);
    h1_the2Lab->GetXaxis()->SetTitle("the2Lab (GeV)");
    h1_the2Lab->GetYaxis()->SetTitle("Events");
    h1_the2Lab->SetLineColor(1);
    h1_the2Lab->Draw();
    
    c4->cd(2);
    h2_the1LabVSpatheta->SetTitle(text);
    // h2_the1LabVSpatheta->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_the1LabVSpatheta->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_the1LabVSpatheta->GetXaxis()->SetTitleSize(0.05);
    h2_the1LabVSpatheta->GetYaxis()->SetTitleSize(0.05);
    h2_the1LabVSpatheta->GetXaxis()->SetTitle("patheta (deg)");
    h2_the1LabVSpatheta->GetYaxis()->SetTitle("the1Lab (deg)");
    h2_the1LabVSpatheta->SetLineColor(1);
    h2_the1LabVSpatheta->Draw("colz");
    
    c4->cd(3);
    h2_the2LabVSpatheta->SetTitle(text);
    // h2_the2LabVSpatheta->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_the2LabVSpatheta->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_the2LabVSpatheta->GetXaxis()->SetTitleSize(0.05);
    h2_the2LabVSpatheta->GetYaxis()->SetTitleSize(0.05);
    h2_the2LabVSpatheta->GetXaxis()->SetTitle("patheta (deg)");
    h2_the2LabVSpatheta->GetYaxis()->SetTitle("the2Lab (deg)");
    h2_the2LabVSpatheta->SetLineColor(1);
    h2_the2LabVSpatheta->Draw("colz");
    
    c4->cd(4);
    h2_the2LabVSthe1Lab->SetTitle(text);
    // h2_the2LabVSthe1Lab->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_the2LabVSthe1Lab->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_the2LabVSthe1Lab->GetXaxis()->SetTitleSize(0.05);
    h2_the2LabVSthe1Lab->GetYaxis()->SetTitleSize(0.05);
    h2_the2LabVSthe1Lab->GetXaxis()->SetTitle("the1Lab (deg)");
    h2_the2LabVSthe1Lab->GetYaxis()->SetTitle("the2Lab (deg)");
    h2_the2LabVSthe1Lab->SetLineColor(1);
    h2_the2LabVSthe1Lab->Draw("colz");
    
    ymin=0;
    ymax=h1_E1diff->GetMaximum()*1.5;
    c4->cd(5);
    h1_E1diff->SetTitle(text);
    // h1_E1diff->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_E1diff->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_E1diff->GetXaxis()->SetTitleSize(0.05);
    h1_E1diff->GetYaxis()->SetTitleSize(0.05);
    h1_E1diff->GetXaxis()->SetTitle("(E1 - E1meas)/E1");
    h1_E1diff->GetYaxis()->SetTitle("Events");
    h1_E1diff->SetLineColor(1);
    h1_E1diff->Draw();
    
    ymin=0;
    ymax=h1_E2diff->GetMaximum()*1.5;
    c4->cd(6);
    h1_E2diff->SetTitle(text);
    // h1_E2diff->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_E2diff->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_E2diff->GetXaxis()->SetTitleSize(0.05);
    h1_E2diff->GetYaxis()->SetTitleSize(0.05);
    h1_E2diff->GetXaxis()->SetTitle("(E2 - E2meas)/E2");
    h1_E2diff->GetYaxis()->SetTitle("Events");
    h1_E2diff->SetLineColor(1);
    h1_E2diff->Draw();
    
    TCanvas *c5 = new TCanvas("c5", "c5",200,10,1000,1000);
    c5->Divide(2,3);
    
    ymin=0;
    ymax=h1_pi0ratio->GetMaximum()*1.5;
    c5->cd(1);
    h1_pi0ratio->SetTitle(text);
    // h1_pi0ratio->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_pi0ratio->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_pi0ratio->GetXaxis()->SetTitleSize(0.05);
    h1_pi0ratio->GetYaxis()->SetTitleSize(0.05);
    h1_pi0ratio->GetXaxis()->SetTitle("pi0mass meas/gen");
    h1_pi0ratio->GetYaxis()->SetTitle("Events");
    h1_pi0ratio->SetLineColor(1);
    h1_pi0ratio->Draw();
    
    c5->cd(2);
    h2_pi0ratioVSE1->SetTitle(text);
    // h2_pi0ratioVSE1->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_pi0ratioVSE1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1->GetYaxis()->SetTitle("pi0mratio");
    h2_pi0ratioVSE1->SetLineColor(1);
    h2_pi0ratioVSE1->Draw("colz");
    
    
    h2_pi0ratioVSE1->FitSlicesY();
    TH1D *h2_pi0ratioVSE1_0 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_0");
    TH1D *h2_pi0ratioVSE1_1 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_1");
    TH1D *h2_pi0ratioVSE1_2 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_2");
    
    c5->cd(3);
    h2_pi0ratioVSE1_2->Draw();
    
    c5->cd(4);
    xmin = 1;
    xmax=7;
    ymin=0.8;
    ymax=1.2;
    h2_pi0ratioVSE1_1->Fit(Ecorr,"","",xmin,xmax);
    h2_pi0ratioVSE1_1->SetTitle(text);
    // h2_pi0ratioVSE1_1->GetXaxis()->SetRangeUser(xmin,xmax);
    h2_pi0ratioVSE1_1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_1->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_1->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_1->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_1->GetYaxis()->SetTitle("pi0mratio");
    h2_pi0ratioVSE1_1->SetLineColor(1);
    h2_pi0ratioVSE1_1->Draw();
    
    
    c5->cd(5);
    h2_pi0ratioVSE1_1->Draw();
    
    text.Form("p_{#pi}=%.2f, E_{thresh}=%.2f, a=%.3f, b=%.3f, iter=%d",paLab0,Ethresh,BCALa,BCALb,iter);
    TCanvas *c6 = new TCanvas("c6", "c6",200,10,1000,1000);
    c6->Divide(2,3);
    
    ymin=0;
    ymax=h1_pi0ratio_corr1->GetMaximum()*1.5;
    c6->cd(1);
    h1_pi0ratio_corr1->SetTitle(text);
    // h1_pi0ratio_corr1->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_pi0ratio_corr1->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_pi0ratio_corr1->GetXaxis()->SetTitleSize(0.05);
    h1_pi0ratio_corr1->GetYaxis()->SetTitleSize(0.05);
    h1_pi0ratio_corr1->GetXaxis()->SetTitle("pi0mass meas/gen");
    h1_pi0ratio_corr1->GetYaxis()->SetTitle("Events");
    h1_pi0ratio_corr1->SetLineColor(1);
    h1_pi0ratio_corr1->Draw();
    
    c6->cd(2);
    h2_pi0ratioVSE1_corr1->SetTitle(text);
    // h2_pi0ratioVSE1_corr1->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_pi0ratioVSE1_corr1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_corr1->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr1->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr1->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_corr1->GetYaxis()->SetTitle("pi0mratio");
    h2_pi0ratioVSE1_corr1->SetLineColor(1);
    h2_pi0ratioVSE1_corr1->Draw("colz");
    
    
    h2_pi0ratioVSE1_corr1->FitSlicesY();
    TH1D *h2_pi0ratioVSE1_corr1_0 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_corr1_0");
    TH1D *h2_pi0ratioVSE1_corr1_1 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_corr1_1");
    TH1D *h2_pi0ratioVSE1_corr1_2 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_corr1_2");
    
    c6->cd(3);
    h2_pi0ratioVSE1_corr1_2->Draw();
    
    c6->cd(4);
    xmin = 1;
    xmax=7;
    ymin=0.8;
    ymax=1.2;
    h2_pi0ratioVSE1_corr1_1->Fit(Ecorr,"","",xmin,xmax);
    h2_pi0ratioVSE1_corr1_1->SetTitle(text);
    // h2_pi0ratioVSE1_corr1_1->GetXaxis()->SetRangeUser(xmin,xmax);
    h2_pi0ratioVSE1_corr1_1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_corr1_1->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr1_1->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr1_1->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_corr1_1->GetYaxis()->SetTitle("pi0mratio corr(E2)");
    h2_pi0ratioVSE1_corr1_1->SetLineColor(1);
    h2_pi0ratioVSE1_corr1_1->Draw();
    
    
    c6->cd(5);
    h2_pi0ratioVSE1_corr1_0->Draw();
    
    text.Form("p_{#pi}=%.2f, E_{thresh}=%.2f, a=%.3f, b=%.3f, iter=%d check",paLab0,Ethresh,BCALa,BCALb,iter);
    TCanvas *c7 = new TCanvas("c7", "c7",200,10,1000,1000);
    c7->Divide(2,3);
    
    ymin=0;
    ymax=h1_pi0ratio_corr2->GetMaximum()*1.5;
    c7->cd(1);
    h1_pi0ratio_corr2->SetTitle(text);
    // h1_pi0ratio_corr2->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_pi0ratio_corr2->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_pi0ratio_corr2->GetXaxis()->SetTitleSize(0.05);
    h1_pi0ratio_corr2->GetYaxis()->SetTitleSize(0.05);
    h1_pi0ratio_corr2->GetXaxis()->SetTitle("pi0mass meas/gen");
    h1_pi0ratio_corr2->GetYaxis()->SetTitle("Events");
    h1_pi0ratio_corr2->SetLineColor(1);
    h1_pi0ratio_corr2->Draw();
    
    c7->cd(2);
    h2_pi0ratioVSE1_corr2->SetTitle(text);
    // h2_pi0ratioVSE1_corr2->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_pi0ratioVSE1_corr2->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_corr2->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr2->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr2->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_corr2->GetYaxis()->SetTitle("pi0mratio");
    h2_pi0ratioVSE1_corr2->SetLineColor(1);
    h2_pi0ratioVSE1_corr2->Draw("colz");
    
    
    h2_pi0ratioVSE1_corr2->FitSlicesY();
    TH1D *h2_pi0ratioVSE1_corr2_0 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_corr2_0");
    TH1D *h2_pi0ratioVSE1_corr2_1 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_corr2_1");
    TH1D *h2_pi0ratioVSE1_corr2_2 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_corr2_2");
    
    c7->cd(3);
    h2_pi0ratioVSE1_corr2_2->Draw();
    
    c7->cd(4);
    xmin = 1;
    xmax=7;
    ymin=0.8;
    ymax=1.2;
    h2_pi0ratioVSE1_corr2_1->Fit(Ecorr,"","",xmin,xmax);
    h2_pi0ratioVSE1_corr2_1->SetTitle(text);
    // h2_pi0ratioVSE1_corr2_1->GetXaxis()->SetRangeUser(xmin,xmax);
    h2_pi0ratioVSE1_corr2_1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_corr2_1->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr2_1->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_corr2_1->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_corr2_1->GetYaxis()->SetTitle("pi0mratio corr(E1,E2)");
    h2_pi0ratioVSE1_corr2_1->SetLineColor(1);
    h2_pi0ratioVSE1_corr2_1->Draw();
    
    
    c7->cd(5);
    h2_pi0ratioVSE1_corr2_0->Draw();
    
     
    text.Form("p_{#pi}=%.2f, E_{thresh}=%.2f, a=%.3f, b=%.3f, iter=%d %s #DeltaE2",paLab0,Ethresh,BCALa,BCALb,iter,tag.Data());
    TCanvas *c8 = new TCanvas("c8", "c8",200,10,1000,1000);
    c8->Divide(2,3);
    
    ymin=0;
    ymax=h1_pi0ratio_DE2->GetMaximum()*1.5;
    c8->cd(1);
    h1_pi0ratio_DE2->SetTitle(text);
    // h1_pi0ratio_DE2->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_pi0ratio_DE2->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_pi0ratio_DE2->GetXaxis()->SetTitleSize(0.05);
    h1_pi0ratio_DE2->GetYaxis()->SetTitleSize(0.05);
    h1_pi0ratio_DE2->GetXaxis()->SetTitle("pi0mass meas/gen");
    h1_pi0ratio_DE2->GetYaxis()->SetTitle("Events");
    h1_pi0ratio_DE2->SetLineColor(1);
    h1_pi0ratio_DE2->Draw();
    
    c8->cd(2);
    h2_pi0ratioVSE1_DE2->SetTitle(text);
    // h2_pi0ratioVSE1_DE2->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_pi0ratioVSE1_DE2->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_DE2->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_DE2->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_DE2->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_DE2->GetYaxis()->SetTitle("pi0mratio");
    h2_pi0ratioVSE1_DE2->SetLineColor(1);
    h2_pi0ratioVSE1_DE2->Draw("colz");
    
    
    h2_pi0ratioVSE1_DE2->FitSlicesY();
    TH1D *h2_pi0ratioVSE1_DE2_0 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_DE2_0");
    TH1D *h2_pi0ratioVSE1_DE2_1 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_DE2_1");
    TH1D *h2_pi0ratioVSE1_DE2_2 = (TH1D*) gROOT->FindObject("h2_pi0ratioVSE1_DE2_2");
    
    c8->cd(3);
    h2_pi0ratioVSE1_DE2_2->Draw();
    
    c8->cd(4);
    xmin = 1.5;
    xmax=7;
    ymin=0.8;
    ymax=1.2;
    h2_pi0ratioVSE1_DE2_1->Fit(Ecorr,"","",xmin,xmax);
    h2_pi0ratioVSE1_DE2_1->SetTitle(text);
    // h2_pi0ratioVSE1_DE2_1->GetXaxis()->SetRangeUser(xmin,xmax);
    h2_pi0ratioVSE1_DE2_1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_pi0ratioVSE1_DE2_1->GetXaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_DE2_1->GetYaxis()->SetTitleSize(0.05);
    h2_pi0ratioVSE1_DE2_1->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_pi0ratioVSE1_DE2_1->GetYaxis()->SetTitle("pi0mratio #DeltaE2");
    h2_pi0ratioVSE1_DE2_1->SetLineColor(1);
    h2_pi0ratioVSE1_DE2_1->Draw();
    
    
    c8->cd(5);
    h2_pi0ratioVSE1_DE2_0->Draw();
    
     TCanvas *c9 = new TCanvas("c9", "c9",200,10,1000,1000);
     c9->Divide(2,2);
    
    c9->cd(1);
    h2_pi0ratioVSE1_1->Draw();
    
    c9->cd(2);
    h2_pi0ratioVSE1_corr1_1->Draw();
    
    c9->cd(3);
    h2_pi0ratioVSE1_corr2_1->Draw();
    
    c9->cd(4);
    h2_pi0ratioVSE1_DE2_1->Draw();
    
    
     TCanvas *c10 = new TCanvas("c10", "c10",200,10,1000,1000);
     c10->Divide(2,2);
    
    gPad->SetGridx();
    gPad->SetGridy();
    
    c10->cd(1);
    gPad->SetGridx();
    gPad->SetGridy();
    h2_pi0ratioVSE1_1->Draw();
    
    c10->cd(2);
    gPad->SetGridx();
    gPad->SetGridy();
    h2_pi0ratioVSE1_corr1_1->Draw();
    
    c10->cd(3);
    gPad->SetGridx();
    gPad->SetGridy();
    h2_pi0ratioVSE1_corr2_1->Draw();
    
    
    xmin=1.5;
    xmax=6;
    c10->cd(4);
    gPad->SetGridx();
    gPad->SetGridy();
    TH2D* h2_pi0ratioVSE1_DE2_1_copy = (TH2D*)h2_pi0ratioVSE1_DE2_1->Clone("h2_pi0ratioVSE1_DE2_1_copy");
    // inverse->FixParameter(1,(DeltaE2_min+DeltaE2_max)/2.);    // fix energy
    // h2_pi0ratioVSE1_DE2_1_copy->Fit(inverse,"","",xmin,xmax);
    inverse2->FixParameter(2,(DeltaE2_min+DeltaE2_max)/2.);    // fix energy
    h2_pi0ratioVSE1_DE2_1_copy->Fit(inverse2,"","",xmin,xmax);
    
    // output nonlinear parameters
    
    TF1 *fitcorr = h2_pi0ratioVSE1_DE2_1->GetFunction("Ecorr");   // make sure parameters are from plot with E2_corr ONLY
    q0 = fitcorr->GetParameter(0);
    q1 = fitcorr->GetParameter(1);
    q2 = fitcorr->GetParameter(2);
    q3 = fitcorr->GetParameter(3);
    q4 = fitcorr->GetParameter(4);
    q5 = fitcorr->GetParameter(5);
    q6 = fitcorr->GetParameter(6);
    q7 = fitcorr->GetParameter(7);
    q8 = fitcorr->GetParameter(8);
    
    outlist << q0 << " " << q1 << " " << q2 << " " << q3 << " " <<  q4 << " " << q5 << " " << q6 << " " << q7 << " " << q8 << endl;

    outlist.close();
    
    TString outfile;
    outfile.Form(filename+"_iter%d_%.0f_%.0f_%.0f_%.0f",iter,1000*paLab0,1000*Ethresh,1000*BCALa,1000*BCALb);
    // c1->SaveAs(outfile+".pdf(");
    c2->SaveAs(outfile+".pdf(");
    c3->SaveAs(outfile+".pdf");
    c4->SaveAs(outfile+".pdf");
    c5->SaveAs(outfile+".pdf");
    c6->SaveAs(outfile+".pdf");
    c7->SaveAs(outfile+".pdf");
    c8->SaveAs(outfile+".pdf");
    c9->SaveAs(outfile+".pdf");
    c10->SaveAs(outfile+".pdf)");
    
}