// set colors/fills via parameter

static Double_t step[2] =  { 0.01, 0.001};
static Double_t vstart[2] = { 0.01, 0.01};

static Double_t G_E_S = 0;
static Double_t G_M_S = 0;
static Double_t eG_E_S = 0;
static Double_t eG_M_S = 0;
static Double_t MinLikli = 0;

static const Double_t PI=3.1415927;
static const Int_t NMEAS=8;
// new tilted SAMPLE treated like others
//                           
static Double_t slp[NMEAS]; 
static Double_t intr[NMEAS];
static Double_t eintr[NMEAS];
static Double_t emeas[NMEAS];
static Int_t dset[NMEAS];
static Int_t npts;
static Int_t dind[NMEAS]={-1};
static Int_t dcol[NMEAS];
static Int_t dfill[NMEAS];

static TMinuit *minuit;


//#define PRELIM

Int_t LoadOne(Int_t choice) {
  // choice 1: SAMPLE
  //        2: A4
  //        3: HAPPEX-H '04
  //        4: HAPPEX-He '04
  //        5: HAPPEX-H '05
  //        6: HAPPEX-He '05
  if (choice==1) { // SAMPLE
    dind[1] = npts;
    dset[npts] = 1;
    slp[npts] =  -1.668;   // from betsy
    //intr[npts] =  0.590; // From betsy
    //emeas[npts] = 0.280; // From betsy
    emeas[npts] = 0.335; // from review, with slope added
    intr[npts] =  0.617; // From review, with slope added
//     slp[npts] =  -1.642; // from doug
//     intr[npts] =  0.651; // From doug
//     emeas[npts] = 0.318; // From doug, eintr = .612
    //dcol[npts] = 8;
    dcol[npts] = 49;
    dfill[npts] = 3005;
  } if (choice==2) { // A4
    dind[2] = npts;
    dset[npts] = 2;
    slp[npts] =  -0.106;    
    intr[npts] =  0.071;
    emeas[npts] = 0.036;
    //dcol[npts] = 49;
    dcol[npts] = 8;
    dfill[npts] = 3006;
  } if (choice==3) { // HAPPEX-H '04
    dind[3] = npts;
    dset[npts] = 3;
    slp[npts] = -0.0802;    
    intr[npts] = 0.030;
    emeas[npts] = 0.028;
    dcol[npts] = 2;
        dfill[npts] = 3007;
#ifdef PRELIM
    slp[npts] = -0.080;    
    intr[npts] = 0.0339;
    emeas[npts] = 0.030;
#endif
  } if (choice==4) {  // HAPPEX-He '04
    dind[4] = npts;
    dset[npts] = 4;
    slp[npts] =  0.00;    
    intr[npts] = -0.039;
    emeas[npts] = 0.043;
    dcol[npts] = 9;
    dfill[npts] = 3004;
#ifdef PRELIM
    slp[npts] =  0.00;    
    intr[npts] = -0.019;
    emeas[npts] = 0.0485;
#endif
  } if (choice==5) {  // HAPPEX-H '05
    dind[5] = npts;
    dset[npts] = 5;
    slp[npts] = -0.080;    
    intr[npts] = 0.0;
    emeas[npts] = 0.010;
    dcol[npts] = 9;
    dfill[npts] = 3001;
  } if (choice==6) {  // HAPPEX-He '05
    // HAPPEX-He
    dind[6] = npts;
    dset[npts] = 6;
    slp[npts] =  0.00;    
    intr[npts] = 0.0;
    emeas[npts] = 0.015;
    dcol[npts] = 2;
    dfill[npts] = 3002;
    // 3% measurement of 7.8 ppm asym, 
    // with comparable systematic error = 0.015 on GES
  } if (choice==7) { // G0
    dind[7] = npts;
    dset[npts] = 7;
    slp[npts] = -0.080;    
    intr[npts] = 0.0485;
    emeas[npts] = 0.031;
    dcol[npts] = 7;
    dfill[npts] = 3014;
  } 
  return ++npts;
}

