absCal.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////
//                                                            //
//    Absolute calibratrion parameters                        //
//                                                            //
//    Created by:  Luke Vinton < July/2015                    //
//    Modified by: Diana Patricia Mendez > July 2015 ...      //
//                                                            //
//    Output:   cMCdataND.pdf     - Corrected response/cm     //
//                                = (PECorr/cm) = MEU         //
//              cMCND_truth.pdf   - Simulated dE/dx (MeV/cm)  //
//              cMCdataND_cal.pdf - Calibrated dE/dx (MeV/cm) //
//                                                            //
//  Calibrated dE/dx is multiplied by the conversion factor   //
//                                                            //
////////////////////////////////////////////////////////////////

#include <iostream>
#include <vector>

#include <TFile.h>
#include <TPad.h>
#include <TH1.h>
#include <TH2.h>
#include <TH2F.h>
#include <TH2D.h>
#include <TProfile2D.h>
#include <TH1D.h>
#include <TH1F.h>
#include <TH3D.h>
#include <TCanvas.h>
#include <TVector.h>
#include <TColor.h>
#include <TLegend.h>
#include <TF1.h>
#include <TLatex.h>
#include <TMath.h>
#include <TStyle.h>
#include <TROOT.h>
#include <TPad.h>
#include <TFile.h>
#include <TTree.h>
#include <TSystem.h>
#include <TColor.h>
#include "TChain.h"

#include "Utilities/rootlogon.C"

// canvas size
#define c1_x 700
#define c1_y 500 

#define c2_x 1000
#define c2_y 500




double Landau(TH1D* h, double& err);
double Gaus(TH1D* h, double& err, bool isTruth);
double median(TH1D* h);
std::vector<double> quantiles( TH1D* h );

double errAoverB(double& err, double A, double B, double Aerr, double Berr);
double GetHistMax(TH1* h0, TH1* h1, TH1* h2, TH1* h3);


