DedxGamma.C

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cstring>

#include "TFile.h"
#include "TChain.h"
#include "TH2.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TProfile.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGaxis.h"


TH1D* MakeRatio(TH1D* num, TH1D* denom,const char* name)
{

  TH1D* ratio = new TH1D(name,";;Data/MC",num->GetNbinsX(),num->GetXaxis()->GetXbins()->GetArray());

  ratio->GetXaxis()->SetTitle(num->GetXaxis()->GetTitle());

  ratio->Divide(num,denom);

  ratio->SetLineColor(denom->GetLineColor());
  ratio->SetMarkerColor(denom->GetMarkerColor());
  ratio->SetMarkerSize(denom->GetMarkerSize());

  return ratio;

}


TCanvas* MakeSplitRatio(TCanvas* top, TCanvas* bot,const char* name)
{

  TCanvas* can = new TCanvas(name,name,600,550);

  double pixels = 18.0;
  double labPix = 15.0;
  double labOff = 5.0;

  can->SetBottomMargin(0);
  can->Divide(1,2,0,0);
  TVirtualPad* can_1 = can->cd(1);
  can_1->SetPad(0,0.4,1,1);
  can_1->SetFillStyle(0);
  top->DrawClonePad();
  can_1->SetRightMargin(0.05);
  can_1->SetTopMargin(0.1);
  can_1->SetBottomMargin(0.005);
  can_1->SetLeftMargin(0.13);
  can_1->Update();
  double textSizeFac = 1.0/(fabs(can_1->YtoAbsPixel(can_1->GetY2()) - can_1->YtoAbsPixel(can_1->GetY1())));
  double marginSize = 80;
  double textsizeA = pixels/(can_1->GetAbsHNDC()*can_1->GetWh());
  double textsizeAl = labPix/(can_1->GetAbsHNDC()*can_1->GetWh());
  double textsizeAo = labOff/(can_1->GetWw() * can_1->GetAbsWNDC());

  TList* primlist = can_1->GetListOfPrimitives();
  for(int ilist = 0; ilist < primlist->GetSize(); ++ilist){
    TObject* prim = primlist->At(ilist);
    if(prim->InheritsFrom("TH1")){
      TH1* primhist = (TH1*)prim;
      primhist->GetYaxis()->SetTitleSize(textsizeA);
      primhist->GetYaxis()->SetLabelSize(textsizeAl);
      primhist->GetYaxis()->SetLabelOffset(textsizeAo);
      primhist->GetYaxis()->SetTitleOffset(0.035 / textsizeA);   
      primhist->GetXaxis()->SetNoExponent(kTRUE);
    }
  }

  can_1->Update();

  TVirtualPad* can_2 = can->cd(2);
  can_2->SetPad(0,0,1,0.4);
  bot->DrawClonePad();
  can_2->SetTopMargin(0.025);
  can_2->SetBottomMargin(0.22);
  can_2->SetRightMargin(0.05);
  can_2->SetLeftMargin(0.13);
  can_2->Update();

  double textsize = pixels/(can_2->GetAbsHNDC()*can_2->GetWh());
  double textsizel = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
  double textsizeo = labOff/(can_2->GetWw() * can_2->GetAbsWNDC());
  double textsizex = pixels/(can_2->GetAbsHNDC()*can_2->GetWh());
  double textsizelx = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
  textsizeAl = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
  double textsizeox = labOff/(can_2->GetWw() * can_2->GetAbsWNDC());

  TList* primlist_2 = can_2->GetListOfPrimitives();
  for(int ilist = 0; ilist < primlist_2->GetSize(); ++ilist){
    TObject* prim = primlist_2->At(ilist);
    if(prim->InheritsFrom("TH1")){
      TH1* primhist = (TH1*)prim;
      primhist->GetYaxis()->SetTitleSize(textsize);
      primhist->GetYaxis()->SetLabelSize(textsizel);
      primhist->GetYaxis()->SetLabelOffset(textsizeAo);
      primhist->GetXaxis()->SetTitleSize(textsizex);
      primhist->GetXaxis()->SetLabelSize(textsizelx);
      primhist->GetXaxis()->SetLabelOffset(textsizeox);
      primhist->GetYaxis()->SetTitleOffset(0.035 / textsize);  
      primhist->GetXaxis()->SetTitleOffset(0.9);
      primhist->GetXaxis()->SetNoExponent(kTRUE);
    }
  }

  can_2->Update();

  return can;

}

