plot_bbf_overlay.C

plot_bbf_overlay.C

{

  TString lax_version = "1.2.3";
  
  const int nHists = 4;
  TString hname = "hbkg_py0_1.00_rf0_1.00";
  TString hname_append[nHists] = {"_data", "", "_debug1", "_debug2"};
  TString htitle[nHists] = {"Data", "MC (new, ef8ae8a)", "MC (debug1)", "MC (debug2)"};
  enum henum {data, mc, mc_old};
  int icolor[nHists] = {kBlue, kBlack, kRed, kGreen+2};
  int istyle[nHists] = {1, 1, 1, 1};
  
  TH2F *hist[nHists];
  TH1D *hist_proj[nHists];
  
  TString infile_name[nHists] = {
    "data/Rn220_SR1_pax6.8.0_hax2.1.1_lax1.2.3_cs1LT200.root",
    "ERBackgroundModel_SR1_From_Combined180104Fit.root",
    //"ERBackgroundModel_SR1_From_Combined180104Fit_2e672f0.root"
    "ForDebug1.root",
    "ForDebug2.root",
  };

  // Load MC
  for (int ihist=1; ihist<nHists; ihist++) {
    TFile *infile = new TFile(infile_name[ihist]);

    if (ihist>=mc_old) hname = "hmc";
    
    hist[ihist] = (TH2F*)infile->Get(hname);
    hist[ihist]->SetName(hname + hname_append[ihist]);
    hist[ihist]->SetTitle(";cs1 (pe); log10(cs2_bottom)");
    cout << "Loaded: " << infile->GetName() << " " << hist[ihist]->Integral() << endl;
  }

  // Load data
  TFile *infile = new TFile(infile_name[data]);
  TTree *tree = (TTree*)infile->Get("tree");
  hist[data] = (TH2F*)hist[mc]->Clone();
  hist[data]->SetName(hname + hname_append[data]);
  hist[data]->Reset();
  tree->Project(hist[data]->GetName(), "log10(cs2_bottom):cs1");

  // Project to 1D
  double low_cs1 = 40;
  double hi_cs1 = 60;
  for (int ihist=0; ihist<nHists; ihist++) {

    int lowbin = hist[ihist]->GetXaxis()->FindBin(low_cs1);
    int hibin = hist[ihist]->GetXaxis()->FindBin(hi_cs1) - 1;

    // Normalize MC to data
    if (ihist != data) hist[ihist]->Scale(hist[data]->Integral()/hist[ihist]->Integral());

    TString proj_opt = "";
    if (ihist == data) proj_opt = "e";
    
    hist_proj[ihist] = hist[ihist]->ProjectionY("_py", lowbin, hibin, proj_opt);
    hist_proj[ihist]->SetLineWidth(2);
    hist_proj[ihist]->SetTitle(Form("Rn220 ER Calibration: %.0f < cs1 < %.0f pe", low_cs1, hi_cs1));
    hist_proj[ihist]->GetYaxis()->SetTitle("Number of events");

  }

  // Draw
  TLegend *leg = new TLegend(0.63, 0.55, 0.92, 0.85);
  leg->SetFillColor(0);

  bool isDrawn = 0;
  for (int ihist=nHists-1; ihist>=0; ihist--) {

    if (ihist==data) {
      hist_proj[ihist]->SetMarkerStyle(8);
      hist_proj[ihist]->SetMarkerSize(0.8);
    }
    
    hist_proj[ihist]->SetMarkerColor(icolor[ihist]);
    hist_proj[ihist]->SetLineColor(icolor[ihist]);
    hist_proj[ihist]->SetLineStyle(istyle[ihist]);
    
    if (isDrawn==0) {
      hist_proj[ihist]->Draw();
      hist_proj[ihist]->GetYaxis()->SetRangeUser(0, 220);
      hist_proj[ihist]->GetXaxis()->SetRangeUser(3, 3.8);
      isDrawn = 1;
    }
    
    else
      hist_proj[ihist]->Draw("same");

    TString leg_opt;
    if (ihist==data) leg_opt = "p";
    else leg_opt = "l";
    leg->AddEntry(hist_proj[ihist], htitle[ihist], leg_opt);
  }

  leg->Draw();
  
  TString cname = Form("images/ER_cs1_%d_%d.png", (int)low_cs1, (int)hi_cs1);
  c1->Print(cname);
}