int main(){

  //This constants refer to the MEU value for mc and data
  //This should be changed to get the right calibrated plot
  double const mccal = 0.0451381;
  double const datacal = 0.0442841;

  gROOT->ForceStyle();

  //Using TChain avoids using merge and then having an only input file
  //which sometimes do not work due to size and nunber or files to merge
  //specially for FD
  std::cout << "Making MC chain ..." << std::endl;
  TChain* tNDMC = new TChain("muondedxana/fTree");
  tNDMC->Add("mc/*.root");

  double PECorr_MC, path_MC, trueE_MC, x_MC;
  double planes_MC, xplanes_MC, yplanes_MC;
  int cont_MC, nhits_MC,hitId_MC;
  
  //Declaring branches saves time compared to the Draw function
  std::cout << "Setting MC branches ..." <<std::endl;
  tNDMC->SetBranchAddress("PECorr", &PECorr_MC);
  tNDMC->SetBranchAddress("path"  , &path_MC);
  tNDMC->SetBranchAddress("trueE" , &trueE_MC);
  tNDMC->SetBranchAddress("x"     , &x_MC);


  std::cout << "Declaring MC histograms ..." <<std::endl;
  TH1D* hMCND_PECorr = new TH1D("hMCND_PECorr", "hMCND_PECorr", 1000, 0, 500);
  TH1D* hMCND_trueE  = new TH1D("hMCND_trueE", "hMCND_trueE", 1000, 0, 50);
  TH1D* hMCND_cal    = new TH1D("hMCND_cal", "hMCND_cal", 60, 0, 6);

  std::cout << "Filling MC histos ..." <<std::endl;
  unsigned int nEntries_MC = tNDMC->GetEntries();
  for(unsigned int i=0; i<nEntries_MC; i++){
    tNDMC->GetEntry(i);

    if( x_MC>100 && x_MC<200){
      if(PECorr_MC>0){
        hMCND_PECorr ->Fill(PECorr_MC/path_MC);
        hMCND_cal ->Fill(mccal*PECorr_MC/path_MC);
      }
      if(trueE_MC>0){
        hMCND_trueE ->Fill(trueE_MC/path_MC);
      }
    }
    
  }
  


  std::cout << "Making data chain ..." <<std::endl;
  TChain* tNDdata = new TChain("muondedxana/fTree");
  tNDdata->Add("data/*.root");                                                 

  double PECorr_data, path_data, trueE_data, x_data;
  double planes_data, xplanes_data, yplanes_data;
  int cont_data, nhits_data,hitId_data;

  std::cout << "Setting data branches ..." <<std::endl;
  tNDdata->SetBranchAddress("PECorr", &PECorr_data);
  tNDdata->SetBranchAddress("path"  , &path_data);
  tNDdata->SetBranchAddress("trueE" , &trueE_data);
  tNDdata->SetBranchAddress("x"     , &x_data);

  std::cout << "Declaring data histos ..." <<std::endl;

  TH2D* hdataND_PECorrX = new TH2D("hdataND_PECorrX", "hdataND_PECorrX", 160, -800, 800, 30, 0, 60);
  TH1D* hdataND_PECorr  = new TH1D("hdataND_PECorr","hdataND_PECorr", 1000, 0, 500);
  TH1D* hdataND_cal     = new TH1D("hdataND_cal","hdataND_cal", 60, 0, 6);

  std::cout << "Filling data histos ..." <<std::endl;

  unsigned int nEntries_data = tNDdata->GetEntries();
  for(unsigned int i=0; i<nEntries_data; i++){
    tNDdata->GetEntry(i);

    //    if(hitId_data!=2) continue;
    // if(!PassActivityTrigger(planes_data, xplanes_data, yplanes_data, nhits_data, cont_data)) continue;
    
    if(x_data>100 && x_data<200)
      if(PECorr_data>0){
        hdataND_PECorrX ->Fill(x_data,PECorr_data/path_data);
        hdataND_PECorr  ->Fill(PECorr_data/path_data);
        hdataND_cal     ->Fill(datacal*PECorr_data/path_data);

      }
    
  }
  


    
  //landau
  TH1D* hdataND_PECorr_landau = (TH1D*)hdataND_PECorr->Clone("hdataND_PECorr_landau");
  TH1D* hMCND_PECorr_landau   = (TH1D*)hMCND_PECorr  ->Clone("hMCND_PECorr_landau");
  TH1D* hMCND_trueE_landau    = (TH1D*)hMCND_trueE   ->Clone("hMCND_trueE_landau");
  //gaus
  TH1D* hdataND_PECorr_gaus   = (TH1D*)hdataND_PECorr->Clone("hdataND_PECorr_gaus");
  TH1D* hMCND_PECorr_gaus     = (TH1D*)hMCND_PECorr  ->Clone("hMCND_PECorr_gaus");
  TH1D* hMCND_trueE_gaus      = (TH1D*)hMCND_trueE   ->Clone("hMCND_trueE_gaus");

  //set data lines to black:
  hdataND_PECorr        -> SetLineColor(1);
  hdataND_PECorr_landau -> SetLineColor(1);
  hdataND_PECorr_gaus   -> SetLineColor(1);

  //set histo range
  hdataND_PECorr -> GetXaxis() -> SetRangeUser(0,100);
  hMCND_PECorr   -> GetXaxis() -> SetRangeUser(0,100);
  hMCND_trueE    -> GetXaxis() -> SetRangeUser(0,6);
  hdataND_cal -> GetXaxis() -> SetRangeUser(0,6);
  hMCND_cal   -> GetXaxis() -> SetRangeUser(0,6);


  double mpvErr_dataND, mpvErr_MCND, mpvErr_trueE_MCND;
  double mpv_dataND = Landau(hdataND_PECorr_landau, mpvErr_dataND);  
  double mpv_MCND = Landau(hMCND_PECorr_landau, mpvErr_MCND);  
  double mpv_trueE_MCND = Landau(hMCND_trueE_landau, mpvErr_trueE_MCND);

  double gausErr_dataND, gausErr_MCND, gausErr_trueE_MCND;
  double gausMean_dataND = Gaus(hdataND_PECorr_gaus, gausErr_dataND, false);  
  double gausMean_MCND = Gaus(hMCND_PECorr_gaus, gausErr_MCND, false);  
  double gausMean_trueE_MCND = Gaus(hMCND_trueE_gaus, gausErr_trueE_MCND, true);

  double mean_dataND = hdataND_PECorr->GetMean(1);
  double mean_MCND = hMCND_PECorr->GetMean(1);
  double mean_trueE_MCND = hMCND_trueE->GetMean(1);
  double meanErr_dataND = hdataND_PECorr->GetRMS()/(pow(hdataND_PECorr->GetEntries(), 0.5));
  double meanErr_MCND = hMCND_PECorr->GetRMS()/(pow(hMCND_PECorr->GetEntries(), 0.5));
  double meanErr_trueE_MCND = hMCND_trueE->GetRMS()/(pow(hMCND_trueE->GetEntries(), 0.5));

  double Median_dataND;
  double pointFive = 0.5;
  hdataND_PECorr->GetQuantiles(1,&Median_dataND,&pointFive);
  double Median_MCND;
  hMCND_PECorr->GetQuantiles(1,&Median_MCND,&pointFive);
  double Median_trueE_MCND;
  hMCND_trueE->GetQuantiles(1,&Median_trueE_MCND,&pointFive);


  // print MEU values from all techniques along with errors
  std::cout<< " ----------- Printing MEU values: --------" <<std::endl;
  std::cout<< "*** Mean: *** " <<std::endl;
  std::cout<< "dataND: " << mean_dataND << " +/- " << meanErr_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << mean_MCND << " +/- " << meanErr_MCND << " (stat.)" <<std::endl;
  std::cout<< "trueE_MCND: " << mean_trueE_MCND << " +/- " << meanErr_trueE_MCND << " (stat.)" <<std::endl;

  std::cout<< "*** Landau: *** " <<std::endl;
  std::cout<< "dataND: " << mpv_dataND << " +/- " << mpvErr_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << mpv_MCND << " +/- " << mpvErr_MCND << " (stat.)" <<std::endl;
  std::cout<< "trueE_MCND: " << mpv_trueE_MCND << " +/- " << mpvErr_trueE_MCND << " (stat.)" <<std::endl;

  std::cout<< "*** Gaus: *** " <<std::endl;
  std::cout<< "dataND: " << gausMean_dataND << " +/- " << gausErr_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << gausMean_MCND << " +/- " << gausErr_MCND << " (stat.)" <<std::endl;
  std::cout<< "trueE_MCND: " << gausMean_trueE_MCND << " +/- " << gausErr_trueE_MCND << " (stat.)" <<std::endl;



  //abs. scales with diff. techniques
  //landau:
  double abs_landau_dataND = mpv_trueE_MCND/mpv_dataND;
  double abs_landau_MCND = mpv_trueE_MCND/mpv_MCND;
  double absErr_landau_dataND, absErr_landau_MCND;
  errAoverB(absErr_landau_dataND, mpv_trueE_MCND, mpv_dataND, mpvErr_trueE_MCND, mpvErr_dataND);
  errAoverB(absErr_landau_MCND, mpv_trueE_MCND, mpv_MCND, mpvErr_trueE_MCND, mpvErr_MCND);
  //gaus:
  double abs_gaus_dataND = gausMean_trueE_MCND/gausMean_dataND;
  double abs_gaus_MCND = gausMean_trueE_MCND/gausMean_MCND;
  double absErr_gaus_dataND, absErr_gaus_MCND;
  errAoverB(absErr_gaus_dataND, gausMean_trueE_MCND, gausMean_dataND, gausErr_trueE_MCND, gausErr_dataND);
  errAoverB(absErr_gaus_MCND, gausMean_trueE_MCND, gausMean_MCND, gausErr_trueE_MCND, gausErr_MCND);
  //mean:
  double abs_mean_dataND = mean_trueE_MCND/mean_dataND;
  double abs_mean_MCND = mean_trueE_MCND/mean_MCND;
  double absErr_mean_dataND, absErr_mean_MCND;
  errAoverB(absErr_mean_dataND, mean_trueE_MCND, mean_dataND, meanErr_trueE_MCND, meanErr_dataND);
  errAoverB(absErr_mean_MCND, mean_trueE_MCND, mean_MCND, meanErr_trueE_MCND, meanErr_MCND);
  //median: 
  // Need to find technique to calculate error on the median...
  double abs_median_dataND = Median_trueE_MCND/Median_dataND;
  double abs_median_MCND = Median_trueE_MCND/Median_MCND;

  // Print Absolute scale values along with statistical errors
  std::cout<< "------------- Printing Abs. scale values: ----------- " <<std::endl;
  std::cout<< "*** Mean: *** " <<std::endl;
  std::cout<< "dataND: " << abs_mean_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << abs_mean_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;
  std::cout<< "*** Median: *** " <<std::endl;
  std::cout<< "dataND: " << abs_median_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << abs_median_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;

  std::cout<< "*** Landau: *** " <<std::endl;
  std::cout<< "dataND: " << abs_landau_dataND << " +/- " << absErr_landau_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << abs_landau_MCND << " +/- " << absErr_landau_MCND << " (stat.)" <<std::endl;
  std::cout<< "*** Gaus: *** " <<std::endl;
  std::cout<< "dataND: " << abs_gaus_dataND << " +/- " << absErr_gaus_dataND << " (stat.)" <<std::endl;
  std::cout<< "MCND: " << abs_gaus_MCND << " +/- " << absErr_gaus_MCND << " (stat.)" <<std::endl;




  std::ofstream values("absCal.txt");

  //printing values in created file
  values<< " ----------- Printing MEU values: --------" <<std::endl;
  values<< "*** Mean: *** " <<std::endl;
  values<< "dataND: " << mean_dataND << " +/- " << meanErr_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << mean_MCND << " +/- " << meanErr_MCND << " (stat.)" <<std::endl;
  values<< "trueE_MCND: " << mean_trueE_MCND << " +/- " << meanErr_trueE_MCND << " (stat.)" <<std::endl;

  values<< "*** Landau: *** " <<std::endl;
  values<< "dataND: " << mpv_dataND << " +/- " << mpvErr_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << mpv_MCND << " +/- " << mpvErr_MCND << " (stat.)" <<std::endl;
  values<< "trueE_MCND: " << mpv_trueE_MCND << " +/- " << mpvErr_trueE_MCND << " (stat.)" <<std::endl;

  values<< "*** Gaus: *** " <<std::endl;
  values<< "dataND: " << gausMean_dataND << " +/- " << gausErr_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << gausMean_MCND << " +/- " << gausErr_MCND << " (stat.)" <<std::endl;
  values<< "trueE_MCND: " << gausMean_trueE_MCND << " +/- " << gausErr_trueE_MCND << " (stat.)" <<std::endl;

  values<< "------------- Printing Abs. scale values: ----------- " <<std::endl;
  values<< "*** Mean: *** " <<std::endl;
  values<< "dataND: " << abs_mean_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << abs_mean_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;
  values<< "*** Median: *** " <<std::endl;
  values<< "dataND: " << abs_median_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << abs_median_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;

  values<< "*** Landau: *** " <<std::endl;
  values<< "dataND: " << abs_landau_dataND << " +/- " << absErr_landau_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << abs_landau_MCND << " +/- " << absErr_landau_MCND << " (stat.)" <<std::endl;
  values<< "*** Gaus: *** " <<std::endl;
  values<< "dataND: " << abs_gaus_dataND << " +/- " << absErr_gaus_dataND << " (stat.)" <<std::endl;
  values<< "MCND: " << abs_gaus_MCND << " +/- " << absErr_gaus_MCND << " (stat.)" <<std::endl;

  values.close();


  ////////////     plotting     ///////////

  TCanvas *cDataND = new TCanvas("cDataND","cDataND",c1_x,c1_y);
  cDataND->cd();
  hdataND_PECorr->SetTitle("");
  hdataND_PECorr->GetXaxis()->SetTitle("Calibrated detector response/cm");
  hdataND_PECorr->GetYaxis()->SetTitle("Hits in track window");
  CenterTitles(hdataND_PECorr);
  hdataND_PECorr->Draw();
  Preliminary();

  TCanvas *cDataND_landau = new TCanvas("cDataND_landau","cDataND_landau", c1_x, c1_y);
  cDataND_landau->cd();
  hdataND_PECorr_landau->GetXaxis()->SetTitle("Calibrated detector response/cm");
  hdataND_PECorr_landau->SetTitle("");
  CenterTitles(hdataND_PECorr_landau);
  hdataND_PECorr_landau->Draw();
  Preliminary();

  
  TCanvas *cDataND_gaus = new TCanvas("cDataND_gaus","cDataND_gaus",c1_x ,c1_y);
  cDataND_gaus->cd();
  hdataND_PECorr_gaus->GetXaxis()->SetTitle("Calibrated detector response/cm");
  hdataND_PECorr_gaus->SetTitle("");
  CenterTitles(hdataND_PECorr_gaus);
  hdataND_PECorr_gaus->Draw();
  Preliminary();


  TCanvas *cMCND = new TCanvas("cMCND","cMCND", c1_x, c1_y);
  cMCND -> cd();
  hMCND_PECorr->SetTitle("");
  hMCND_PECorr->GetXaxis()->SetTitle("dE/dx (MeV/cm)");
  hMCND_PECorr->GetYaxis()->SetTitle("Hits in track window");
  hMCND_PECorr->SetLineColor(2);
  CenterTitles(hMCND_PECorr);
  hMCND_PECorr->Draw();
  Simulation();

  cMCND->Print("cMCND.pdf");


  TCanvas *cMCdataND = new TCanvas("cMCdataND","cMCdataND", c1_x, c1_y);
  cMCdataND -> cd();
  hMCND_PECorr->SetTitle("");
  hMCND_PECorr->GetXaxis()->SetTitle("Corrected response/cm (PECorr/cm)");
  hMCND_PECorr->GetYaxis()->SetTitle("Hits in track window");

  hMCND_PECorr->SetLineColor(2);
  CenterTitles(hMCND_PECorr);
  hMCND_PECorr->Draw();
  hdataND_PECorr->Scale(hMCND_PECorr->Integral() / hdataND_PECorr->Integral());
  // Call below is commented out because I cannot find a definition for
  // SetTH1Max anywhere.  BJR
  //SetTH1Max(hMCND_PECorr, hdataND_PECorr);
  hdataND_PECorr->SetTitle("");
  CenterTitles(hdataND_PECorr);
  hdataND_PECorr->Draw("same");
  Preliminary();
  
  TLegend lMCdataND(0.6,0.55,0.85,0.7);  
  lMCdataND.AddEntry(hdataND_PECorr,"FD Data","l");
  lMCdataND.AddEntry(hMCND_PECorr, "FD MC","l");
  lMCdataND.SetBorderSize(0);
  lMCdataND.SetTextFont(42);
  lMCdataND.DrawClone();

  cMCdataND->Print("cMCdataND.pdf");

  
  TCanvas *cMCND_landau = new TCanvas("cMCND_landau","cMCND_landau", c1_x, c1_y);
  cMCND_landau -> cd();
  hMCND_PECorr_landau->SetTitle("");
  hMCND_PECorr_landau->GetXaxis()->SetTitle("Calibrated detector response/cm");
  hMCND_PECorr_landau->SetLineColor(2);
  CenterTitles(hMCND_PECorr_landau);
  hMCND_PECorr_landau->Draw();
  Simulation();

  
  TCanvas *cMCND_gaus = new TCanvas("cMCND_gaus","cMCND_gaus",c1_x,c1_y);
  cMCND_gaus -> cd();
  hMCND_PECorr_gaus->SetTitle("");
  hMCND_PECorr_gaus->GetXaxis()->SetTitle("Calibrated detector response/cm");
  hMCND_PECorr_gaus->SetLineColor(2);
  CenterTitles(hMCND_PECorr_gaus);
  hMCND_PECorr_gaus->Draw();
  Simulation();


  TCanvas *cMCND_truth = new TCanvas("cMCND_truth","cMCND_truth",c1_x,c1_y);
  cMCND_truth -> cd();
  gStyle->SetOptStat(0);
  hMCND_trueE->SetTitle("");
  hMCND_trueE->GetXaxis()->SetTitle("Simulated dE/dx (MeV/cm)");
  hMCND_trueE->GetYaxis()->SetTitle("Hits in track window");
  hMCND_trueE->SetLineColor(2);
  CenterTitles(hMCND_trueE);
  hMCND_trueE->Draw();
  Simulation();

  TLegend lTruthND(0.5,0.55,0.8,0.65);  
  lTruthND.AddEntry(hMCND_trueE,"FD MC","l");
  lTruthND.SetBorderSize(0);
  lTruthND.SetTextFont(42);
  lTruthND.DrawClone();

  cMCND_truth->Print("cMCND_truth.pdf");


  TCanvas *cMCdataND_cal = new TCanvas("cMCdataND_cal","cMCdataND_cal", c1_x, c1_y);
  cMCdataND_cal -> cd();
  hMCND_cal->SetTitle("");
  hMCND_cal->GetXaxis()->SetTitle("Calibrated dE/dx (MeV/cm)");
  hMCND_cal->GetYaxis()->SetTitle("Hits in track window");
  
  hMCND_cal->SetLineColor(2);
  CenterTitles(hMCND_cal);
  hMCND_cal->SetTitle("");
  hMCND_cal->Draw();
  hdataND_cal->Scale(hMCND_cal->Integral() / hdataND_cal->Integral());
  // Call below is commented out because I cannot find a definition for
  // SetTH1Max anywhere.  BJR
  //SetTH1Max(hMCND_cal, hdataND_cal);
  CenterTitles(hdataND_cal);
  hdataND_cal->Draw("same");
  Preliminary();
  
  TLegend lMCdataND_cal(0.5,0.55,0.8,0.7);  
  lMCdataND_cal.AddEntry(hdataND_cal,"FD Data","l");
  lMCdataND_cal.AddEntry(hMCND_cal, "FD MC","l");
  lMCdataND_cal.SetBorderSize(0);
  lMCdataND_cal.SetTextFont(42);
  lMCdataND_cal.DrawClone();

  cMCdataND_cal->Print("cMCdataND_cal.pdf");

  TCanvas *cMCND_truth_landau = new TCanvas("cMCND_truth_landau","cMCND_truth_landau", c1_x, c1_y);
  cMCND_truth_landau -> cd();
  hMCND_trueE_landau->SetTitle("");
  hMCND_trueE_landau->GetXaxis()->SetTitle("MeV/cm");
  hMCND_trueE_landau->SetLineColor(2);
  CenterTitles(hMCND_trueE_landau);
  hMCND_trueE_landau->Draw();
  Simulation();

  
  TCanvas *cMCND_truth_gaus = new TCanvas("cMCND_truth_gaus","cMCND_truth_gaus", c1_x, c1_y);
  cMCND_truth_gaus -> cd();

  hMCND_trueE_gaus->SetTitle("");
  hMCND_trueE_gaus->GetXaxis()->SetTitle("MeV/cm");
  hMCND_trueE_gaus->SetLineColor(2);
  CenterTitles(hMCND_trueE_gaus);
  hMCND_trueE_gaus->Draw();
  Simulation();


}