Int_t LoadData(Int_t choice) {
  // choice 1: all Ges/Gms at Q2=0.1Gev
  //        2: just HAPPEX '04
  //        3: projected HAPPEX '05
  //        4: just sample and A4
  //        5: project HAPPEX '05, with values matching current best fit
  //        6: all Ges/Gms at Q2=0.1Gev with contour
  Int_t ind;
  if (choice==1 || choice==6) {
    npts =0;
    npts = LoadOne(1); // SAMPLE
    npts = LoadOne(4); // HAPPEX-He '04
    npts = LoadOne(2); // A4
    npts = LoadOne(3); // HAPPEX-H '04


    for (Int_t ipt=0 ; ipt<npts; ipt++) {
      eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
    }
    return 1;

  } else  if (choice==2) {
    npts =0;
    npts = LoadOne(3); // HAPPEX-H '04
    npts = LoadOne(4); // HAPPEX-He '04
    for (Int_t ipt=0 ; ipt<npts; ipt++) {
      eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
    }
    return 1;
  } else if (choice==3) {
    // 3: projected HAPPEX '05
    npts =0;
    npts = LoadOne(5); // HAPPEX-H '05
    npts = LoadOne(6); // HAPPEX-He '05
    for (Int_t ipt=0 ; ipt<npts; ipt++) {
      eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
    }
    return 1;
  } else if (choice==4) {
    //        4: just sample and A4
    // A4
    npts =0;
    npts = LoadOne(1); // SAMPLE
    npts = LoadOne(2); // A4
    for (Int_t ipt=0 ; ipt<npts; ipt++) {
      eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
    }
    return 1;
  } else if (choice==5) {
    // 5: project HAPPEX '05, with values matching current best fit
    npts =0;
    npts = LoadOne(3); // HAPPEX-H '04
    npts = LoadOne(4); // HAPPEX-He '04
    npts = LoadOne(1); // SAMPLE
    npts = LoadOne(2); // A4
    for (Int_t ipt=0 ; ipt<npts; ipt++) {
      eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
    }
    return 1;
  } else if (choice==7) {
    npts =0;
    npts = LoadOne(1); // SAMPLE
    npts = LoadOne(4); // HAPPEX-He '04
    npts = LoadOne(2); // A4
    npts = LoadOne(7); // G0
    npts = LoadOne(3); // HAPPEX-H '04

    for (Int_t ipt=0 ; ipt<npts; ipt++) {
      eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
    }
    return 1;

  }

  return 0;  
}


void FF_ALL(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, int iflag){

  Double_t total = 0;
  
  Double_t gm = x[0];
  Double_t ge = x[1];

  if (iflag==1) {
    // do nothing special
  }

  //   Double_t slp[npt] =  {-0.080,-0.0,  -0.106}; 
  //   Double_t intr[npt] = {0.00,  0.0, 0.074};
  //   Double_t eintr[npt] ={0.015, 0.015,0.036};
  
  for (Int_t ip =0; ip<npts; ip++) {
//     Double_t dsq = (ge-intr[ip])**2 
//       + (slp[ip]**2)*gm**2 - 2*slp[ip]*gm*(ge-intr[ip]);
//     dsq /= (1+slp[ip]**2);
//     dsq /= eintr[ip]**2;
    Double_t x0 = (gm + slp[ip]*(ge - intr[ip]))/(1+slp[ip]*slp[ip]);
    Double_t dsq = (gm-x0)*(gm-x0) + 
      (ge - intr[ip] - slp[ip]*x0)*(ge - intr[ip] - slp[ip]*x0);
    dsq /= emeas[ip]**2;
    total+=dsq;
  }
    
  if (iflag ==3) {
    // do nothing, allow numeric differentiation
  }
    
  f = total;
  
}

