Plotting_Data2DataComp.C

Go to the documentation of this file.

/*******************************************************************
 * A SCRIPT TO READ THE OUTPUT ROOT FILE OF Analyse_Data2DataComp.C,
 * WILL PLOT SPECTRA FROM BOTH DATA SETS AND THE RATIO OF THE TWO.
 * OUTPUT CSV FILE CREATED WITH KEY QUANTITIES SUCH AS THE FRACTIONAL
 * CHANGES IN THE MEANS ETC.
 * USAGE:
 * Plotting_Data2DataComp.C <Input root file containing spectra>
 *                          <Directory to output pdfs>
 *                          <String to append outfile names>
 *                          <int color for first spectra>
 *                          <int color for second spectra>
 ********************************************************************/

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Analysis/Plots.h"

#include <fstream>
#include "TCanvas.h"
#include "TFile.h"  
#include "TGraph.h" 
#include "TH1.h"    
#include "TH2.h"    
#include "TH2D.h"    
#include "TMath.h"  
#include "TGaxis.h" 
#include "TMultiGraph.h"
#include "TLegend.h"    
#include "TLatex.h"     
#include "TStyle.h"     
#include "THStack.h"    
#include "TPaveText.h"
#include "TLine.h"


using namespace ana;

double getFractionalError(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
{
  double fractionalError = 0.;

  double numerator      = mean1 - mean0; 
  double denominator    = mean0;  
  double numeratorVar   = meanErr0*meanErr0 + meanErr1*meanErr1;
  double denominatorVar = meanErr0*meanErr0;

  std::cout << numerator << " " << denominator << " " << numeratorVar << " " << denominatorVar << std::endl;

  return std::sqrt((numerator/denominator)*(numerator/denominator)*(numeratorVar/(numerator*numerator)+denominatorVar/(denominator*denominator)))*100.;
}

std::string is2SigmaBandsTouching(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
{
  std::pair<double,double> largeMean;
  std::pair<double,double> smallMean;

  if(mean0 > mean1)
  {
    largeMean = {mean0, meanErr0};    
    smallMean = {mean1, meanErr1};
  }
  else
  {
    largeMean = {mean1, meanErr1};    
    smallMean = {mean0, meanErr0};
  }

  double largeMeanShifted = largeMean.first - 2.*largeMean.second;
  double smallMeanShifted = smallMean.first + 2.*smallMean.second;

  return (largeMeanShifted <= smallMeanShifted ? "True" : "False");
}

void Plotting_Data2DataComp(std::string s_InFile = "", std::string s_OutDir = "", std::string s_OutName = "",
                                                 unsigned int color1 = 1, unsigned int color2 = 2) {

  std::cout << "\n\n\nMAKING THE PLOTS...\n\n\n";

  std::ofstream f_MeansFile;
  f_MeansFile.open(s_OutDir+"meansFile_"+s_OutName+".csv");

  TFile f_In(s_InFile.c_str());
  char s_DirName [256];
  char s_FileName[256];
  unsigned int nQuant = 4;
  unsigned int nSpec  = 41;

  for(unsigned int S = 0; S <= nSpec; ++S) {
    for(unsigned int Q = 0; Q <= nQuant; ++Q) {

      sprintf(s_DirName, "dir_Var%d_Q%d_DataSet0",S,Q);
      std::unique_ptr<Spectrum> spec0 = Spectrum::LoadFrom(f_In, s_DirName);
      sprintf(s_DirName, "dir_Var%d_Q%d_DataSet1",S,Q);
      std::unique_ptr<Spectrum> spec1 = Spectrum::LoadFrom(f_In, s_DirName);
      
      TCanvas *c  = new TCanvas("c", "c", 1000, 800);      
      TPad    *p1 = new TPad("p1", "p1", 0, 0, 1, 1);
      TPad    *p2 = new TPad("p2", "p2", 0, 0, 1, 1);
      SplitCanvas(0.3,p1,p2);

      p1->cd();
      //p1->SetLogy();
      gStyle->SetOptStat(0);

      TH1 *h0 = spec0->ToTH1(spec1->POT());
      TH1 *h1 = spec1->ToTH1(spec1->POT());

      h0->SetLineWidth(2);
      h0->SetLineColor(color1);

      h1->SetLineWidth(2);
      h1->SetLineColor(color2);

      h0->GetXaxis()->SetTitle("");
      h0->GetXaxis()->SetLabelSize(0);
      h1->GetXaxis()->SetTitle("");
      h1->GetXaxis()->SetLabelSize(0);

      h0->Draw();
      h1->Draw("SAME");
      
      p2->cd();
      Ratio rat(*spec1, *spec0);
      std::cout << "DIVINING " << spec1->GetLabels()[0] << " BY " << spec0->GetLabels()[0] << std::endl; 

      TH1* hRat = rat.ToTH1();
      hRat->SetAxisRange(0.9,1.1,"Y");
      hRat->Draw();

      double x[2] = {hRat->GetXaxis()->GetXmin(), hRat->GetXaxis()->GetXmax()};
      double y[2] = {1.0, 1.0};
      TGraph* l_Xaxis = new TGraph(2, x, y);
      l_Xaxis->SetLineColor(kGray);
      l_Xaxis->SetLineWidth(3);
      l_Xaxis->SetLineStyle(2);
      l_Xaxis->Draw("L SAME");

      sprintf (s_FileName, "Var%d_Q%d_%s.pdf", S, Q+1, s_OutName.c_str());
      c->Print((s_OutDir+(std::string)s_FileName).c_str());

      std::cout << "\n\nMeans for:\t" << hRat->GetXaxis()->GetTitle() << " for quant "      << Q+1 << "\n"
                << "old data:   \t"   << h0->GetMean() << " +/- "     << h0->GetMeanError() <<        "\n"
                << "keepup data:\t"   << h1->GetMean() << " +/- "     << h1->GetMeanError() <<        "\n\n";

      f_MeansFile << s_FileName << "=" << hRat->GetXaxis()->GetTitle() << " for quant " << Q+1 << "=" 
                  << h0->GetMean() << "=" << h0->GetMeanError() << "="
                  << h1->GetMean() << "=" << h1->GetMeanError() << "="
                  << (h1->GetMean()-h0->GetMean())/h0->GetMean()*100. << "=" 
                  << getFractionalError   (h0->GetMean(), h0->GetMeanError(), h1->GetMean(), h1->GetMeanError()) << "=" 
                  << is2SigmaBandsTouching(h0->GetMean(), h0->GetMeanError(), h1->GetMean(), h1->GetMeanError()) << std::endl;

      delete l_Xaxis;
      delete c;
    }
  }

  f_MeansFile.close();

  std::cout << "\n\n\nFINISHED MAKING PLOTS!\n\n\n";

}