//landau
double Landau(TH1D* h, double& err){
  int bin1 = h -> FindFirstBinAbove(h -> GetMaximum()/2);
  int bin2 = h -> FindLastBinAbove(h -> GetMaximum()/2);
  double xLow = h -> GetBinCenter(bin1);
  double xHigh = h -> GetBinCenter(bin2);

  //fit with landau
  TF1 *f1 = new TF1("f1","landau",xLow,xHigh);
  h -> Fit(f1,"","",xLow,xHigh);
  
  double mpv = f1->GetParameter(1);
  err = f1->GetParError(1);
  std::cout<< "err: " << err <<std::endl;
  return mpv;
}


double Gaus(TH1D* h, double& err, bool isTruth){

  Double_t parameters[3];
  parameters[0] = h->GetBinContent(h->GetMaximumBin());
  parameters[1] = h->GetBinCenter(h->GetMaximumBin());
  parameters[2] = h->GetRMS();

  // only care about finding peak position, use gaussian truncated about mode  
  //Double_t bound_size = (isTruth) ? 0.2 : 0.5; // narrower for truth, *not* MC reco... don't mess this up.

  Double_t bound_size;
  if(isTruth) bound_size = 0.2;
  else bound_size = 0.5;

  std::cout<<" ****** bound size: " << bound_size <<std::endl;

  //Double_t lowerBound = max(0,parameters[1] - bound_size*h->GetRMS());
  Double_t lowerBound;
  if(parameters[1] - bound_size*h->GetRMS() > 0)
    lowerBound  = parameters[1] - bound_size*h->GetRMS();
  else lowerBound = 0;

  Double_t upperBound = parameters[1] + bound_size*h->GetRMS();

  TF1* fitfunction = new TF1("fitfunction","[0]*TMath::Gaus(x,[1],[2])",lowerBound,upperBound);
  fitfunction->SetParameters(parameters);

  // fitting                                                          
  h->Fit(fitfunction,"RQ");

  // retrieve values
  err = fitfunction->GetParError(1);
  return fitfunction->GetParameter(1);
}

