const int maxF = 4;
TH1F* bkgdFrac[maxF]; //f's for background correction [tau,Z,EEMC,QCD]
TH1F* beta; //1 - sum(f's) -> polarization dilution factor
TH1F* betaSyst;
TH1F* betaSystUp;
TH1F* betaSystDn;
TH1F* zFracSecondEEMC;

ofstream outSummary;
TFile *outF;

TCanvas *cc;
TCanvas *dd;

const int maxEta = 8;
//storage of Second EEMC background histos (data and Z MC) for mirror etaBin
TH1F* mirrorSecondEEMC[2][maxEta];
TH1F* mirrorZSecondEEMC[2][maxEta];
const int mirrorEtaBin[maxEta] = {4,3,2,1,0,0,0,8};

void plBkgdEtaDep(int year=2012){
  
  gStyle->SetOptStat(0);
  gStyle->SetOptDate(0);

  outSummary.open(Form("bkgdSummary_%d.txt",year));

  TString outFile = "betaBkgdModel_"; outFile+=year; outFile+=".hist.root";
  TString name[maxF]={"tau","zee","end","qcd"};

  for(int charge=-1; charge<=1; charge+=2) {
    char Q = 'N'; if(charge>0) Q = 'P';
    //if(charge > 0) continue;

    //create output histo file
    if(charge==-1) outF = new TFile(outFile,"RECREATE"); //create new file
    else if(charge==1) outF = new TFile(outFile,"UPDATE"); //write P to same file as N

    //write summary text file
    outSummary<<"W"<<Q<<" background summary: yields 25<ET<50 GeV"<<endl;
    outSummary<<"starPhysEtaBin, rawYield, qcdBkgd, secondEEMC, zeeBkgd, wTauYield, totalBkgd, wYield, beta +/- err, +syst, -syst"<<endl;

    //get Second EEMC background histos (data and Z MC) for mirror etaBin
    for(int ieta=0; ieta<maxEta; ieta++){
      if(ieta>=4 && ieta<=6) continue;
      plBkgdEtaDepX(charge,Form("Eta%d",ieta+1),ieta+1,true,year);
    }

    //initialize background fraction histograms
    for(int iF=0; iF<maxF; iF++){
      TString title=Form("F_%s%c_%d",name[iF].Data(),Q,year);
      bkgdFrac[iF] = new TH1F(title,title,maxEta,0.5,maxEta+0.5);
    }
    TString titleBeta=Form("beta%c_%d",Q,year);
    beta = new TH1F(titleBeta,Form("%s; starPhysEta bin; #beta",titleBeta.Data()),maxEta,0.5,maxEta+0.5);
    TString titleBetaSyst=Form("betaSyst%c_%d",Q,year);
    betaSyst = new TH1F(titleBetaSyst,titleBetaSyst,maxEta,0.5,maxEta+0.5);
    TString titleBetaUp=Form("betaSystUp%c_%d",Q,year);
    betaSystUp = new TH1F(titleBetaUp,titleBetaUp,maxEta,0.5,maxEta+0.5);
    TString titleBetaDn=Form("betaSystDn%c_%d",Q,year);
    betaSystDn = new TH1F(titleBetaDn,titleBetaDn,maxEta,0.5,maxEta+0.5);
    TString titleZfrac=Form("zFracSecondEEMC%c_%d",Q,year);
    zFracSecondEEMC = new TH1F(titleZfrac,titleZfrac,maxEta,0.5,maxEta+0.5);
    
    cc = new TCanvas("cc","cc",1200,800);
    cc->Divide(4,2);
    dd = new TCanvas("dd","dd",1200,800);
    dd->Divide(4,2);

    for(int ieta=0; ieta<maxEta; ieta++){
      if(ieta>=4 && ieta<=6) continue;

      cout<<endl<<"----------------------------------"<<endl;
      cout<<"  Background for W"<<Q<<": Eta"<<ieta+1<<endl;
      cout<<"----------------------------------"<<endl<<endl;
      plBkgdEtaDepX(charge,Form("Eta%d",ieta+1),ieta+1,false,year);
    }
    outSummary<<endl<<endl;

    // add hard coded endcap values from separate macro
    if(charge==1) {
      beta->SetBinContent(7,0.898);
      beta->SetBinError(7,0.030);
    }
    else if(charge==-1) {
      beta->SetBinContent(7,0.903);
      beta->SetBinError(7,0.016);
    }

    //write histos to output file
    if(outF->IsOpen()){
      outF->cd();
      for(iF=0; iF<maxF; iF++){
	bkgdFrac[iF]->Write();
      }
      beta->Write();
      zFracSecondEEMC->Write();
    }
    
    outF->Close();

    cc->Print(Form("bkgdET_%d_%c.png",year,Q));
    cc->Print(Form("bkgdET_%d_%c.ps",year,Q));
    dd->Print(Form("beta_%d_%c.png",year,Q));
    dd->Print(Form("beta_%d_%c.ps",year,Q));
  }

  outSummary.close();
  
  return;
}


