plot_estimate_energy.C

Go to the documentation of this file.

//----------------------------------------
// Making Energy Estimator Plots for 
// Inclusive numu CC xsection Analysis.
//
// Author : Biswaranjan Behera
// email  : bbehera@fnal.gov
//----------------------------------------

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SystShifts.h"

#include <vector>
#include <string>
#include <memory>

#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TProfile.h"

#include "Utilities/rootlogon.C"

using namespace ana;

enum MuonFitMode
{
  kMuonFitAct,
  kMuonFitActAndCat,
  kMuonFitCat,
};

TF1 * linposfit = new TF1("linposfit", "[0]+[1]*x", 0., 15.);

void MuonFit(Spectrum *plot1, string zaxisname, 
             string label,
             MuonFitMode fitmode,
             string outDir=""){
  
  double pot = 11E20;
  TH2* hist  = plot1->ToTH2(pot);
  hist->SetTitle("");
  switch (fitmode)
  {
    case kMuonFitAct:
      hist->SetTitle("Active Region Only");
      hist->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
      hist->GetYaxis()->SetTitle("True Muon Energy (GeV)");
      break;
    case kMuonFitActAndCat:
    if (fitmode == kMuonFitActAndCat){
      hist->SetTitle("Muons Entering Catcher");
      hist->GetXaxis()->SetTitle("Reco Muon Track Length in Active (m)");
      hist->GetYaxis()->SetTitle("True Active Region Muon Energy (GeV)");  
      hist->GetYaxis()->SetTitleSize(0.048);
      break;
    case kMuonFitCat:
      hist->SetTitle("Muons Entering Catcher");
      hist->GetXaxis()->SetTitle("Reco Muon Track Length in Catcher (m)");
      hist->GetYaxis()->SetTitle("True Catcher Region Muon Energy (GeV)");
      break;
    default:
      throw std::runtime_error("Invalid fit mode");
    }
  }
  
  hist->GetZaxis()->SetTitle(zaxisname.c_str());
  static const int nn = 150; //Number of points on graph - should be equal to or smaller than number of bins in x
  float xxx [nn];
  float exxl[nn];
  float exxh[nn];
  float yyx [nn];
  float eyxl[nn];
  float eyxh[nn];
    
  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
    // make projection in Y view
    TH1D* py       = hist->ProjectionY("py",iw,iw,"");
    int maxbin     = py->GetMaximumBin();        //Bin in slice with most entries
    double maxcont = py->GetBinContent(maxbin);  //Height of the bin with most entries
    
    if (maxcont == 0){continue; }    //If no entries in the slice, don't make a graph point
    
    double totalEntries = py->GetEntries();
    if (totalEntries < 30){continue; }    //If not many entries in the slice, don't make a graph point
    
    double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
    
    //Find the location of width at half max
    int yyhighbin = maxbin+1;
    int yylowbin  = maxbin-1;
    
    // find upper bin of half max
    for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
      yyhighbin = ihigh;
      if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
    }
    
    // find lower bin of half max
    for(int ilow = maxbin-1; ilow > 0; --ilow){
      yylowbin = ilow;
      if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
    }
    
    float multiplier  = 1.5; //This defines how many times wider than width at half max I want to fit over
    int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
    int veryLowBin  = maxbin - multiplier*(maxbin-yylowbin);
    
    double lowSide  = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
    double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
    
    xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
    
    TF1 *f1 = new TF1("customGaus","gaus",0,5);
    py->Fit("customGaus","OQ","",lowSide,highSide); //O and Q mean do not plot result, quiet mode
    double mean = f1->GetParameter(1);
    double error = f1->GetParError(1);
    double chi = f1->GetChisquare();
    double ndf = f1->GetNDF();
    
    //    if (ndf > 0) {chiValues->Fill(chi/ndf);}
    //else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
    yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
    //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
    exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
    exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
    //Setting the error on the mean reported by the fit as error in y
    eyxl[iw-1] = error;
    eyxh[iw-1] = error;
    
  }//End of loop over bins in 2D histogram
  CenterTitles(hist);
  // TProfile * profile = hist->ProfileX();

  //->GetZaxis()->CenterTitle();
  //Making graph from points
  TGraphAsymmErrors* customProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
  // TH1D* customProfile = profile->ProjectionX();
  customProfile->GetXaxis()->SetNoExponent(kTRUE);
  customProfile->GetXaxis()->CenterTitle();
  customProfile->GetYaxis()->CenterTitle();
  switch (fitmode){
    case kMuonFitAct:
      customProfile->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
      customProfile->GetYaxis()->SetTitle("True Muon Energy (GeV)");
      customProfile->GetXaxis()->SetLimits(0,15);
      customProfile->GetHistogram()->SetMaximum(5);
      customProfile->GetHistogram()->SetMinimum(0);  
      // customProfile->SetMaximum(5);
      // customProfile->SetMinimum(0);  
      break;
    case kMuonFitActAndCat:
      customProfile->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
      customProfile->GetYaxis()->SetTitle("True Muon E - True Muon E in Catcher (GeV)");
      customProfile->GetXaxis()->SetLimits(0,15);  
      customProfile->GetHistogram()->SetMaximum(5);
      customProfile->GetHistogram()->SetMinimum(0);
      // customProfile->SetMaximum(5);
      // customProfile->SetMinimum(0);
      customProfile->GetYaxis()->SetTitleSize(0.048);
      break;
    case kMuonFitCat:
      customProfile->GetXaxis()->SetTitle("Reco Muon Track Length in Catcher (m)");
      customProfile->GetYaxis()->SetTitle("True Muon Energy in Catcher (GeV)");
      customProfile->GetXaxis()->SetLimits(0,3);
      customProfile->GetHistogram()->SetMaximum(3);
      customProfile->GetHistogram()->SetMinimum(0); 
      // customProfile->SetMaximum(3);
      // customProfile->SetMinimum(0); 
      break;
    default:
      throw std::runtime_error("Invalid fit mode");
  }
  
  //CenterTitles(customProfile);
  customProfile->SetTitle("");
  
  //Overlays 2d plot and graph
  TCanvas* c1 = new TCanvas("c1","c1");
  customProfile->SetMarkerColor(18);
  customProfile->SetLineColor(18);
  hist->Draw("colz");
  customProfile->Draw("P");
  // customProfile->Draw("e same");
  Simulation();
  c1->SetRightMargin(0.15);
  c1->SetLeftMargin(0.15);

  if (!outDir.empty()){
    c1->SaveAs((outDir+"/muon_" + label + ".png").c_str());
    c1->SaveAs((outDir+"/muon_" + label + ".pdf").c_str());
  }
  
  //Overlays 2d loqz plot and graph
  TCanvas* c2 = new TCanvas("c2","c2");
  c2->SetLogz();
  customProfile->SetMarkerColor(13);
  customProfile->SetLineColor(13);
  hist->Draw("colz");
  customProfile->Draw("P");
  // customProfile->Draw("e same");
  Simulation();
  c2->SetRightMargin(0.15);
  c2->SetLeftMargin(0.15);
  if (!outDir.empty()){
    c2->SaveAs((outDir+"/muon_logz" + label + ".png").c_str());
    c2->SaveAs((outDir+"/muon_logz" + label + ".pdf").c_str());
  }

  //Just the graph
  TCanvas* c3 = new TCanvas("c3","c3");
  customProfile->SetMarkerColor(1);
  customProfile->SetLineColor(1);
  customProfile->Draw("AP");
  // customProfile->Draw("e same");
  c3->SetRightMargin(0.15);
  switch (fitmode){
    case kMuonFitAct:
      customProfile->Fit("pol3", "","",0.2, 12.0);
      break;

    case kMuonFitActAndCat:
      customProfile->Fit("pol3", "","",3.75, 12.2);
      break;

    case kMuonFitCat:
      customProfile->Fit("pol1", "","",0.02, 2.3);
      break;

    default:
      throw std::runtime_error("Invalid fit mode");
  }
    
  c3->SetLeftMargin(0.15);

  if (!outDir.empty()){
    c3->SaveAs((outDir+ "/fit_" + label + ".png").c_str());
    c3->SaveAs((outDir+ "/fit_" + label + ".pdf").c_str());
  }
}