double median(TH1D* h){
  return quantiles(h)[2];
}

std::vector<double> quantiles( TH1D* h ){
  //=================================                                             
  double q[5];
  double probs[5] = {0.025, 0.16, 0.5, 1 - 0.16, 0.975 };
  h->GetQuantiles( 5, q, probs );
  std::vector<double> r(5);
  for (int i=0; i<5; i++)
    {

      r[i] = q[i];
    }
  return r;
}

double errAoverB(double& err, double A, double B, double Aerr, double Berr){

  double C = A/B;
  //std::cout<< "C: " << C <<std::endl;

  err = std::sqrt( std::pow((Aerr/B),2) + std::pow((Berr*A)/(B*B),2) );
  //double err = std::sqrt(std::pow(C,2)) * std::sqrt( std::pow((Berr/B),2) + std::pow((Aerr/A),2));

  //std::cout<< "err: " << err <<std::endl;
  return err;

}


double GetHistMax(TH1* h0, TH1* h1, TH1* h2, TH1* h3){
  double max[4];
  max[0] = h0->GetMaximum();
  max[1] = h1->GetMaximum();
  max[2] = h2->GetMaximum();
  max[3] = h3->GetMaximum();

  double hiy = 0;
  for(int i=1;i<4;i++){
    if(max[i]>hiy)
      hiy = max[i];
  }
  return (hiy * 1.1);
}

int absCal(){
  return main();
}