void FindTruncate(TH1D* hist, double& low, double& high)
{
  //std::cout<<"BEGINING OF TRUNCATE FUNCATION!"<<std::endl;
  double percentCut = 0.05;
  double histMin = hist->GetXaxis()->GetXmin();
  double histMax = hist->GetXaxis()->GetXmax();
  int minbin = hist->GetXaxis()->FindFixBin(histMin);
  int maxbin = hist->GetXaxis()->FindFixBin(histMax)-1;
  double integral = hist->Integral(minbin,maxbin);
  //std::cout<<"Total integral is: "<<integral<<std::endl;
  //if (integral == 363){
  //TCanvas* cs = new TCanvas("cs","cs");
  //hist->Draw();
  //cs->Print("test.pdf");
  //}
  double lowCut = integral*percentCut;
  double highCut = integral*(1-percentCut);

  int lowbin = 0;
  int highbin = 0;
  bool foundLow = false;
  double ongoingIntegral = 0.0;

  for(int ibin = minbin; ibin <= maxbin; ++ibin){
    double bincont = hist->GetBinContent(ibin);
    //std::cout<<"bin: "<<ibin<<" with content: "<<bincont<<std::endl;
    if (!foundLow){
      ongoingIntegral += bincont;
      //std::cout<<" ongoingIntegral: "<<ongoingIntegral<<std::endl;
      if (ongoingIntegral > lowCut){
        lowbin = ibin;
        foundLow = true;
        //std::cout<<"Found Low! Set lowbin : "<<lowbin<<std::endl;
      }
    }
    else{
      if ((ongoingIntegral<highCut)&&((ongoingIntegral+bincont)>=highCut)){
        ongoingIntegral += bincont;
        highbin = ibin;
        //std::cout<<"Found High! Set highbin : "<<highbin<<std::endl;
      }
      else if (ongoingIntegral<highCut){
        ongoingIntegral += bincont;
        //std::cout<<" ongoingIntegral: "<<ongoingIntegral<<std::endl;
      } 
    }//Already found low, now looking for high bin
    //Now need to set low, high to exclude the lowbin and highbin
    low = hist->GetXaxis()->GetBinUpEdge(lowbin);
    high = hist->GetXaxis()->GetBinLowEdge(highbin);
  }
  return;
}

double Mode(TH1D* hist)
{

  double modex = hist->GetXaxis()->GetBinCenter(hist->GetMaximumBin());
  //std::cout<<"In MODE with Modex: "<<modex<<std::endl;
  //SL -- adding this line, because I think it needs something ....
  return modex;
}

double ModeRestrict(TH1D* hist, double minx, double maxx)
{

  double mode = Mode(hist);
  //std::cout<<"In mode restrict with Mode: "<<mode<<" maxx: "<<maxx<<" minx: "<<minx<<std::endl;
  if( mode > maxx || mode < minx){
    // loop over bins between maxx and minx to find the mode in the restricted range;
    //std::cout<<"HERE HERE HERE!!!"<<std::endl;
    int minxbin = hist->GetXaxis()->FindFixBin(minx);
    int maxxbin = hist->GetXaxis()->FindFixBin(maxx);
    double modeval = hist->GetBinContent(minxbin);
    int modebin = minxbin;
//     std::cout<<"minxbin: "<<minxbin<<" maxxbin: "<<maxxbin<<" modeval: "<<modeval<<std::endl;
    for(int ibin = minxbin+1; ibin <= maxxbin; ++ibin){
      double bincont = hist->GetBinContent(ibin);
//       std::cout<<"ibin: "<<ibin<<" bincont: "<<bincont<<std::endl;
      if(bincont > modeval){ 
        modebin = ibin;
        modeval = bincont;
        //std::cout<<"modeval: "<<modeval<<" modebin: "<<modebin<<std::endl;
      }
    }
    mode = hist->GetXaxis()->GetBinCenter(modebin);
  }

  return mode;

}