Int_t minFF_Likli(){


  Double_t minos_gms_p,minos_gms_m,minos_gms,gcc_gms;
  Double_t minos_ges_p,minos_ges_m,minos_ges,gcc_ges;
  //  gROOT->LoadMacro("FF_ALL.C");
  //  TMinuit *minuit = new TMinuit(2);
  minuit = new TMinuit(2);
  minuit->SetFCN(FF_ALL);
  
  Double_t arglist[10];
  Int_t ierflg = 0;
  
  arglist[0] = 1;
  minuit->mnexcm("SET ERR", arglist ,1,ierflg);

  // Set starting values and step sizes for parameters
  minuit->mnparm(0, "GMs", vstart[0], step[0], 0,0,ierflg);
  minuit->mnparm(1, "GEs", vstart[1], step[1], 0,0,ierflg);

// Now ready for minimization step
  arglist[0] = 500;
  arglist[1] = 1.;

  Double_t df[2] = {0};
  Double_t Chisq;
  minuit->Eval(2,df,Chisq,vstart,1);

  cout << Chisq << endl;

  FF_ALL(2,df,Chisq,vstart,2);
  cout << endl;
  cout << Chisq << endl;
  cout << endl;

  minuit->mnexcm("MIGRAD", arglist ,2,ierflg);

  cout << "****  Calling Minos ****" << endl;
  arglist[1]=0;
  minuit->mnexcm("MINOS", arglist ,2,ierflg);

// Print results
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);

   minuit->GetParameter(0,G_M_S,eG_M_S);
   minuit->GetParameter(1,G_E_S,eG_E_S);
   MinLikli = amin;

   //   minuit->mnprin(4,amin);

   cout << "*********************************" << endl;
   cout << Form(" Gms = %6.3f +/- %5.3f " ,G_M_S,eG_M_S) << endl;
   cout << Form(" Ges = %6.3f +/- %5.3f " ,G_E_S,eG_E_S) << endl;
   cout << "  Minimum Liklihood = " << MinLikli << endl;
   cout << "*********************************" << endl;

}

Double_t FFLikelihood(Double_t * Gs, Double_t *pars){

  Double_t gm = Gs[0];
  Double_t ge = Gs[1];


  // calc chisquared
  // new tilted SAMPLE treated like others
  //  Double_t Chisq_Min = 0.483483;
  
  Double_t Likli =0;
  for (Int_t ip =0; ip<npts; ip++) {
//     Double_t dsq = (ge-intr[ip])**2 
//       + (slp[ip]**2)*gm**2 - 2*slp[ip]*gm*(ge-intr[ip]);
//     dsq /= (1+slp[ip]**2);
//     dsq /= eintr[ip]**2;
    Double_t x0 = (gm + slp[ip]*(ge - intr[ip]))/(1+slp[ip]*slp[ip]);
    Double_t dsq = (gm-x0)*(gm-x0) + 
      (ge - intr[ip] - slp[ip]*x0)*(ge - intr[ip] - slp[ip]*x0);
    dsq /= emeas[ip]**2;
    Likli+=dsq;
  }

  Likli -= MinLikli;

  return Likli;
}