// ********************************************************
// This code makes the necessary root histos for the 
// background subtraction and systematic uncertainty calculations
// Set charge == 1 for postive charges
// Set charge == -1 for negative charges
// Set charge == 0 for the sum of both (not tested)
// ********************************************************

void plBkgdEtaDepX(int charge=1, char* etaBin="Eta8", int starPhysEtaBin=8, bool secondEEMC=false, int year) {

  // ******************************************************
  // Read in the data set and all the histograms needed to do the actual calculation
  // ******************************************************

  TString dataDir="/star/u/devika/wAnalysis/2012/makeALplots/testwithrun13test/"; dataDir+=year; dataDir+="/";
  string dataName;
  if(year==2012) dataName = "53eta.wana.hist.root";
  else if(year==2011) dataName = "run11long_P11id.wana.hist.root";

  TFile *f1 = new TFile(Form("%s%s",dataDir.Data(),dataName));
  f1->cd(etaBin); //needed for eta dependent

  // get the signal w/ and w/o the EEMC in the algo
  if (charge == 1) {
    TH1F *signal = (TH1F*)gDirectory->Get("pos_muclustpTbal_wE");
    TH1F *signal_wo_eemc = (TH1F*)gDirectory->Get("pos_muclustpTbal_noE");
  } else if (charge == -1) {
    TH1F *signal = (TH1F*)gDirectory->Get("neg_muclustpTbal_wE");
    TH1F *signal_wo_eemc = (TH1F*)gDirectory->Get("neg_muclustpTbal_noE");
  } else if (charge == 0) {
    TH1F *signal = (TH1F*)gDirectory->Get("muclustPtBal");
    TH1F *signal_wo_eemc = (TH1F*)gDirectory->Get("muclustPtBalnoE");
  }

  // ******************************************************
  // Read in all the MC data sets for background subtractions and studies
  // ******************************************************

  TFile *MC_fs[6];
  //  TString mcDir="/star/institutions/mit/stevens4/freezer/2012-W-PRL-9.12.13/GPC/"; mcDir+=year; mcDir+="/";
  TString mcDir="/star/u/devika/wAnalysis/2012/makeALplots/testwithrun13test/"; mcDir+=year; mcDir+="/";

  MC_fs[0] = new TFile(Form("%sWplus_enu.wana.hist.root",mcDir.Data())); // W+ -> e++nu
  MC_fs[1] = new TFile(Form("%sWminus_enu.wana.hist.root",mcDir.Data())); // W- -> e-+nu
  MC_fs[2] = new TFile(Form("%sWplus_taunu.wana.hist.root",mcDir.Data())); // W+ -> tau+nu
  MC_fs[3] = new TFile(Form("%sWminus_taunu.wana.hist.root",mcDir.Data())); // W- -> tau+nu
  MC_fs[4] = new TFile(Form("%sZ_eplus_eminus_inter.wana.hist.root",mcDir.Data())); // Z/gamma* -> e+ e-

  //get eta dependent signal and background from MC
  TH1F *MC_dists_raw[5][3];
  for (int i=0; i<5; i++) {
    MC_fs[i]->cd(etaBin);
    if (charge == 1) {
      MC_dists_raw[i][0] = (TH1F*)gDirectory->Get("pos_muclustpTbal_wE");
      MC_dists_raw[i][1] = (TH1F*)gDirectory->Get("pos_muclustpTbal_noE");
    } else if (charge == -1) {
      MC_dists_raw[i][0] = (TH1F*)gDirectory->Get("neg_muclustpTbal_wE");
      MC_dists_raw[i][1] = (TH1F*)gDirectory->Get("neg_muclustpTbal_noE");
    } else if (charge == 0) {
      MC_dists_raw[i][0] = (TH1F*)gDirectory->Get("muclustPtBal");
      MC_dists_raw[i][1] = (TH1F*)gDirectory->Get("muclustPtBalnoE");
    }
    MC_dists_raw[i][2] = (TH1F*)gDirectory->Get("muclustPtBal_bckgrd");
  }

  float WplusXsec  = 98.5;
  float WminusXsec = 31.3;
  float ZgammaXsec = 23.9;
  float WplusEvents, WminusEvents, ZgammaEvents;
  float dataLumi;
  if(year==2012) { 
    WplusEvents = 109051; WminusEvents = 34936.;
    // WplusEvents = 80000; WminusEvents = 20000.; 
    ZgammaEvents = 25126.; dataLumi = 77.4; //From Jamie's numbers 77.4 for Run 12 (638 runs)
    //ZgammaEvents = 7000.; dataLumi = 40.4;
  }
  else if(year==2011) { 
    WplusEvents = 11919; WminusEvents = 3816.; 
    ZgammaEvents = 2749.; dataLumi = 9.4; //From Jamie's numbers 9.4 for Run 11 longitudinal (? runs)
  }

  //pure lumi from pythia to scale MC
  float lumi[5] = {WplusEvents/WplusXsec,WminusEvents/WminusXsec,WplusEvents/WplusXsec,WminusEvents/WminusXsec,ZgammaEvents/ZgammaXsec};

  //scale MC to match lumi of the dataset
  float lumi_fact[6];
  for (int i=0; i<5; i++) {lumi_fact[i] = dataLumi/lumi[i];} 

  for (int i=0; i<5; i++) {
    for (int j=0; j<3; j++) {
      MC_dists_raw[i][j]->Scale(lumi_fact[i]);
    }
  }

  // Repack the MC histograms to mesh with the odd staggered binning that is being used
  char str[200];
  TH1F *MC_dists_repack[5][3];
  for (int i=0; i<5; i++) {
    MC_dists_repack[i][0]=(TH1F*) MC_dists_raw[i][0]->Clone(); 
    MC_dists_repack[i][1]=(TH1F*) MC_dists_raw[i][1]->Clone();
    MC_dists_repack[i][2]=(TH1F*) MC_dists_raw[i][2]->Clone(); 
  }

  // Make second endcap background distributions for the MC
  for (int i=0; i<5; i++) {
    MC_dists_repack[i][1]->Add(MC_dists_repack[i][0],-1);
    MC_dists_raw[i][1]->Add(MC_dists_raw[i][0],-1.);
  }
  TH1F *zeemc = MC_dists_raw[4][1]->Clone();

  // Make second endcap background distribution for the data
  TH1F *eemc_bkgd = signal_wo_eemc->Clone();
  eemc_bkgd->Add(signal,-1.); //this is the '2nd endcap' bkgd

  // first pass of macro sets second EEMC background to use for mirror etaBin
  if(secondEEMC) {
    if(charge>0) {
      mirrorSecondEEMC[0][starPhysEtaBin-1] = eemc_bkgd;
      mirrorZSecondEEMC[0][starPhysEtaBin-1] = zeemc;
    }
    else {
      mirrorSecondEEMC[1][starPhysEtaBin-1] = eemc_bkgd;
      mirrorZSecondEEMC[1][starPhysEtaBin-1] = zeemc;
    }
    return;
  }

  // Now eemc_bkgd and zeemc is for the mirror etaBin
  int bin = starPhysEtaBin-1;
  eemc_bkgd = charge > 0 ? mirrorSecondEEMC[0][mirrorEtaBin[bin]-1] : mirrorSecondEEMC[1][mirrorEtaBin[bin]-1];
  zeemc = charge > 0 ? mirrorZSecondEEMC[0][mirrorEtaBin[bin]-1] : mirrorZSecondEEMC[1][mirrorEtaBin[bin]-1];

  // **********************************************
  // **********************************************
  // Repack all necessary histograms
  // **********************************************
  // **********************************************

  //all histos rebinned to match binning starting at 25 GeV
  TH1F *zsig = MC_dists_raw[4][0]->Clone();
  TH1F *zback = MC_dists_raw[4][2]->Clone();
  zsig->Add(zeemc,-1.); // remove second EEMC Zs from the Z MC in the signal
  for(int i=0; i<zsig->GetNbinsX(); i++){
    if(zsig->GetBinContent(i)<0.)  zsig->SetBinContent(i,0.);
  }

  TH1F *wsig;
  if(charge>0) wsig = (TH1F*)MC_dists_raw[0][0]->Clone();
  else if(charge<0) wsig = (TH1F*)MC_dists_raw[1][0]->Clone();
  else {
    wsig = (TH1F*)MC_dists_raw[0][0]->Clone();
    wsig->Add(MC_dists_raw[1][0]);
  }

  TH1F *tausig;
  if(charge>0) tausig = (TH1F*)MC_dists_raw[2][0]->Clone();
  else if(charge<0) tausig = (TH1F*)MC_dists_raw[3][0]->Clone();
  else {
    tausig = (TH1F*)MC_dists_raw[2][0]->Clone();
    tausig->Add((TH1F*)MC_dists_raw[3][0]);
  }

  // Subtract second endcap background from the signal distribution
  TH1F *signal_final = signal->Clone();
  signal_final->Add(eemc_bkgd,-1); 
  // Subtract off MC Z background from the signal distribution
  signal_final->Add(zsig,-1.);
  // Subtract off MC tau background from the signal distribution
  signal_final->Add(tausig,-1);

  // First get the "nominal" QCD background shape (charge summed)
  f1->cd(etaBin); //switch back to data file for gDirectory
  TH1F *bkgd_shape1 = (TH1F*)gDirectory->Get("muclustPtBal_bckgrd");
  TH1F *bkgd_shape_nom = (TH1F*)bkgd_shape1->Clone();

  TH1F *eemc_bkgd2 = (TH1F*)eemc_bkgd->Clone(); 
  TH1F *tausig_bkgd2 = (TH1F*)tausig->Clone();
  TH1F *zsig_bkgd2 = (TH1F*)zsig->Clone();
  TH1F *zeemc_bkgd2 = (TH1F*)zeemc->Clone(); 
  TH1F *zback_bkgd2 = (TH1F*)zback->Clone(); 

  TH1F *signal_final2 = signal_final->Clone(); 
  signal_final2->SetLineColor(2); signal_final2->SetLineWidth(2.*signal_final2->GetLineWidth());
  TH1F *signal2 = (TH1F*)signal->Clone(); 
  TH1F *wsig2 = (TH1F*)wsig->Clone(); 

  //use W MC for signal in normalization window
  TH1F *signal_in_norm_region = (TH1F*)wsig->Clone();

  // ******************************************************************************
  // ******************************************************************************
  // Do the normalization of the W signal to the background for the nominal shape
  // ******************************************************************************
  // ******************************************************************************
  
  TH1F *bkgd_shape_unnorm = (TH1F*)bkgd_shape_nom->Clone();
  // subtract off the remaining Z events from the background shape
  bkgd_shape_unnorm->Add(zback_bkgd2,-1.);
  // make sure background shape is non-negative (caused by Z subtraction)
  for (int j=1; j<=bkgd_shape_unnorm->GetNbinsX(); j++) {
    if (bkgd_shape_unnorm->GetBinContent(j) < 0) bkgd_shape_unnorm->SetBinContent(j,0.);
  }

  // calculate the normalization factor
  int minBin = bkgd_shape_unnorm->FindBin(14.01);
  int maxBin = bkgd_shape_unnorm->FindBin(15.99);
  float normt = 0, normb = 0.;
  for (int k=minBin; k<=maxBin; k++) {
    if (bkgd_shape_unnorm->GetBinContent(k) > 0) {
      normt += signal_final2->GetBinContent(k)-signal_in_norm_region->GetBinContent(k);
      normb += bkgd_shape_unnorm->GetBinContent(k);
    }
  }

  //determine normalization for background shape
  float new_bkgd_norm;
  if (normb > 0 && normt > 0) {
    float norm = normt/normb;
    bkgd_shape_unnorm->Scale(norm);
    new_bkgd_norm=norm;
  }
  else if(normb > 0 && normt < 0){
    float norm = 0;
    bkgd_shape_unnorm->Scale(norm);
    new_bkgd_norm=norm;
  }
 
  //final QCD background
  TH1F *new_bkgd;
  new_bkgd = (TH1F*)bkgd_shape_unnorm->Clone();
  new_bkgd->SetName("new_bkgd");


  // ******************************************************
  // ******************************************************
  // Calculate all the 60 shapes for the background systematic uncertainty
  // ******************************************************
  // ******************************************************

  // get the various background shapes made by using different parameters for sPtBal cut
  TH2F *bkgd_hists_from_file;
  TH2F *bkgd_hists_from_file2;
  char str[200];
  sprintf(str,"pos_failAwaySide_Awayside_pt_bin_%d",0);
  bkgd_hists_from_file = (TH2F*)gDirectory->Get(str);
  sprintf(str,"neg_failAwaySide_Awayside_pt_bin_%d",0);
  bkgd_hists_from_file2 = (TH2F*)gDirectory->Get(str);
  bkgd_hists_from_file->Add(bkgd_hists_from_file2);

  // make array of TH1F's from the 2D histogram above
  const int maxSteps = 81;
  TH1F *bkgd_hists[maxSteps];
  for (int i=0; i<maxSteps; i++) {
    sprintf(str,"bkgd_hist_%d",i);
    bkgd_hists[i] = new TH1F(str,str,bkgd_hists_from_file->GetNbinsX(),0,100);
    for (int k=1; k<=100; k++) 
      bkgd_hists[i]->SetBinContent(k,bkgd_hists_from_file->GetBinContent(k,i+1));
  }

  // loop over the possible combinations of normalizationi ranges and background shapes
  const int maxNorm = 10;
  int normLow = bkgd_hists[0]->FindBin(14.01);
  int normHigh[maxNorm]; //upper end of normalization region in ET (16, 17, 18, 19, 20)
  for(int i=0; i<maxNorm; i++) normHigh[i] = bkgd_hists[0]->FindBin(15.99+i*0.5); 
  TH1F *new_bkgd_hists[maxSteps][maxNorm];
  float new_bkgd_norm[maxSteps][maxNorm];
  
  for (int i=0; i<maxSteps; i++) { //loop over possible sPtBalance cuts
    for(int j=0; j<maxNorm; j++){
      TH1F *bkgd_shape_unnorm;
      bkgd_shape_unnorm = (TH1F*)bkgd_hists[i]->Clone();
      
      // subtract off the Z contamination to background shape
      bkgd_shape_unnorm->Add(zback_bkgd2,-1.);
    
      // make sure background shape is non-negative (caused by Z subtraction)
      for (int m=1; m<=100; m++) 
	if (bkgd_shape_unnorm->GetBinContent(m) < 0) {bkgd_shape_unnorm->SetBinContent(m,0.);}
    
      // calculate the normalization factor for 1 bin
      float normt = 0, normb = 0.;
      for (int k=15; k<=normHigh[j]; k++) {
	if (bkgd_shape_unnorm->GetBinContent(k) > 0) {
	  normt += signal_final2->GetBinContent(k)-signal_in_norm_region->GetBinContent(k);
	  normb += bkgd_shape_unnorm->GetBinContent(k);
	}
      }
      if (normb > 0 && normt > 0) {
	float norm = normt/normb;
	bkgd_shape_unnorm->Scale(norm);
      }
      else if(normb > 0 && normt < 0){
	float norm = 0;
	bkgd_shape_unnorm->Scale(norm);
      }
      
      // save each of the background distributions for syst study later
      new_bkgd_hists[i][j] = (TH1F*)bkgd_shape_unnorm->Clone();
      new_bkgd_norm[i][j] = (normb > 0 && normt > 0) ? normt/normb : 0.;
    }
  }

#if 0
  // plot all the QCD background distributions
  signal_final2->SetStats(kFALSE);
  if (charge == 1) {
    signal_final2->SetTitle("W+ Background Shapes");
  } else if (charge == -1) {
    signal_final2->SetTitle("W- Background Shapes");
  }
  signal_final2->GetXaxis()->SetRangeUser(0.,70.);
  signal_final2->GetXaxis()->SetTitle("2x2 Cluster E_{T} (GeV)");
  signal_final2->Draw();

  //eemc_bkgd->Draw("same"); eemc_bkgd->SetLineColor(kRed);
  //zsig->Draw("same"); zsig->SetLineColor(kGreen);
  //tausig->Draw("same"); tausig->SetLineColor(kMagenta);
  
  for (int i=0; i<maxSteps; i++) {
    for(int j=0; j<maxNorm; j++){
      new_bkgd_hists[i][j]->Draw("same");
    }
  }
  new_bkgd->SetLineColor(4);
  new_bkgd->Draw("same");
  return;
#endif

  // ************************************************
  // Calculate all background numbers and their 
  // systematic uncertainties (and spit them out to histograms)
  // ************************************************

  // First get the simple numbers (backgrounds and their statistical uncertainties)
  float tau_norm = lumi_fact[2];
  float Z_norm = lumi_fact[4];
  float bkgd_sum = 0.;
  float signal_sum = 0.;
  float raw_sum = 0.;
  float QCD_sum = 0., tau_sum = 0., eemc_sum = 0.;
  float QCD_raw_sum = 0.;
  float zsig_sum = 0., zeemc_sum = 0.,zback_sum = 0.;
  //get counts in ET>25 && ET<50 for different contributions
  int minBin25 = signal_final->FindBin(25.01);
  int maxBin50 = signal_final->FindBin(49.99);
  for (int i=minBin25; i<=maxBin50; i++) { 
    bkgd_sum += new_bkgd->GetBinContent(i);
    bkgd_sum += tausig_bkgd2->GetBinContent(i);
    bkgd_sum += eemc_bkgd2->GetBinContent(i);
    QCD_raw_sum += bkgd_shape_nom->GetBinContent(i);
    QCD_sum += new_bkgd->GetBinContent(i);
    tau_sum += tausig_bkgd2->GetBinContent(i);
    eemc_sum += eemc_bkgd2->GetBinContent(i);
    signal_sum += signal_final2->GetBinContent(i);
    raw_sum += signal2->GetBinContent(i);
    zsig_sum += zsig_bkgd2->GetBinContent(i);
    zeemc_sum += zeemc_bkgd2->GetBinContent(i);
    zback_sum += zback_bkgd2->GetBinContent(i);
  }
  cout << "The total background for 25 < ET < 50 is " << bkgd_sum+zsig_sum << endl;
  cout << "QCD = " << QCD_sum << ", tau = " << tau_sum << ", eemc = " << eemc_sum << ", and Z = " << zsig_sum << endl;
  cout << "Raw = " << raw_sum << endl;
  cout << "W Signal (w/o tau) = " << signal_sum-QCD_sum << endl;
  cout << "Z in sig = " << zsig_sum << endl;
  cout << "Z in eemc = " << zeemc_sum << endl;
  cout << "Z in back = " << zback_sum << endl;
  cout << "QCD raw in back = " << QCD_raw_sum << endl; 
  cout << "The QCD stat unc. is " << sqrt(new_bkgd_norm*QCD_sum) << endl;

  float QCD_stat = new_bkgd_norm*sqrt(QCD_sum/new_bkgd_norm);
  float tau_stat = tau_norm*sqrt(tau_sum/tau_norm);
  float tau_syst = 0.13*tau_sum;
  cout << "The tau stat unc. is " << tau_stat << " syst is " << tau_syst << " and total is " << sqrt(tau_stat*tau_stat + tau_syst*tau_syst) << endl;
  float eemc_stat = feldmanCousins(eemc_sum);
  float eemc_syst = 0.13*zeemc_sum;
  cout << "The eemc stat unc. is " << eemc_stat << " syst is " << eemc_syst << " and total is " << sqrt(eemc_stat*eemc_stat + eemc_syst*eemc_syst) << endl;
  float zsig_stat = (zsig_sum>0) ? Z_norm*sqrt(zsig_sum/Z_norm) : 0; 
  float zsig_syst = 0.13*zsig_sum;
  cout << "The Z stat unc. is " << zsig_stat << " syst is " << zsig_syst << " and total is " << sqrt(zsig_stat*zsig_stat + zsig_syst*zsig_syst) << endl;
   
  //find fs and beta
  float f_tau = tau_sum/raw_sum;
  float f_z = zsig_sum/raw_sum;
  float f_eemc = eemc_sum/raw_sum;
  float f_QCD = QCD_sum/raw_sum;
  float beta_sum = 1. - f_z - f_eemc - f_QCD;
  //find error on fs and beta
  float ferr_tau = (tau_sum>0) ? f_tau*sqrt(tau_stat*tau_stat/tau_sum/tau_sum + 1./raw_sum) : 0;
  float ferr_z = (zsig_sum>0) ? f_z*sqrt(zsig_stat*zsig_stat/zsig_sum/zsig_sum + 1./raw_sum) : 0;
  float ferr_eemc = (eemc_sum>0) ? f_eemc*sqrt(eemc_stat*eemc_stat/eemc_sum/eemc_sum + 1./raw_sum) : 0;
  float ferr_QCD = (QCD_sum>0) ? f_QCD*sqrt(QCD_stat*QCD_stat/QCD_sum/QCD_sum + 1./raw_sum) : 0;
  float beta_stat = sqrt( ferr_z*ferr_z + ferr_eemc*ferr_eemc + ferr_QCD*ferr_QCD );
  
  bkgdFrac[0]->SetBinContent(starPhysEtaBin,f_tau);
  bkgdFrac[0]->SetBinError(starPhysEtaBin,ferr_tau);
  bkgdFrac[1]->SetBinContent(starPhysEtaBin,f_z);
  bkgdFrac[1]->SetBinError(starPhysEtaBin,ferr_z);
  bkgdFrac[2]->SetBinContent(starPhysEtaBin,f_eemc);
  bkgdFrac[2]->SetBinError(starPhysEtaBin,ferr_eemc);
  bkgdFrac[3]->SetBinContent(starPhysEtaBin,f_QCD);
  bkgdFrac[3]->SetBinError(starPhysEtaBin,ferr_QCD);
  zFracSecondEEMC->SetBinContent(starPhysEtaBin,(eemc_sum>0) ? zeemc_sum/eemc_sum : 0);
  
  cout << "f_tau = " << f_tau << endl;
  cout << "f_Z = " << f_z << endl;
  cout << "f_eemc = " << f_eemc << endl;
  cout << "f_QCD = " << f_QCD << endl;

  // ********** Old method max-extent uncertainty **********
  // Now go through all the background shapes and find the
  // high and low in each ET bin to give the maximum extent uncertainty

  TH1F *QCD_syst_high_err = new TH1F("QCD_syst_high_err","QCD_syst_high_err",signal_final->GetNbinsX(),0.,100.);
  TH1F *QCD_syst_low_err = new TH1F("QCD_syst_low_err","QCD_syst_low_err",signal_final->GetNbinsX(),0.,100.);

  TH1F *low_bkgd = new_bkgd->Clone();
  low_bkgd->SetName("low_bkgd");
  TH1F *high_bkgd = new_bkgd->Clone();
  high_bkgd->SetName("high_bkgd");

  TH2F *minCounter = new TH2F("minCounter","minCounter; ;shape index",signal_final->GetNbinsX(),0,100.,maxSteps,0,maxSteps);
  TH2F *maxCounter = new TH2F("maxCounter","maxCounter; ;shape index",signal_final->GetNbinsX(),0,100.,maxSteps,0,maxSteps);

  float low_sum = 0.;
  float high_sum = 0.;
  int minBin = signal_final->FindBin(14.01);
  int maxBin = signal_final->FindBin(99.99);
  for (int i=minBin; i<=maxBin; i++) {
    float high = 0.;
    float low = 10000.;
    int imin=61; int imax=61;
    for (int j=0; j<maxSteps; j++) {
      for (int k=0; k<maxNorm; k++) {
	if (new_bkgd_hists[j][k]->GetBinContent(i) < low) {
	  if (new_bkgd_hists[j][k]->GetBinContent(i) >= 0) {
	    low = new_bkgd_hists[j][k]->GetBinContent(i); imin=j+k*20;
	  }
	}
	if (new_bkgd_hists[j][k]->GetBinContent(i) > high) {
	  high = new_bkgd_hists[j][k]->GetBinContent(i); imax=j+k*20;
	}
      } // end of k-loop
    } // end of j-loop
    
    maxCounter->Fill(i,imax);
    minCounter->Fill(i,imin);

    // calculate the sum
    low_bkgd->SetBinContent(i,0.);
    if ((low != 10000) && (new_bkgd->GetBinContent(i)-low > 0)) {
      if ((i>=minBin25)&&(i<=maxBin50)) {low_sum += low;}
      low_bkgd->SetBinContent(i,low);
    }
    if ((i>=minBin25)&&(i<=maxBin50)) {high_sum += high;} 
    high_bkgd->SetBinContent(i,high);

    // set the bin-by-bin error too
    if ((low != 10000) && (new_bkgd->GetBinContent(i)-low > 0)) {
      QCD_syst_low_err->SetBinContent(i,new_bkgd->GetBinContent(i)-low);
    } else {
      QCD_syst_low_err->SetBinContent(i,0.);
    }
    QCD_syst_high_err->SetBinContent(i,high-new_bkgd->GetBinContent(i));
  } // end of i-loop 

  cout << "QCD shape sys. unc. calc************" << endl;
  cout << "The QCD low sum = " << low_sum << endl;
  cout << "The QCD high sum = " << high_sum << endl; 

  cout << "The QCD low error = " << QCD_sum-low_sum << endl;
  cout << "The QCD high error = " << high_sum-QCD_sum << endl; 
  
  // ********** New method for systematic variation of beta with shape/norm region **********
  // get beta for each of the possible shape/normalizations combinations
  float beta_width = 0.025;
  if(charge>=0) beta_width = 0.015;
  float minBeta = beta_sum - beta_width;
  float maxBeta = beta_sum + beta_width;
  TH1F* hBetaPDF = new TH1F("hBetaPDF","hBetaPDF",50,minBeta,maxBeta);
  for (int j=0; j<maxSteps; j++) {
    for (int k=1; k<maxNorm; k++) {
      float new_bkgd_sum = 0;
      
      for (int i=minBin25; i<=maxBin50; i++) 
	new_bkgd_sum += new_bkgd_hists[j][k]->GetBinContent(i);
      
      // fill many values for beta here
      //float weight = (new_bkgd_sum>0.) ? 1./(new_bkgd_sum*new_bkgd_norm[j][k]) : 0; //weight by QCD statistical error
      hBetaPDF->Fill(1. - f_z - f_eemc - new_bkgd_sum/raw_sum); 
    }    
  }
  
  dd->cd(starPhysEtaBin);
  hBetaPDF->Draw();
  TLine* ln = new TLine(beta_sum,0,beta_sum,hBetaPDF->GetMaximum());
  ln->SetLineColor(kBlue); ln->Draw("same");
  //return;
  
  // average beta from varying QCD shape and normalization
  cout << "QCD systematic error from beta RMS = " << hBetaPDF->GetRMS() << endl;
  beta->SetBinContent(starPhysEtaBin,hBetaPDF->GetMean());
  beta->SetBinError(starPhysEtaBin,sqrt(beta_stat*beta_stat + hBetaPDF->GetRMS()*hBetaPDF->GetRMS()));
  betaSyst->SetBinContent(starPhysEtaBin,hBetaPDF->GetRMS());

  // fill QCD systematic ********* old method ***********
  float f_QCDmax = high_sum/raw_sum;
  float f_QCDmin = low_sum/raw_sum;
  float beta_max = 1. - f_z - f_eemc - f_QCDmin;
  float beta_min = 1. - f_z - f_eemc - f_QCDmax;
  betaSystUp->SetBinContent(starPhysEtaBin,beta_max);
  betaSystDn->SetBinContent(starPhysEtaBin,beta_min);

  outSummary<<Form("%d, %.00f, %.02f, %02d, %.02f, %.02f, %.02f, %.02f, %.03f +/- %.03f, %.03f, %.03f",starPhysEtaBin,raw_sum,QCD_sum,int(eemc_sum),zsig_sum,tau_sum,QCD_sum+eemc_sum+zsig_sum,raw_sum-(QCD_sum+eemc_sum+zsig_sum),beta_sum,beta_stat,beta_max-beta_sum,beta_sum-beta_min)<<endl;

  //*******************************************
  //********* plot the jacobian peak **********
  //*******************************************
  cc->cd(starPhysEtaBin);

  // repack histos so start at 14 GeV with 4 GeV wide bins
  TH1F* zsig_bkgd_repack = new TH1F("zsig_bkgd_repack","zsig_bkgd_repack",50,2.,102.);
  TH1F* wsig_repack = new TH1F("wsig_repack","wsig_repack",50,2.,102.);
  TH1F* tausig_bkgd_repack = new TH1F("tausig_repack","tausig_repack",50,2.,102.);
  TH1F* eemc_bkgd_repack = new TH1F("eemc_bkgd_repack","eemc_bkgd_repack",50,2.,102.);
  TH1F* new_bkgd_repack = new TH1F("new_bkgd_repack","new_bkgd_repack",50,2.,102.);
  TH1F* signal2_repack = new TH1F("signal2_repack","signal2_repack",50,2.,102.);

#if 1
  zsig_bkgd2->Rebin(2);
  wsig2->Rebin(2);
  tausig_bkgd2->Rebin(2);
  eemc_bkgd2->Rebin(2);
  new_bkgd->Rebin(2);
  signal2->Rebin(2);
#endif

  for (int k=1; k<=50; k++) {
    zsig_bkgd_repack->SetBinContent(k,zsig_bkgd2->GetBinContent(2*k+1)+zsig_bkgd2->GetBinContent(2*k+2));
    wsig_repack->SetBinContent(k,wsig2->GetBinContent(2*k+1)+wsig2->GetBinContent(2*k+2));
    tausig_bkgd_repack->SetBinContent(k,tausig_bkgd2->GetBinContent(2*k+1)+tausig_bkgd2->GetBinContent(2*k+2));
    eemc_bkgd_repack->SetBinContent(k,eemc_bkgd2->GetBinContent(2*k+1)+eemc_bkgd2->GetBinContent(2*k+2));
    new_bkgd_repack->SetBinContent(k,new_bkgd->GetBinContent(2*k+1)+new_bkgd->GetBinContent(2*k+2));
    signal2_repack->SetBinContent(k,signal2->GetBinContent(2*k+1)+signal2->GetBinContent(2*k+2));
  }

  //4 GeV bins
  THStack *hs = new THStack("hs",";E_{T}^{e} (GeV);Counts");
  int rebin=1;
  zsig_bkgd_repack->SetFillColor(kViolet+9); zsig_bkgd_repack->Rebin(rebin); hs->Add(zsig_bkgd_repack);
  tausig_bkgd_repack->SetFillColor(49); tausig_bkgd_repack->Rebin(rebin); hs->Add(tausig_bkgd_repack);
  eemc_bkgd_repack->SetFillColor(30); eemc_bkgd_repack->Rebin(rebin); hs->Add(eemc_bkgd_repack);
  new_bkgd_repack->SetFillColor(kCyan-10);  new_bkgd_repack->Rebin(rebin); hs->Add(new_bkgd_repack);
  wsig_repack->SetLineColor(kRed); wsig_repack->SetLineWidth(2); wsig_repack->SetLineStyle(2);
  wsig_repack->Rebin(rebin);  wsig_repack->GetXaxis()->SetRangeUser(10,70);
  hs->Add(wsig_repack); 

  signal2_repack->Rebin(rebin); signal2_repack->GetXaxis()->SetRangeUser(10,70);
  signal2_repack->SetTitle(Form("%s: %s; E_{T} (GeV)",signal2->GetTitle(),etaBin)); 

  float max;
  if(year==2012) {
    max = 250;  
    if(charge==1 && starPhysEtaBin==8) max=800;
    if(charge==-1 && starPhysEtaBin<8) max=135;
    if(charge==-1 && starPhysEtaBin==8) max=420;
  }
  else if(year==2011) {
    float max = 50;  
    if(charge==1 && starPhysEtaBin==8) max=110;
    if(charge==-1 && starPhysEtaBin<8) max=35;
    if(charge==-1 && starPhysEtaBin==8) max=80;
  }

  signal2_repack->SetMaximum(max);

  tx=new TLatex(60,max*0.4,"W+");   
  if(charge==-1) tx=new TLatex(60,max*0.4,"W-"); 
  tx->SetTextAlign(21); tx->SetTextSize(0.1);
  
  gPad->SetGridx(false); gPad->SetGridy(false); 
  signal2_repack->Draw("he");
  hs->Draw("same");

  TLegend *leg = new TLegend(0.5,0.63,0.89,0.89);
  leg->SetEntrySeparation(0.01);
  leg->SetMargin(0.2);
  leg->SetFillColor(0);
  leg->SetTextSize(0.04);
  leg->SetLineColor(kWhite);
  leg->AddEntry(signal2_repack," STAR 2012 Data","el");
  leg->AddEntry(wsig_repack," W #rightarrow e #nu MC","l");
  leg->AddEntry(new_bkgd_repack," Data-driven QCD","f");
  leg->AddEntry(eemc_bkgd_repack," Second EEMC","f");
  leg->AddEntry(tausig_bkgd_repack," W #rightarrow #tau #nu MC","f");
  leg->AddEntry(zsig_bkgd_repack," Z #rightarrow ee MC","f");
  leg->Draw();
  tx->Draw("same");
  
  // write histos to output file for pub quality plot
    TDirectory *dir;
    if(charge==1) dir = outF->mkdir(Form("pubPos%d",starPhysEtaBin));
    else if(charge==-1) dir = outF->mkdir(Form("pubNeg%d",starPhysEtaBin));
    dir->cd();
    signal2_repack->Write();
    wsig_repack->Write();
    new_bkgd_repack->Write();
    eemc_bkgd_repack->Write();
    tausig_bkgd_repack->Write();
    zsig_bkgd_repack->Write();

  return;  
}


// 1-sigma CL upper bound for poisson process with small number of observed events 
// from Feldman-Cousins paper -> arXiv:physics/9711021
float feldmanCousins(int x){

  if(x==0) return 1.29;
  else if(x==1) return 1.75;
  else if(x==2) return 2.25;
  else if(x==3) return 2.30;
  else if(x==4) return 2.78;
  else if(x==5) return 2.81;
  else if(x==6) return 3.28;
  else if(x==7) return 3.30;
  else if(x==8) return 3.32;
  else if(x==9) return 3.79;

  return sqrt(x);
}