void  GausFit(TH1D* hist, double& mean, double& err, double minrange, double maxrange, double lowsig, double highsig, bool useTruncate)
{
   std::cout<<"        We are inside GausFit! "<<std::endl;
  std::cout<<"        we have minrange = "<<minrange<<" and maxrange = "<<maxrange<<std::endl;
  std::cout<<"        we have lowsig = "<<lowsig<<" and highsig = "<<highsig<<std::endl;
  double mode = ModeRestrict(hist,minrange,maxrange);
  std::cout<<"        we have mode = "<<mode<<std::endl;
  double fitminx = -1.0;
  double fitmaxx = -1.0;
  if (useTruncate){
    fitminx = lowsig;
    fitmaxx = highsig;
  }
  else{
    fitminx = mode-lowsig;
    fitmaxx = mode+highsig;
  }
  //std::cout<<"        we have fitminx = "<<fitminx<<" and fitmaxx = "<<fitmaxx<<std::endl;
  if(fitminx < minrange){ fitminx = minrange; }
  if(fitmaxx > maxrange){ fitmaxx = maxrange; }
  //std::cout<<"        we have fitminx = "<<fitminx<<" and fitmaxx = "<<fitmaxx<<std::endl;  
  TF1* fgaus = new TF1("gaus","gaus(0)",fitminx,fitmaxx);
  std::cout<<"fitminx: "<<fitminx<<" fitmaxx: "<<fitmaxx<<" mode: "<<mode<<std::endl;
  fgaus->SetParameters(hist->GetMaximum(),mode,(fitmaxx-fitminx)/2.0);
  int minbin = hist->GetXaxis()->FindFixBin(fitminx);
  int maxbin = hist->GetXaxis()->FindFixBin(fitmaxx);
  double integral = hist->Integral(minbin,maxbin);
  mean = -1;
  err = 0.000001;
  if(integral > 100){
    std::cout<<"        integral is "<<integral<<std::endl;
    if (integral == 320){
      TCanvas* cs = new TCanvas("cs","cs");
      hist->Fit(fgaus,"RIQ");
      cs->Print("test.pdf");
    }
    int fitstatus =  hist->Fit(fgaus,"RI");
    std::cout<<"And our fit status is: "<<fitstatus<<std::endl;
    if(fitstatus == 0){
      mean = fgaus->GetParameter(1);
      err = fgaus->GetParError(1);
      if(mean > maxrange || mean < minrange){
        mean = -1;
        err = 0.000001;
      }
      //err = fgaus->GetParameter(2);
    }
  }

   std::cout<<"mean: "<<mean<<" err: "<<err<<std::endl;

}