Int_t Impact(int choice=1, int print=0) {
  // choice 1: all Ges/Gms at Q2=0.1Gev
  //        2: just HAPPEX '04
  //        3: projected HAPPEX '05
  //        4: just sample and A4
  //        5: project HAPPEX '05, with values matching current best fit
  //        6: all Ges/Gms at Q2=0.1Gev with contour
  // if choice involve contours, prepare liklihood function by finding minimum
  // if choice ==5, first modify central value for 2005 results

  //  Bool_t useLabels = kFALSE;
  Bool_t useLabels = kTRUE;
  Bool_t useTheory = kFALSE;
  Bool_t drawZeroPt = kTRUE;
  LoadData(choice);
  for (Int_t ipt=0; ipt<npts; ipt++) {
    cout << slp[ipt] << "  " << intr[ipt] << "   " << eintr[ipt] << endl;
  }

  Double_t gesrange=.185, gmsrange=1.501;


  TF2* ftry;
  Int_t n_cont = 0;
  TGraph *gEcont[10];
  if (choice==5 || choice==6 || choice==7) {
    minFF_Likli();  
    if (choice==5) {
      // overwrite HAPPEX central value with new center, and proposed precision
      intr[dind[3]] = G_E_S - slp[0]*G_M_S;
      emeas[dind[3]]= 0.010;
      eintr[dind[3]] = emeas[dind[3]] / cos( atan( slp[dind[3]] ) );
      cout << Form(" HAPPEX-H result: ges + 0.08*gms = %6.3f",intr[dind[3]]) 
	   << endl;
      intr[dind[4]] = G_E_S;    Double_t gzed[2] = {0.0,0.0};
    TString zeropt = Form("Likelihood of zero, zero = %7.3f", FFLikelihood(gzed,gzed));
    cout << zeropt << endl;

      emeas[dind[4]]= 0.015;
      eintr[dind[4]] = emeas[dind[4]] / cos( atan( slp[dind[4]] ) );
      minFF_Likli();  
    }

    ftry = new TF2("ftry",FFLikelihood,
		   -1*gmsrange,gmsrange,-1*gesrange,gesrange,0);

    //1.0   // "1-sigma"
    //2.30  // 68.27%
    //4.0   // "2-sigma"
    //4.61  // 90%
    //5.99  // 95%
    //9.21  // 99%
    Int_t nmark = 500;
    minuit->SetErrorDef(5.99);
    //    minuit->SetErrorDef(1.0);
    minuit->SetPrintLevel(0);
    gEcont[n_cont] = (TGraph*)minuit->Contour(nmark,0,1);
    gEcont[n_cont]->SetLineWidth(3);
    gEcont[n_cont]->SetLineStyle(1);
    gEcont[n_cont]->SetLineColor(1);
    n_cont++;



//     minuit->SetErrorDef(4.0);
//     minuit->SetPrintLevel(0);
//     gEcont[n_cont] = (TGraph*)minuit->Contour(nmark,0,1);
//     gEcont[n_cont]->SetLineWidth(2);
//     gEcont[n_cont]->SetLineStyle(2);
//     gEcont[n_cont]->SetLineColor(1);
//     n_cont++;

//     npts =LoadOne(7);
//     for (Int_t ipt=0 ; ipt<npts; ipt++) {
//       eintr[ipt] = emeas[ipt] / cos( atan( slp[ipt] ) );
//     }
//     minuit->SetErrorDef(5.99);
//     gEcont[n_cont] = (TGraph*)minuit->Contour(nmark,0,1);
//     gEcont[n_cont]->SetLineWidth(3);
//     gEcont[n_cont]->SetLineStyle(1);
//     gEcont[n_cont]->SetLineColor(2);
//     n_cont++;

  }
  
  gStyle->SetCanvasColor(kWhite);
  gStyle->SetCanvasBorderSize(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderSize(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadLeftMargin(0.12);
  gStyle->SetPadBottomMargin(0.12);

  TCanvas *canv = new TCanvas("canv","Ges/Gms plot",600,600);
  canv->SetFillStyle(4000);
  canv->SetFillColor(0);
  //  canv->SetGrid();
  TH1F* cFrame = canv->DrawFrame(-gmsrange,-gesrange,gmsrange,gesrange);
  canv->SetTicks(1);

  if (useLabels) {
    TLatex *tex = new TLatex(-gmsrange-0.4,-0.01,"G^{s}_{E}");
    tex->SetTextSize(0.07);
    tex->SetTextFont(2);
    tex->Draw();
    
    tex = new TLatex(-0.1,-gesrange-0.04,"G^{s}_{M}");
    tex->SetTextSize(0.07);
    tex->SetTextFont(2);
    tex->Draw();
  }

  //   // Ges Gms paramaters;
  Double_t slope, intercept, inter_err;
  const Int_t n=50;

  for (Int_t is = 0; is<npts; is++) {
    slope = slp[is];
    intercept = intr[is];
    inter_err = eintr[is];

    Double_t onegms[n] = {0};
    Double_t oneges[n] = {0};
    Double_t onegesmin[n] = {0};
    Double_t onegesmax[n] = {0};
    Double_t onegesmin_p[n] = {0};
    Double_t onegesmax_p[n] = {0};
    for (Int_t i=0;i<n;i++) {
      onegms[i] = -gmsrange + i*(2*gmsrange/(n-1));
      oneges[i] = intercept + slope*onegms[i];
      onegesmin[i] = oneges[i] - inter_err;
      onegesmax[i] = oneges[i] + inter_err;
    }
    
    TGraph *onegrmin   = new TGraph(n,onegms,onegesmin);
    TGraph *onegrmax   = new TGraph(n,onegms,onegesmax);
    TGraph *onegrmin_p = new TGraph(n,onegms,onegesmin_p);
    TGraph *onegrmax_p = new TGraph(n,onegms,onegesmax_p);
    TGraph *onegr      = new TGraph(n,onegms,oneges);
    TGraph *onegrshade = new TGraph(2*n);
    TGraph *onegrshade_p = new TGraph(2*n);
    
    for(Int_t i=0;i<n;i++) {
      onegrshade->SetPoint(i,onegms[i],onegesmax[i]);
      onegrshade->SetPoint(n+i,onegms[n-i-1],onegesmin[n-i-1]);
      onegrshade_p->SetPoint(i,onegms[i],onegesmax_p[i]);
      onegrshade_p->SetPoint(n+i,onegms[n-i-1],onegesmin_p[n-i-1]);
    }
    int onecolor = dcol[is];
    int onefill = dfill[is];
    
    onegrshade->SetFillStyle(onefill);
    onegrshade->SetFillColor(onecolor);
    onegrshade->Draw("F");
    onegrmin->SetLineColor(onecolor+100);
    onegrmax->SetLineColor(onecolor+100);
//     onegrmin->SetLineColor(onecolor);
//     onegrmax->SetLineColor(onecolor);
    onegrmax->SetLineWidth(2);
    onegrmin->SetLineWidth(2);
    onegrmin->Draw("L");
    onegrmax->Draw("L");

    onegr->SetLineStyle(2);
    onegr->SetLineWidth(4);
    onegr->SetLineColor(onecolor+100);
    //    onegr->SetLineColor(onecolor);
    onegr->SetMarkerColor(onecolor);
    onegr->SetMarkerStyle(21);
    onegr->Draw("C");


    if (useLabels) {
      switch (dset[is]): {

      case 1:
	TLatex *samptex = new TLatex(0.51,0.155,"SAMPLE with");
	TLatex *samptex2 = new TLatex(0.51,0.135,"G_{A} calculation");
	samptex->SetTextSize(0.035);
	samptex->SetTextColor(onecolor+100);
	//	samptex->SetTextColor(onecolor);
	samptex->Draw();
	samptex2->SetTextSize(0.04);
	samptex2->SetTextColor(onecolor+100);
	//	samptex2->SetTextColor(onecolor);
	samptex2->Draw();
	break; 

      case 2:
	TLatex *a4tex  = new TLatex(-0.40,0.155,"A4");
	a4tex->SetTextSize(0.04);
	a4tex->SetTextColor(onecolor+100);
	//	a4tex->SetTextColor(onecolor);
	a4tex->SetTextAngle(-40);
	a4tex->Draw();
	break; 

      case 3: 
	Double_t hytexy = slope*(-1.4) + intercept - inter_err*1.0 - 0.018;
	TLatex *hytex = 
	  new TLatex(-1.4,hytexy,"HAPPEX-H ('04)");
#ifdef PRELIM
	hytex = new TLatex(-1.3,hytexy,"HAPPEX-H (preliminary)");	
#endif
	hytex->SetTextSize(0.04);
	hytex->SetTextColor(onecolor+100);
	//	hytex->SetTextColor(onecolor);
	hytex->SetTextAngle(-32.5);
	hytex->Draw();    
	break;

      case 4: 
	TLatex *hetex = new TLatex(-1.4,0.01,"HAPPEX-^{4}He ('04)");
#ifdef PRELIM
	hetex = new TLatex(-0.8,-0.045,
			   "#splitline{HAPPEX-^{4}He}{(preliminary)}");
#endif
	hetex->SetTextSize(0.04);
	//	hetex->SetTextColor(onecolor);
	hetex->SetTextColor(onecolor+100);
	hetex->Draw();    
	break;

      case 5: 
	Double_t hytexy = slope*(-1.3) + intercept - inter_err*1.0 - 0.015;
	TLatex *hytex = new TLatex(-1.4,hytexy,"HAPPEX-H (projected)");
	hytex->SetTextSize(0.04);
	hytex->SetTextColor(onecolor+100);
	//	hytex->SetTextColor(onecolor);
	hytex->SetTextAngle(-30);
	hytex->Draw();    
	break;

      case 6: 
	TLatex *hetex = new TLatex(-0.8,-0.045,
			   "#splitline{HAPPEX-^{4}He}{(projected)}");
	hetex->SetTextSize(0.04);
	hetex->SetTextColor(onecolor+100);
	//	hetex->SetTextColor(onecolor);
	hetex->Draw();    
	break;

      case 7: 
	Double_t g0texy = slope*(-1.3) + intercept + inter_err*1.0 - 0.02;
	TLatex *g0tex = 
	  new TLatex(-1.25,g0texy,"G0 (extrapolated)");
	g0tex->SetTextSize(0.04);
	g0tex->SetTextColor(onecolor+100);
	//	g0tex->SetTextColor(onecolor);
	g0tex->SetTextAngle(-33);
	g0tex->Draw();    
	break;

      }
    }

  }

  for (Int_t ic=0; ic<n_cont; ic++) {
    gEcont[ic]->Draw("l");
  }
  TLatex *ctxt95 = new TLatex(0.85,0.03,"95% C.L.");
  ctxt95->SetTextSize(0.04);
  ctxt95->SetTextFont(22);
  ctxt95->SetTextColor(1);
  ctxt95->Draw();
  
  if (ftry) {
    TMarker* bestPt = new TMarker(G_M_S,G_E_S,20);
    bestPt->Draw();
  }

  if (drawZeroPt) {
    TMarker* zeroPt = new TMarker(0.0,0.0,20);
    zeroPt->SetMarkerSize(1.2);
    zeroPt->Draw();

    Double_t gzed[2] = {0.0,0.0};
    TString strzeropt = Form("Likelihood of zero, zero = %7.3f", FFLikelihood(gzed,gzed));
    cout << strzeropt << endl;
  }


  //  if (useTheory) {

    // Draw theory points on impact plot
    
    const Int_t n=11;
    Double_t mean = 0.-0.046;
    Double_t mean_err = 0.019;
    
    Double_t lat_gms[n], lat_ges[n];
    Double_t lat_gmsmin[n], lat_gmsmax[n];
    for (Int_t i=0;i<n;i++) {
      lat_gms[i] = mean;
      lat_ges[i] = -gesrange + i*(2*gesrange/(n-1));
      lat_gmsmin[i] = lat_gms[i] - mean_err;
      lat_gmsmax[i] = lat_gms[i] + mean_err;
    }
    
    // Draw Lattice QCD band, Leinweber et al.
    TGraph *lat_grmin   = new TGraph(n,lat_gmsmin,lat_ges);
    TGraph *lat_grmax   = new TGraph(n,lat_gmsmax,lat_ges);
    TGraph *lat_gr      = new TGraph(n,lat_gms,lat_ges);
    TGraph *lat_grshade = new TGraph(2*n);
    
    for(Int_t i=0;i<n;i++) {
      lat_grshade->SetPoint(i,lat_gmsmax[i],lat_ges[i]);
      lat_grshade->SetPoint(n+i,lat_gmsmin[n-i-1],lat_ges[n-i-1]);
    }
    lat_grshade->SetFillStyle(0);
    lat_grshade->SetFillColor(6+100);
    //    lat_grshade->SetFillColor(6);
    lat_grshade->SetLineColor(6);
    lat_grmin->SetLineColor(6);
    lat_grmax->SetLineColor(6);
    lat_grshade->Draw("F");
     lat_grmin->Draw("L");
     lat_grmax->Draw("L");
      
    lat_gr->SetLineStyle(2);
    lat_gr->SetLineColor(6);
    lat_gr->SetLineWidth(2);
    lat_gr->Draw("C");
    
    TLatex *lein = new TLatex(-0.7,-0.115,"Leinweber ");
    lein->SetTextSize(0.04);
    lein->SetTextFont(22);
    lein->SetTextColor(6);
    lein->Draw();
    TLatex *lein2 = new TLatex(-0.7,-0.13,"et al. (2005)");
    lein2->SetTextSize(0.04);
    lein2->SetTextFont(32);
    lein2->SetTextColor(6);
    lein2->Draw();


  if (useTheory) {
    const Int_t ntheory = 6;
    // // ntheory = 0 => Weigel, Nambu-Jona-Lasinio model w/ soliton picture
    // ntheory = 0 => Silva, SU(3) chiral quark soliton model, kaon values
    // ntheory = 1 => Lyubovitskij, perturbative chiral quark model at one loop
    //    // ntheory = 3 => Hemmert, heavy baryon chiral perturbation theory
    // ntheory = 2 => Lewis, quenched lattice QCD w/ quenched ChiPT
    // ntheory = 3 => Hammer, dispersion relations and KK-bar continuum
    // ntheory = 4 => Park, Skyrme model
    // ntheory = 5 => Drechsel, dispersion-theoretical fit

    Double_t gms[ntheory] = { 0.078, -0.036,  0.050, -0.500, -0.036, -0.310 }; //0.155,0.235
    Double_t ges[ntheory] = { 0.045,  0.004,  0.015, -0.172, -0.061, -0.084 }; //0.039,-0.021,
    Double_t gmserr[ntheory] = {  0.012,      0,  0.060,  0.170,  0,      0.03};
    Double_t geserr[ntheory] = {  0.017,      0,  0.005,  0.007,  0,      0.013};

//     Double_t gmserr[ntheory] = { 0, 0.012,      0,  0.452, 0.060,  0.170 };
//     Double_t geserr[ntheory] = { 0, 0.017,      0,  0.038, 0.005,  0.007 };
   
//     TGraph *g_weigel = new TGraph(1,&(gms[0]),&(ges[0]));
//     g_weigel->SetMarkerSize(1.5);
//     g_weigel->SetMarkerStyle(20);
//     TLatex *weigel = new TLatex(0.2,0.030,"[14]");
//     weigel->SetTextSize(0.04);

//     TGraph *g_silva = new TGraph(1,&(gms[1]),&(ges[1]));
//     g_silva->SetMarkerSize(1.5);
//     g_silva->SetMarkerStyle(21);
    TLatex *silva = new TLatex(0.12,0.035,"[19]");
    silva->SetTextSize(0.04);

//     TGraph *g_lyub = new TGraph(1,&(gms[2]),&(ges[2]));
//     g_lyub->SetMarkerSize(1.8);
//     g_lyub->SetMarkerStyle(22);
    TLatex *lyub = new TLatex(-0.3,0.0,"[20]");
    lyub->SetTextSize(0.04);

//     TLatex *hemmert = new TLatex(0.35,-0.035,"[17]");
//     hemmert->SetTextSize(0.04);
    TLatex *lewis = new TLatex(0.15,0.01,"[21]");
    lewis->SetTextSize(0.04);
    TLatex *lein = new TLatex(-0.15,-0.14,"[22]");
    lein->SetTextSize(0.04);
    TLatex *hammer = new TLatex(-0.75,-0.165,"[18]");
    hammer->SetTextSize(0.04);
    TLatex *park = new TLatex(-0.30,-0.065,"[16]");
    park->SetTextSize(0.04);
    TLatex *drechsel = new TLatex(-0.55,-0.10,"[17]");
    drechsel->SetTextSize(0.04);

    TGraphErrors *theory = new TGraphErrors(ntheory,gms,ges,gmserr,geserr);
    theory->SetMarkerStyle(21);
    theory->SetMarkerSize(1);
    theory->SetMarkerColor(1);
    theory->SetLineWidth(2);
    theory->Draw("P");  
    
    //    g_weigel->Draw("p");
    //    weigel->Draw();
    //    g_silva->Draw("p");
    silva->Draw();
    //    g_lyub->Draw("p");
    lyub->Draw();
    //    hemmert->Draw();
    lewis->Draw();
    lein->Draw();
    hammer->Draw();
    park->Draw();
    drechsel->Draw();
    

    //    TLegend *leg = new TLegend(0.15,0.15,0.4,0.4);
    //    leg->Add(weigel,20,"p");
    //    leg->Draw();
  }
 
  if (useLabels) {
    //    TLatex *qtex = new TLatex(-1.25,-0.14,"Q^{2} = 0.1 GeV^{2}");
    TLatex *qtex = new TLatex(-1.34,-0.164,"Q^{2} = 0.1 GeV^{2}");
    qtex->SetTextFont(2);
    qtex->SetTextSize(0.045);
    qtex->Draw();
  }

  if(useLabels && (choice==2 || choice==3 || choice==5)) {
    TLatex *titletex;
    TString sTitle;

    switch (choice): {
    case 2: 
      sTitle = "#splitline{HAPPEX 2004}{}";
      titletex = new TLatex(-0.1,0.09,sTitle);
      titletex->SetTextSize(0.06);
      break;

    case 3:
      sTitle = "#splitline{After 2005 run:}{#splitline{Anticipated results}{(centered at zero)}}";
      titletex = new TLatex(-0.3,0.09,sTitle);
      titletex->SetTextSize(0.06);
      break;

    case 5:
      sTitle = "#splitline{After 2005 run:}{#splitline{Anticipated results}{(centered at current best fit)}}";
      titletex = new TLatex(-0.23,-0.09,sTitle);
      titletex->SetTextSize(0.04);
      break;

    }
    titletex->Draw();
  }



  canv->Update();
  //  return;
  TString canvname;
  switch (choice): {
  case 1:
    canvname = "2004_current_Ges_Gms";
    break;

  case 2:
    canvname = "2004_HAPPEX_ONLY_Ges_Gms";
    break;

  case 3:
    canvname = "2005_HAPPEX_ONLY_Ges_Gms";
    break;

  case 4:
    canvname = "No_HAPPEX_Ges_Gms";
    break;

  case 5:
    canvname = "2005_special_contour";
    break;

  case 6:
    canvname = "2004_current_contour";
    break;

  case 7:
    canvname = "2004_ALL";
    break;
 
  }

  canv->cd(0);

  if (!useLabels) canvname += "_nl"; 
  if (print==1) {
    canv->Print(""+canvname+".eps");
    canv->Print(""+canvname+".gif");
  }

}