///////////////////////////////////////////////////////////////////////////////
//  Oct 2003 P.Murat
//  plot results of PDF error analysis - so far CTEQ6
//  
///////////////////////////////////////////////////////////////////////////////

TH1F*  gH1;
TH2F*  gH2;

//_____________________________________________________________________________
void get_pdf_err(TString Process, double Q[]) {

  const char module    [   ] = "TauMcAna";
  char       get_hist_c[200], filename[200];

  const char* work_dir = gSystem->Getenv("SRT_LOCAL");
  sprintf(get_hist_c,"%s/plot/get_hist.C",work_dir);

  TH2F    *h2;
  TH1F    *h1;
  TString hist_name;

  if (!gInterpreter->IsLoaded(get_hist_c)) gInterpreter->LoadMacro(get_hist_c);

  char format[200];

  for (int i=0; i<41; i++) { 
    
    if (Process == "wenu") {
      strcpy(format,"hist/dev_240/pythia_ana/wenu/wenu.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h2   = (TH2F*) get_hist(filename,module,"genp/ele_nl_vs_nt_1");
      Q[i] = h2->GetBinContent(2,2);
    }
    else if (Process == "wmunu") {
      strcpy(format,"hist/dev_240/pythia_ana/wmunu/wmunu.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h2   = (TH2F*) get_hist(filename,module,"genp/muo_nl_vs_nt_1");
      Q[i] = h2->GetBinContent(2,2);
    }
    else if (Process == "wtaunu") {
      strcpy(format,"hist/dev_240/pythia_ana/wtaunu/wtaunu.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h1   = get_hist(filename,module,"genp/tau_nt_1");
      Q[i] = h1->GetBinContent(2);
    }
    else if (Process == "zee") {
      strcpy(format,"hist/dev_240/pythia_ana/zee/zee.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h1   = get_hist(filename,module,"genp/ele_nl_vs_nt_0");
      Q[i] = h1->GetBinContent(3,3)+h1->GetBinContent(2,3);
    }
    else if (Process == "zee_3") {
      strcpy(format,"hist/dev_240/pythia_ana/zee_3/zee.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h1   = get_hist(filename,module,"genp/ele_nl_vs_nt_0");
      Q[i] = h1->GetBinContent(3,3)+h1->GetBinContent(2,3);
    }
    else if (Process == "zmumu") {
      strcpy(format,"hist/dev_240/pythia_ana/zmumu/zmumu.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h1   = get_hist(filename,module,"genp/muo_nl_vs_nt_0");
      Q[i] = h1->GetBinContent(3,3)+h1->GetBinContent(2,3);
    }
    else if (Process == "zmumu_d0") {
      strcpy(format,"hist/dev_240/pythia_ana/zmumu_d0/zmumu.pythia_ana.%i.1000000.root");
      sprintf(filename,format,30200+i);
      h1   = get_hist(filename,module,"genp/muo_nl_vs_nt_0");
      Q[i] = h1->GetBinContent(3,3)+h1->GetBinContent(2,3)+h1->GetBinContent(1,3);
    }
  }
}

//_____________________________________________________________________________
void plot_pdf_err(TString Process) {

  double q[41], qq, sum;

  get_pdf_err(Process,q);
//-----------------------------------------------------------------------------
// now plot histogram with variations
//-----------------------------------------------------------------------------
  TH1F* h1 = new TH1F(Form("pdf_err_%s",Process.Data()),
			   Form("pdf error for %s",Process.Data()),
				100,0.95,1.05);

  for (int i=1; i<=40; i++) {
    qq = q[i]/q[0];
    h1->Fill(qq);
  }

  h1->Draw();

  sum = 0;
  for (int i=1; i<=20; i++) {
    qq = (q[2*i]-q[2*i-1])/q[0];
    sum += qq*qq;
  }
  sum = sqrt(sum)/2.;

  printf(" --- total eror: %10.5f\n",sum);

}


//_____________________________________________________________________________
void get_pdf_err_wz(TString Lepton, double* Qw, double* Qz) {
  // Lepton = "e" or "mu", plot ratio of W/Z acceptances 

  const char module    [   ] = "TauMcAna";
  char       get_hist_c[200], file_w[200], file_z[200];

  const char* work_dir = gSystem->Getenv("SRT_LOCAL");
  sprintf(get_hist_c,"%s/plot/get_hist.C",work_dir);

  TH2F    *hw, *hz;
  TString hist_name;

  if (!gInterpreter->IsLoaded(get_hist_c)) gInterpreter->LoadMacro(get_hist_c);

  TString form_w;
  TString form_z;

  if (Lepton == "e") {
    form_w = "hist/dev_240/pythia_ana/wenu/wenu.pythia_ana.%i.1000000.root";
    form_z = "hist/dev_240/pythia_ana/zee/zee.pythia_ana.%i.1000000.root";
  }
  else if (Lepton == "e_3") {
    form_w = "hist/dev_240/pythia_ana/wenu/wenu.pythia_ana.%i.1000000.root";
    form_z = "hist/dev_240/pythia_ana/zee_3/zee.pythia_ana.%i.1000000.root";
  }
  else if (Lepton == "mu") {
    form_w = "hist/dev_240/pythia_ana/wmunu/wmunu.pythia_ana.%i.1000000.root";
    form_z = "hist/dev_240/pythia_ana/zmumu/zmumu.pythia_ana.%i.1000000.root";
  }
  else if (Lepton == "mu_d0") {
    form_w = "hist/dev_240/pythia_ana/wmunu/wmunu.pythia_ana.%i.1000000.root";
    form_z = "hist/dev_240/pythia_ana/zmumu_d0/zmumu.pythia_ana.%i.1000000.root";
  }

  for (int i=0; i<41; i++) { 
    sprintf(file_w,form_w.Data(),30200+i);
    sprintf(file_z,form_z.Data(),30200+i);
    
    if ((Lepton == "e") || (Lepton == "e_3")) {
//-----------------------------------------------------------------------------
// electron channel, do not cut on missing Et for Z's (use "nt_0" histogram)
//-----------------------------------------------------------------------------
      hw   = (TH2F*) get_hist(file_w,module,"genp/ele_nl_vs_nt_1");
      Qw[i] = hw->GetBinContent(2,2); 

      hz   = (TH2F*) get_hist(file_z,module,"genp/ele_nl_vs_nt_0");
      Qz[i] = hz->GetBinContent(3,3)+hz->GetBinContent(2,3);
    }
    else if ((Lepton == "mu") || (Lepton == "mu_d0")) {
//-----------------------------------------------------------------------------
// muon channel, do not cut on missing Et for Z's (use "nt_0" histogram)
//-----------------------------------------------------------------------------
      hw   = (TH2F*) get_hist(file_w,module,"genp/muo_nl_vs_nt_1");
      Qw[i] = hw->GetBinContent(2,2); 

      hz   = (TH2F*) get_hist(file_z,module,"genp/muo_nl_vs_nt_0");
      Qz[i] = hz->GetBinContent(3,3)+hz->GetBinContent(2,3);
    }
  }
}


//_____________________________________________________________________________
void get_pdf_err_w_tau_e(double* Qwtaunu, double* Qwenu) {
  // plot ratio of W->taunu/W->enu acceptances 

  const char module    [   ] = "TauMcAna";
  char       get_hist_c[200], file_wtaunu[200], file_wenu[200];

  const char* work_dir = gSystem->Getenv("SRT_LOCAL");
  sprintf(get_hist_c,"%s/plot/get_hist.C",work_dir);

  TH2F    *hw, *hz;
  TString hist_name;

  if (!gInterpreter->IsLoaded(get_hist_c)) gInterpreter->LoadMacro(get_hist_c);

  TString form_wtaunu;
  TString form_wenu;

  form_wtaunu = "hist/dev_240/pythia_ana/wtaunu/wtaunu.pythia_ana.%i.1000000.root";
  form_wenu   = "hist/dev_240/pythia_ana/wenu/wenu.pythia_ana.%i.1000000.root";

  for (int i=0; i<41; i++) { 
    sprintf(file_wtaunu,form_wtaunu.Data(),30200+i);
    sprintf(file_wenu  ,form_wenu.Data  (),30200+i);
    
    hwtaunu    = (TH2F*) get_hist(file_wtaunu,module,"genp/tau_nt_1");
    Qwtaunu[i] = hwtaunu->GetBinContent(2); 

    hwenu      = (TH2F*) get_hist(file_wenu,module,"genp/ele_nl_vs_nt_1");
    Qwenu[i]   = hwenu->GetBinContent(2,2);
  }
}


//_____________________________________________________________________________
void plot_pdf_err_wz(TString Lepton) {
  double qw[41], qz[41], q1, q2, qq, sum;

  get_pdf_err_wz(Lepton,qw,qz);
//-----------------------------------------------------------------------------
// now plot histogram with variations
//-----------------------------------------------------------------------------
  TH1F* h1 = new TH1F(Form("pdf_err_wz_%s",Lepton.Data()),
		      Form("ratio of W/Z acceptancies, channel %s",Lepton.Data()),
		      100,0.95,1.05);

  for (int i=1; i<=40; i++) {
    qq = (qw[i]/qz[i])/(qw[0]/qz[0]);
    h1->Fill(qq);
  }
  h1->Draw();

  sum = 0;
  for (int i=1; i<=20; i++) {
    q1 = (qw[2*i  ]/qz[2*i  ])/(qw[0]/qz[0]);
    q2 = (qw[2*i-1]/qz[2*i-1])/(qw[0]/qz[0]);
    qq = q1-q2;
    sum += qq*qq;
    printf(" i, qq^2 = %3i %10.6f \n",i,qq*qq);
  }
  sum = sqrt(sum)/2.;

  printf(" --- total eror: %10.5f\n",sum);

}


//_____________________________________________________________________________
void plot_pdf_err_w_tau_e() {
  // books gH1 and gH2

  double qtau[41], qe[41], q1, q2, qq, sum, dtau, de;

  get_pdf_err_w_tau_e(qtau,qe);
//-----------------------------------------------------------------------------
// now plot histogram with variations
//-----------------------------------------------------------------------------
  gH1 = new TH1F(Form("pdf_err_w_tau_e"),
		 Form("ratio of W->taunu/W->enu acceptancies"),
		 100,0.95,1.05);

  for (int i=1; i<=40; i++) {
    qq = (qtau[i]/qe[i])/(qtau[0]/qe[0]);
    gH1->Fill(qq);
  }


  gH2 = new TH2F(Form("pdf_corr_wtaunu_wenu"),
		 Form("Acceptance: W->enu Vs W->taunu"),
		 100,0.95,1.05,100,0.95,1.05);
  gH2->GetXaxis()->SetTitle("variation of W #rightarrow #tau #nu acceptance");
  gH2->GetYaxis()->SetTitle("variation of W #rightarrow e #nu acceptance");

  for (int i=1; i<=40; i++) {
    dtau = qtau[i]/qtau[0];
    de   = qe  [i]/qe  [0];
    gH2->Fill(dtau,de,1);
  }

  sum = 0;
  for (int i=1; i<=20; i++) {
    q1 = (qtau[2*i  ]/qe[2*i  ])/(qtau[0]/qe[0]);
    q2 = (qtau[2*i-1]/qe[2*i-1])/(qtau[0]/qe[0]);
    qq = q1-q2;
    sum += qq*qq;
  }
  sum = sqrt(sum)/2.;

  printf(" --- total eror: %10.5f\n",sum);

}


//_____________________________________________________________________________
void plot_pdf_corr_wz(TString Lepton) {
  double qw[41], qz[41], dw, dz;

  get_pdf_err_wz(Lepton,qw,qz);
//-----------------------------------------------------------------------------
// now plot histogram with variations
//-----------------------------------------------------------------------------
  TH2F* h2 = new TH2F(Form("pdf_corr_wz_%s",Lepton.Data()),
		      Form("dZ vs Dw for %s",Lepton.Data()),
		      100,0.95,1.05,100,0.95,1.05);
  h2->GetXaxis()->SetTitle(Form("variation of W #rightarrow %s #nu acceptance",
				Lepton.Data()));
  h2->GetYaxis()->SetTitle(Form("variation of Z #rightarrow %s %s acceptance",
			   Lepton.Data(),Lepton.Data()));

  for (int i=1; i<=40; i++) {
    dw = qw[i]/qw[0];
    dz = qz[i]/qz[0];
    h2->Fill(dw,dz,1);
  }
  h2->Draw();
}


//_____________________________________________________________________________
void plot_pdf_err_slide() {

  TCanvas* c = new_slide("c_pdf_err", "PDF errors",2,3);

  TPad* p1 = c->GetPrimitive("p1");

  p1->cd(1);
  plot_pdf_err("wenu");

  p1->cd(2);
  plot_pdf_err("wmunu");

  p1->cd(3);
  plot_pdf_err("wtaunu");

}