TH1D* PeakProjectionY(TH2F* hist,const char* name, double minrange = -1.0, double maxrange = -1.0, double lowsig = 0.001, double highsig = 0.001, bool useTruncate = true)
{
  std::cout<<"In peak projection!"<<std::endl;
  TH1D* peakhist = new TH1D(name,"",hist->GetNbinsX(),hist->GetXaxis()->GetXbins()->GetArray());
  peakhist->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
  peakhist->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
  //std::cout<<"And the axis title is: "<<hist->GetXaxis()->GetTitle()<<std::endl;
  double fitmin = minrange;
  double fitmax = maxrange;
  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
  if(minrange == -1.0){ fitmin = hist->GetXaxis()->GetXmin(); }
  if(maxrange == -1.0){ fitmax = hist->GetXaxis()->GetXmax(); }
  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
  //std::cout<<"Number of x bins: "<<hist->GetNbinsX()<<std::endl;
  for(int ixbin = 1; ixbin <= hist->GetNbinsX(); ++ixbin){
    // get the projection of the 2d hist
    TH1D* projection = hist->ProjectionY("",ixbin,ixbin,"o");
    //std::cout<<"  X bin: "<<ixbin<<" Entries: "<<projection->GetEntries()<<std::endl;
    if(projection->GetEntries() == 0){ continue; }
    // fit the projection to a Gaussian around the peak
    double mean = 0.00175 ,err = 0.001;
    double lowTruncate = -1.0;
    double highTruncate = -1.0;
    if (useTruncate){
      FindTruncate(projection, lowTruncate, highTruncate);
      if ((lowTruncate != -1.0) && (highTruncate != -1.0)){
        lowsig = lowTruncate;
        highsig = highTruncate;
      }
    }
    //std::cout<<"Low truncate: "<<lowTruncate<<" High truncate: "<<highTruncate<<std::endl;
    GausFit(projection,mean,err,fitmin,fitmax,lowsig,highsig,useTruncate);
    //std::cout<<"  Mean: "<<mean<<" Error: "<<err<<std::endl;
    peakhist->SetBinContent(ixbin,mean);
    peakhist->SetBinError(ixbin,err);

  }
  
  return peakhist;

}