void plot_estimate_energy(string infile, string outputDir = ".")
{
  TFile *fin = TFile::Open(infile.c_str(),"READ");
  //------------------------------------------
  //            Accessing Files Here
  //------------------------------------------
  
  Spectrum* MuonE_hist_act        = Spectrum::LoadFrom(fin, "MuonE_hist_act"      ).release();
  Spectrum* MuonE_hist_actAndCat  = Spectrum::LoadFrom(fin, "MuonE_hist_actAndCat").release();
  Spectrum* MuonE_hist_cat        = Spectrum::LoadFrom(fin, "MuonE_hist_cat"      ).release();
  
  double pot = 11E20;

  //------------------------------------------------------
  // Plots for Energy Estimator
  //------------------------------------------------------
 
  MuonFit(MuonE_hist_act, Form("Events/(%2.2E POT)",pot),
          "act", kMuonFitAct,outputDir);

  MuonFit(MuonE_hist_actAndCat, Form("Events/(%2.2E POT)",pot),
          "act_cat", kMuonFitActAndCat,outputDir);

  MuonFit(MuonE_hist_cat, Form("Events/(%2.2E POT)",pot),
          "cat", kMuonFitCat,outputDir);


  // HadEFit(HadE, Form("Events/(%2.2E POT)",pot),
  //      "hade","hade_Logz",
  //      "hade_fit", true);

}// end of void