void DedxGamma()
{

  //const char* filenamedata = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/datadedxVxFineS15-05-22.root";
  //const char* filenamemc = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/mcdedxVxFineS15-05-22.root";

  const char* filenamedata = "/local/nova/users/lein/datadedxVgammaFineS15-05-22.root";
  const char* filenamemc = "/local/nova/users/lein/mcdedxVgammaFineS15-05-22.root";

  const char* filenamedataP = "/local/nova/users/lein/datadedxVgammaFineS15-05-22Proton.root";
  const char* filenamemcP = "/local/nova/users/lein/mcdedxVgammaFineS15-05-22Proton.root";

  TFile* filedata = new TFile(filenamedata,"READ");
  TFile* filemc = new TFile(filenamemc,"READ");

  TFile* filedataP = new TFile(filenamedataP,"READ");
  TFile* filemcP = new TFile(filenamemcP,"READ");

  TH2F* dedxdata = (TH2F*)filedata->Get("DedxVx");
  TH2F* dedxmc = (TH2F*)filemc->Get("DedxVx");

  TH2F* dedxdataP = (TH2F*)filedataP->Get("DedxVx");
  TH2F* dedxmcP = (TH2F*)filemcP->Get("DedxVx");

  //std::cout<<"This is the pointer: "<<dedxdata<<std::endl;

  gStyle->SetOptStat(00000);

  TH1D* dedxpeakdata =  PeakProjectionY(dedxdata,"dedxpeakdata");
  TH1D* dedxpeakmc =  PeakProjectionY(dedxmc,"dedxpeakmc");

  TH1D* dedxpeakdataP =  PeakProjectionY(dedxdataP,"dedxpeakdataP");
  std::cout<<"About to do the one!"<<std::endl;
  TH1D* dedxpeakmcP =  PeakProjectionY(dedxmcP,"dedxpeakmcP");

  //If you want to run WITHOUT the truncation and just limit the fit range to +- 0.001 around the mode, do these two lines 
  //instead of the two lines above
  //TH1D* dedxpeakdata =  PeakProjectionY(dedxdata,"dedxpeakdata",-1.0,-1.0,0.001,0.001,false);
  //TH1D* dedxpeakmc =  PeakProjectionY(dedxmc,"dedxpeakmc",-1.0,-1.0,0.001,0.001,false);


  TCanvas* c0 = new TCanvas("c0","c0");
  dedxpeakdata->Draw();
  dedxpeakmc->Draw("same");
  dedxpeakmcP->Draw("same");
  dedxpeakdataP->Draw("same");
  dedxpeakmc->SetLineColor(kRed);
  dedxpeakmc->SetMarkerColor(kRed);
  dedxpeakmcP->SetLineColor(kViolet);
  dedxpeakmcP->SetMarkerColor(kViolet);
  dedxpeakdataP->SetLineColor(kGreen+2);
  dedxpeakdataP->SetMarkerColor(kGreen+2);
  dedxpeakdata->GetYaxis()->SetRangeUser(0.0,0.006);
  dedxpeakdata->GetXaxis()->SetNoExponent();
  //dedxpeakdata->GetYaxis()->SetNoExponent();
  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
  leg0->AddEntry(dedxpeakdata,"Muon Data","LE");
  leg0->AddEntry(dedxpeakmc,"Muon MC","LE");
  leg0->AddEntry(dedxpeakdataP,"Proton Data","LE");
  leg0->AddEntry(dedxpeakmcP,"Proton MC","LE");
  leg0->Draw("same");
  c0->Print("dedxfine.pdf");
  c0->Print("dedxfine.root");

  TCanvas* c00 = new TCanvas("c00","c00");
  dedxdata->Draw("colz");
  dedxpeakdata->Draw("same");
  dedxdata->GetYaxis()->SetRangeUser(0.0,0.006);
  dedxdata->GetXaxis()->SetNoExponent();
  //dedxdata->GetYaxis()->SetNoExponent();
  c00->Print("dedxDataMother.pdf");

  TCanvas* c000 = new TCanvas("c000","c000");
  //dedxpeakmc->SetLineColor(kBlack);
  dedxmc->Draw("colz");
  dedxpeakmc->Draw("same");
  dedxmc->GetYaxis()->SetRangeUser(0.0,0.006);
  dedxmc->GetXaxis()->SetNoExponent();
  //dedxmc->GetYaxis()->SetNoExponent();
  c000->Print("dedxMCMother.pdf");

  TCanvas* c00P = new TCanvas("c00P","c00P");
  dedxdataP->Draw("colz");
  dedxpeakdataP->Draw("same");
  dedxdataP->GetYaxis()->SetRangeUser(0.0,0.006);
  dedxdataP->GetXaxis()->SetNoExponent();
  //dedxdata->GetYaxis()->SetNoExponent();
  c00P->Print("dedxDataProtonMother.pdf");

  TCanvas* c000P = new TCanvas("c000P","c000P");
  //dedxpeakmc->SetLineColor(kBlack);
  dedxmcP->Draw("colz");
  dedxpeakmcP->Draw("same");
  dedxmcP->GetYaxis()->SetRangeUser(0.0,0.006);
  dedxmcP->GetXaxis()->SetNoExponent();
  //dedxmc->GetYaxis()->SetNoExponent();
  c000P->Print("dedxMCProtonMother.pdf");

  
  int nsteps = 1000;
  double minscale = 0.95;
  double maxscale = 1.05;

  TH1F* chi2match = new TH1F("chi2match",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
  TH1F* chi2matchrestrict = new TH1F("chi2matchrestrict",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);


  for(int ibin = 1; ibin <= chi2match->GetNbinsX(); ++ibin){
    // Get the scale factor
    double scalefactor = chi2match->GetXaxis()->GetBinCenter(ibin);
    //calculate a chi2 between data and scaled mc
    double chi2 = 0;
    double chi2rest = 0;
    for(int jbin = 1; jbin <= dedxpeakdata->GetNbinsX(); ++jbin){
      // get the dist from the end of the track
      double dist = dedxpeakdata->GetXaxis()->GetBinCenter(jbin);
      double expt = dedxpeakdata->GetBinContent(jbin);
      double errdata = dedxpeakdata->GetBinError(jbin);
      double errmc = dedxpeakmc->GetBinError(jbin);
      double err2 = errdata*errdata+errmc*errmc;
      double obsunscale = dedxpeakmc->GetBinContent(jbin);
      double obs = scalefactor*dedxpeakmc->GetBinContent(jbin);
      double diff = expt-obs;
      //chi2+=(diff*diff/expt);
      // remove error cases
      if(expt < 0 || obsunscale < 0){ continue; }
      chi2+=(diff*diff/err2);
      // std::cout<<"dist: "<<dist<<std::endl;
      if(dist>100 && dist < 200){
        chi2rest+=(diff*diff/err2);
      }

    }
//     // make a histogram to store the scaled version of mc dedxpeak histogram
//     TString name;
//     name+=ibin;
//     TH1D* dedxpeakmcscale = new TH1D(name.Data(),"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
//     dedxpeakmcscale->Add(dedxpeakmc);
//     dedxpeakmcscale->Scale(scalefactor);
// dedxpeakmcscale->
//     double chi2 = dedxpeakdata->Chi2Test(dedxpeakmcscale,"CHI2");
//     std::cout<<"scalefactor: "<<scalefactor<<" chi2: "<<chi2<<std::endl;
    chi2match->SetBinContent(ibin,chi2);
    chi2matchrestrict->SetBinContent(ibin,chi2rest);

  }

  TCanvas* c00a = new TCanvas("c00a","c00a");
  chi2match->Draw();
  c00a->Print("chi2fine.pdf");

  double optscale = chi2match->GetXaxis()->GetBinCenter(chi2match->GetMinimumBin());
  std::cout<<"optscale: "<<optscale<<std::endl;

  TCanvas* c000a = new TCanvas("c000a","c000a");
  chi2matchrestrict->Draw();
  c000a->Print("chi2restrictfine.pdf");

  double optscalerestrict = chi2matchrestrict->GetXaxis()->GetBinCenter(chi2matchrestrict->GetMinimumBin());
  std::cout<<"Restricted optscale: "<<optscalerestrict<<std::endl;

  TH1D* dedxpeakmcscale = new TH1D("OptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
  dedxpeakmcscale->Add(dedxpeakmc);
  dedxpeakmcscale->Scale(optscale);

  TCanvas* c2 = new TCanvas("c2","c2");
  dedxpeakdata->Draw();
  dedxpeakmcscale->Draw("same");
  dedxpeakmcscale->SetLineColor(kRed);
  dedxpeakmcscale->SetMarkerColor(kRed);
  dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
  dedxpeakdata->GetXaxis()->SetNoExponent();
  //dedxpeakdata->GetYaxis()->SetNoExponent();
  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
  leg0->AddEntry(dedxpeakdata,"Data","LE");
  leg0->AddEntry(dedxpeakmcscale,"MC","LE");
  leg0->Draw("same");
  c2->Print("dedxoptfine.pdf");

  double nonoptscale = 0.97375; //1.483/1.55;

  TH1D* dedxpeakmcnonoptscale = new TH1D("NonOptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
  dedxpeakmcnonoptscale->Add(dedxpeakmc);
  dedxpeakmcnonoptscale->Scale(nonoptscale);

  TCanvas* c3 = new TCanvas("c3","c3");
  dedxpeakdata->Draw();
  dedxpeakmcnonoptscale->Draw("same");
  dedxpeakmcnonoptscale->SetLineColor(kRed);
  dedxpeakmcnonoptscale->SetMarkerColor(kRed);
  dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
  dedxpeakdata->GetXaxis()->SetNoExponent();
  //dedxpeakdata->GetYaxis()->SetNoExponent();
  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
  leg0->AddEntry(dedxpeakdata,"Data","LE");
  leg0->AddEntry(dedxpeakmcnonoptscale,"MC","LE");
  leg0->Draw("same");
  c3->Print("dedxnonoptfine.pdf");


  TCanvas* c4 = new TCanvas("c4","c4");
  TH1D* dedxratio = MakeRatio(dedxpeakdata,dedxpeakmc,"dedxratio");
  TLine* unitline = new TLine(-1,1,1001,1);
  unitline->SetLineWidth(2);
  unitline->SetLineStyle(7);
  dedxratio->Draw();
  unitline->Draw("same");
  dedxratio->GetYaxis()->SetRangeUser(0.9,1.1);

  TCanvas* c5 = MakeSplitRatio(c0,c4,"c5");
  c5->Print("dedxratio.pdf");

  


}