#include <iostream>
#include <iomanip>
using namespace std;
#include "TH1D.h"
#include "TF1.h"
#include <TSystem.h>
#include <TCanvas.h>
#include <TApplication.h>
#include <TF1.h>
#include <TH1.h>
#include <TAxis.h>
#include <TLegend.h>
#include <TLine.h>
#include <TStyle.h>
#include <TPaveText.h>
#include "Math/IntegratorOptions.h"
#include "Math/MinimizerOptions.h"
#include "TMinuit.h"
#include "TRandom.h"
#include "fitTools.h"
#include "my_fitTools.h"
#include "fw_fitTools.h"

//TH1D *hd ;
//int pseudo(){
int main(){

  TH1D *LLR0 = new TH1D("LLR0","LLR0",300,20,80);
  TH1D *LLsigma = new TH1D("LLsigma","LLsigma",300,0,10);
  TH1D *LLkb = new TH1D("LLkb","LLkb",300,0,0.34);

  double resoR0=2;
  double resosigma=4;

  int dnbin=40;
  
  TH1D *diforiR0= new TH1D("diforiR0","diforiR0",dnbin, -resoR0, resoR0);
  TH1D *diforisigma= new TH1D("diforisigma","diforisigma",dnbin, -resosigma, resosigma);
  TH1D *difminuR0= new TH1D("difminuR0","difminuR0",dnbin, -resoR0, resoR0);
  TH1D *difminusigma= new TH1D("difminusigma","difminusigma",dnbin, -resosigma, resosigma);
  TH1D *diffwR0= new TH1D("diffwR0","diffwR0",dnbin, -resoR0, resoR0);
  TH1D *diffwsigma= new TH1D("diffwsigma","diffwsigma",dnbin, -resosigma, resosigma);


  TH1D *minuoriR0= new TH1D("minuoriR0","minuoriR0",dnbin, -resoR0, resoR0);
  TH1D *minuorisigma= new TH1D("minuorisigma","minuorisigma",dnbin, -resosigma, resosigma);
  TH1D *fworiR0= new TH1D("fworiR0","fworiR0",dnbin, -resoR0, resoR0);
  TH1D *fworisigma= new TH1D("fworisigma","fworisigma",dnbin, -resosigma, resosigma);


  TH1D *minupullR0= new TH1D("minupullR0","minupullR0",dnbin+1,-11,11);
  TH1D *minupullsigma= new TH1D("minupullsigma","minupullsigma",dnbin+1,-11,11);



  
  hd = new TH1D("hd","hd", nbin, bins);
  // hd->SetBinContent(0  , 0          ); 
  // hd->SetBinContent(1  , 4.53902    ); 
  // hd->SetBinContent(2  , 4.70824    ); 
  // hd->SetBinContent(3  , 4.73104    ); 
  // hd->SetBinContent(4  , 4.8716     ); 
  // hd->SetBinContent(5  , 4.94311    ); 
  // hd->SetBinContent(6  , 4.92585    ); 
  // hd->SetBinContent(7  , 5.0351     ); 
  // hd->SetBinContent(8  , 5.25363    ); 
  // hd->SetBinContent(9  , 5.30603    ); 
  // hd->SetBinContent(10 ,  5.44443   ); 
  // hd->SetBinContent(11 ,  5.53738   ); 
  // hd->SetBinContent(12 ,  5.54719   ); 
  // hd->SetBinContent(13 ,  5.76489   ); 
  // hd->SetBinContent(14 ,  5.94505   ); 
  // hd->SetBinContent(15 ,  6.15808   ); 
  // hd->SetBinContent(16 ,  6.43717   ); 
  // hd->SetBinContent(17 ,  6.73969   ); 
  // hd->SetBinContent(18 ,  6.9798    ); 
  // hd->SetBinContent(19 ,  7.25404   ); 
  // hd->SetBinContent(20 ,  7.69753   ); 
  // hd->SetBinContent(21 ,  8.0544    ); 
  // hd->SetBinContent(22 ,  8.73745   ); 
  // hd->SetBinContent(23 ,  9.82025   ); 
  // hd->SetBinContent(24 ,  11.7717   ); 
  // hd->SetBinContent(25 ,  13.485    ); 
  // hd->SetBinContent(26 ,  7.77343   ); 
  // hd->SetBinContent(27 ,  1.36615   ); 
  // hd->SetBinContent(28 ,  1.25825   ); 
  // hd->SetBinContent(29 ,  0.556537  );  
  // hd->SetBinContent(30 ,  0.316932  );  
  // hd->SetBinContent(31 ,  0.348092  );  
  // hd->SetBinContent(32 ,  0.472119  );  
  // hd->SetBinContent(33 ,  0  	    ); 

  // hd->SetBinError(0,   0            );
  // hd->SetBinError(1,   0.41905      );
  // hd->SetBinError(2,   0.480005     );
  // hd->SetBinError(3,   0.256493     );
  // hd->SetBinError(4,   0.334958     );
  // hd->SetBinError(5,   0.368144     );
  // hd->SetBinError(6,   0.214421     );
  // hd->SetBinError(7,   0.272643     );
  // hd->SetBinError(8,   0.353102     );
  // hd->SetBinError(9,   0.34398      );
  // hd->SetBinError(10,   0.345476    );
  // hd->SetBinError(11,   0.289041    );
  // hd->SetBinError(12,   0.26208     );
  // hd->SetBinError(13,   0.212691    );
  // hd->SetBinError(14,   0.203182    );
  // hd->SetBinError(15,   0.310943    );
  // hd->SetBinError(16,   0.360032    );
  // hd->SetBinError(17,   0.489854    );
  // hd->SetBinError(18,   0.461181    );
  // hd->SetBinError(19,   0.639442    );
  // hd->SetBinError(20,   0.510431    );
  // hd->SetBinError(21,   0.426875    );
  // hd->SetBinError(22,   0.232739    );
  // hd->SetBinError(23,   0.310051    );
  // hd->SetBinError(24,   0.274143    );
  // hd->SetBinError(25,   0.569625    );
  // hd->SetBinError(26,   0.410178    );
  // hd->SetBinError(27,   0.215306    );
  // hd->SetBinError(28,   0.229667    );
  // hd->SetBinError(29,   0.107403    );
  // hd->SetBinError(30,   0.0644741   );
  // hd->SetBinError(31,   0.0677504   );
  // hd->SetBinError(32,   0.0930139   );
  // hd->SetBinError(33,   0  	    );

  QuenchedBragg qb; //TF1 Functors
  TF1* fphotons = new TF1("Fitted Light Output",&qb,0, 78.6606, 5);
  fphotons->SetParNames("R_{0}","#sigma_{R}","#Phi_{0}","xOffset","kB");
  fphotons->FixParameter(3,0.0);
  fphotons->SetNpx(100);
  fphotons->SetRange(0, 92.4508);
  
  //R0         20-80
  //sigma       1-9
  //pkb    0.03-0.1

  // input
  double pR0=70;
  double psigma=2.0;
  double pphi0=6e-7;
  double pkb=0.07;

  //initialization
  double initR0;
  double initsigma;
  double initphi0;
  double initkb;

  // output
  double oriR0=0;
  double orisigma=0;
  double oriphi0=0;
  double orikb=0;

  double minuR0=0;
  double minusigma=0;
  double minuphi0=0;
  double minukb=0;

  double fwR0=0;
  double fwsigma=0;
  double fwphi0=0;
  double fwkb=0;
  
  gRandom->SetSeed(13);

  //  int ntest=256;

  int ntest=32;

  
  for(int nt=0; nt<ntest; nt++){
    cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
    cout<<"HHHHHHHHHHHHHHHHHH          "<<nt<<" th         HHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
    cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
    //    cout<<"TRandom = "<<gRandom->Rndm()<<endl;
    // double pR0=70;
    // double psigma=2.0;
    // double pphi0=6e-7;
    // double pkb=0.07;
    //R0         20-80
    //sigma       1-9
    //pkb    0.03-0.08
    
    pR0=20 + 60*gRandom->Rndm(); //20-80
    psigma= 1+2*gRandom->Rndm();  //1-3
    pphi0= 6e-7;
    pkb = 0.04+0.06*gRandom->Rndm();  //0.04-0.10

    cout<<"HHH input parameters HHH"<<endl;
    cout<<"pR0="<<pR0<<endl;
    cout<<"psigma="<<psigma<<endl;
    cout<<"pphi0="<<pphi0<<endl;
    cout<<"pkb="<<pkb<<endl;

    fphotons->SetParameters(pR0,psigma,pphi0,0,pkb);
    //    cout<<"final integral="<<fphotons->Integral(0, 92.4508)<<endl;
    //  int npse=1000000;

    hd->Reset();

    int npse=10000000;  
    
    for(int np=0; np<npse; np++){
      hd->Fill(fphotons->GetRandom());
    }
    //temporally error from real data
    //  hd->Scale(fphotons->Integral(0, 92.4508 /hd->Integral())/2.5);
    hd->Scale(fphotons->Integral(0, 92.4508) / hd->Integral());  
    //  cout<<hd->Integral()<<"   "<<fphotons->Integral(0, 92.4508)<<endl;
    hd->SetBinError(0,   0            );
    hd->SetBinError(1,   0.41905      );
    hd->SetBinError(2,   0.480005     );
    hd->SetBinError(3,   0.256493     );
    hd->SetBinError(4,   0.334958     );
    hd->SetBinError(5,   0.368144     );
    hd->SetBinError(6,   0.214421     );
    hd->SetBinError(7,   0.272643     );
    hd->SetBinError(8,   0.353102     );
    hd->SetBinError(9,   0.34398      );
    hd->SetBinError(10,   0.345476    );
    hd->SetBinError(11,   0.289041    );
    hd->SetBinError(12,   0.26208     );
    hd->SetBinError(13,   0.212691    );
    hd->SetBinError(14,   0.203182    );
    hd->SetBinError(15,   0.310943    );
    hd->SetBinError(16,   0.360032    );
    hd->SetBinError(17,   0.489854    );
    hd->SetBinError(18,   0.461181    );
    hd->SetBinError(19,   0.639442    );
    hd->SetBinError(20,   0.510431    );
    hd->SetBinError(21,   0.426875    );
    hd->SetBinError(22,   0.232739    );
    hd->SetBinError(23,   0.310051    );
    hd->SetBinError(24,   0.274143    );
    hd->SetBinError(25,   0.569625    );
    hd->SetBinError(26,   0.410178    );
    hd->SetBinError(27,   0.215306    );
    hd->SetBinError(28,   0.229667    );
    hd->SetBinError(29,   0.107403    );
    hd->SetBinError(30,   0.0644741   );
    hd->SetBinError(31,   0.0677504   );
    hd->SetBinError(32,   0.0930139   );
    hd->SetBinError(33,   0  	    );



    //  double R0 =hd->GetBinCenter(hd->GetMaximumBin())+10; //mod aa.... I see from back is important.  not  good for these dependence
    //  double R0 = EstimateR0(hd);
    double R0 =hd->GetBinCenter(hd->GetMaximumBin())+15; //mod aa.... I see from back is important.  not  good for these dependence
    double E0 = pow(R0/alpha,1/p)/mp;
    //  double E0 = 0;
    double xOffset = 0.0;
    double sigma_mono = 0.012*pow(R0,0.935);
    double sigma_E0 = 0.002*E0;                                         //in MeV. realistic initial beam energy spread?
    double sigma_0 = pow(sigma_E0*alpha*p,2)*pow(E0,2*p-2);             //in mm, estimated offset due to detector effects and initial beam energy spread
    double sigma = sqrt(pow(sigma_mono,2) + pow(sigma_0,2));
    double phi0 = 1e-7*hd->GetBinContent(hd->GetMaximumBin());  //rough guess
    double fit_start = 0;
    if(R0 < 30) fit_start = (R0-xOffset-7)/1.8;                         // protons hit sensor directly at low beam energies/ranges
    double fit_end = hd->GetXaxis()->GetBinLowEdge(hd->GetMaximumBin()+3); //light back to zero normally 3 sheets after peak
    double kb = 0.07;
    fphotons->SetParameters(R0, sigma, phi0, 0, kb); // i don't know only this work...
    //    cout<<"init para = "<<R0<<" "<<sigma<<" "<<phi0<<" "<<xOffset<<" "<<kb<<endl;


    double initR0=hd->GetBinCenter(hd->GetMaximumBin())+10;
    double initsigma=2;
    double initphi0=20e-7;
    double initkb=0.07;
    
    // Special check 
    //    fphotons->SetParameters(initR0, initsigma, initphi0, 0, initkb); // no it error


    E0 = pow(fphotons->GetParameter(0)/alpha,1/p)/mp;
    cout<<"Initial value -----------------------------------"<<endl;
    cout<<"PreCalc R0    = "<<fphotons->GetParameter(0)<<endl;
    cout<<"PreCalc sigma = "<<fphotons->GetParameter(1)<<endl;
    cout<<"PreCalc phi   = "<<fphotons->GetParameter(2)<<endl;
    //  cout<<"PreCalc off   = "<<fphotons->GetParameter(3)<<endl;  
    cout<<"PreCalc kB    = "<<fphotons->GetParameter(4)<<endl;  
    cout<<"PreCalc E0    = "<<E0<<endl;

    hd->Fit("Fitted Light Output","IRQ0"); //updates roughly 5x faster without I option		// out to reduce the test time
    E0 = pow(fphotons->GetParameter(0)/alpha,1/p)/mp;
    cout<<"original fit results -----------------------------------"<<endl;
    cout<<"Orig Fit R0    = "<<fphotons->GetParameter(0)<<endl;
    cout<<"Orig Fit sigma = "<<fphotons->GetParameter(1)<<endl;
    cout<<"Orig Fit phi   = "<<fphotons->GetParameter(2)<<endl;
    cout<<"Orig Fit kB    = "<<fphotons->GetParameter(4)<<endl;  
    cout<<"Orig  E0       = "<<E0<<endl;
    
    oriR0=fphotons->GetParameter(0);
    orisigma=fphotons->GetParameter(1);
    oriphi0=fphotons->GetParameter(2);
    orikb=fphotons->GetParameter(4);

    
  
    //HH Full Minuit HH//
    const int npar=4;
    int ier=0;
    TMinuit *gMinuit = new TMinuit(npar);
    gMinuit->SetPrintLevel(-1); // -1:no , 0:summary, 1:many
    gMinuit->SetFCN(QuenchedBraggFitFunc);     
    //H set parameter   ( start, step, min, max)H//
    // gMinuit->mnparm(0, "R0"   ,    73,   0.01,     0,   100, ier);
    // gMinuit->mnparm(1, "sigma",     2, 0.001,   0.1,    20, ier);
    // gMinuit->mnparm(2, "phi0" , 1e-07, 1e-09, 1e-08, 1e-06, ier);
    // gMinuit->mnparm(3, "birks",  0.07, 0.001,  0.01,   0.5, ier);

    // gMinuit->mnparm(0, "R0"   ,    71,  0.01,     0,   100, ier);
    // gMinuit->mnparm(1, "sigma",     2, 0.001,   0.1,    20, ier);
    // gMinuit->mnparm(2, "phi0" , 1e-07, 1e-10, 1e-09, 8e-05, ier);
    // gMinuit->mnparm(3, "birks",  0.07, 0.001,  0.01,   0.5, ier);

    //good
    // gMinuit->mnparm(0, "R0"   ,    71,  0.01,     0,   100, ier);
    // gMinuit->mnparm(1, "sigma",     2, 0.001,   0.0,    10, ier);
    // gMinuit->mnparm(2, "phi0" ,10e-07,  1e-9, 1e-07, 20e-07, ier);
    // gMinuit->mnparm(3, "birks",  0.07, 0.001,  0.01,   0.5, ier);



    // gMinuit->mnparm(0, "R0"   ,    71,  0.01,     0,   100, ier);
    // gMinuit->mnparm(1, "sigma",     2, 0.001,   0.0,    20, ier);
    // gMinuit->mnparm(2, "phi0" ,    10, 0.001,     0,   100, ier);
    // gMinuit->mnparm(3, "birks",  0.07, 0.001,  0.01,   0.5, ier);


    
    gMinuit->mnparm(0, "R0"   ,      initR0,  0.01,     initR0-20,      initR0+20, ier);
    gMinuit->mnparm(1, "sigma",   initsigma, 0.001,   initsigma-1,    initsigma+1, ier);
    gMinuit->mnparm(2, "phi0" ,    initphi0, 0.001, initphi0-10e-7,  initphi0+10e-7, ier);
    gMinuit->mnparm(3, "birks",      initkb, 0.001,   initkb-0.03,    initkb+0.03, ier);

    Double_t arg[] = {2000,1.0};
    gMinuit->mnexcm("MIGRAD", arg, npar, ier);
    TString name;
    double value;
    double error;
    double other;
    double eR0, esigma, ephi0, ekb;
    gMinuit->mnpout(0, name,    R0,    eR0, other, other,  ier);
    gMinuit->mnpout(1, name, sigma, esigma, other, other,  ier);
    gMinuit->mnpout(2, name,  phi0,  ephi0, other, other,  ier);
    gMinuit->mnpout(3, name,    kb,    ekb, other, other,  ier);

    E0 = pow(R0/alpha,1/p)/mp;
    cout<<"my fit results -----------------------------------"<<endl;
    cout<<"MyF_Min R0 = "<<R0<<" +/- "<<eR0<<endl;
    cout<<"MyF_Min sigma = "<<sigma<<" +/- "<<esigma<<endl;
    cout<<"MyF_Min phi0 = "<<phi0<<" +/- "<<ephi0<<endl;
    cout<<"MyF_Min kb = "<<kb<<" +/- "<<ekb<<endl;
    cout<<"MyF_Min E0    = "<<E0<<endl;

    minuR0=R0;
    minusigma=sigma;
    minuphi0=phi0;
    minukb=kb;
    minupullR0->Fill((minuR0-pR0)/eR0);
    minupullsigma->Fill((minusigma-psigma)/esigma);
    cout<<(minuR0-pR0)/eR0<<"    minupull="<<minuR0<<"-"<<pR0<<"/"<<eR0<<endl;


    fphotons->SetParameters(R0,sigma,phi0,xOffset,kb);


    //HHHHH fw test HHHHHH//
    double ob[nbin], obe[nbin];
    for(int kk=1; kk<nbin+1; kk++){
      ob[kk-1]=hd->GetBinContent(kk);
      obe[kk-1]=hd->GetBinError(kk);
    }
    fw_binned(ob, obe, R0, sigma, phi0, kb, initR0);
    fwR0=R0;
    fwsigma=sigma;
    fwphi0=phi0;
    fwkb=kb;

    for(int l=0; l<300; l++){
      double tR0=l*(60./300.)+20;
      double lR0=GetLL(ob, obe, tR0,sigma,phi0,kb);
      LLR0->SetBinContent(l+1,lR0);

      double tsigma=l*(10./300.);
      double lsigma=GetLL(ob, obe, R0,tsigma,phi0,kb);
      LLsigma->SetBinContent(l+1,lsigma);

      double tkb=l*(0.34/300.);
      double lkb=GetLL(ob, obe, R0,sigma,phi0,tkb);
      LLkb->SetBinContent(l+1,lkb);
    }



    phi0=phi0*1e-7;
    cout<<"fw fit results -----------------------------------"<<endl;
    cout<<"MyF_Min R0 = "<<R0<<endl;
    cout<<"MyF_Min sigma = "<<sigma<<endl;
    cout<<"MyF_Min phi0 = "<<phi0<<endl;
    cout<<"MyF_Min kb = "<<kb<<endl;
    cout<<"MyF_Min E0    = "<<E0<<endl;
    fphotons->SetParameters(R0,sigma,phi0,xOffset,kb);
    //  cout<<"final integral="<<fphotons->Integral(0, 92.4508)<<endl;

    // R0 = 73.8947;
    // sigma = 2.00352;
    // phi0 = 6.12822e-07;
    // kb = 0.0751327;
    // fphotons->SetParameters(R0,sigma,phi0,xOffset,kb);

  
    // guess 
    //R0    = 73.0366
    //sigma = 0.666118
    //E0    = 95.9484
    //phi   = 1.3485e-06
    //off   = 0
    //kB    = 0.07
    // cout<<" Maybe origianl best fit ------------------------------ exmaple"<<endl;
    // cout<<"  //fit results"<<endl;
    // cout<<"  //R0 = 73.8947"<<endl;
    // cout<<"  //sigma = 2.00352"<<endl;
    // cout<<"  //phi = 6.12822e-07"<<endl;
    // cout<<"  //off = 0"<<endl;
    // cout<<"  //kB = 0.0751327"<<endl;
    // cout<<"  //E0 = 96.5939"<<endl;

    // R     peack position : 0  to 100 // peack position
    // sigma width          : 0.1 to 20
    //  phi  hight          : 1e-07 to 1e-06?   max 14 is 6 e-7   so max/2e-7 around
    // offset nocurrently
    // kB  0.01 to 0.5
    //  fphotons->SetParameters(73.9, 0.1, 6.1e-07, 0, 0.07);   //test
    
    // red
    
    //HHH make diff plots HHH//
    diforiR0->Fill(oriR0-pR0);
    difminuR0->Fill(minuR0-pR0);
    diffwR0->Fill(fwR0-pR0);
    diforisigma->Fill(orisigma-psigma);
    difminusigma->Fill(minusigma-psigma);
    diffwsigma->Fill(fwsigma-psigma);

    minuoriR0->Fill(minuR0-oriR0);
    fworiR0->Fill(fwR0-oriR0);
    minuorisigma->Fill(minusigma-orisigma);
    fworisigma->Fill(fwsigma-orisigma);
    


    //HHH  Draw  HHH// ====================================================
    TPaveText *tpt = new TPaveText(0.2,0.7,0.9,0.9,"ndc");
    tpt->SetBorderSize(0);
    tpt->SetFillColor(0);
    tpt->SetFillStyle(0);
    string text= "Input R0="+to_string(pR0,1)+", Sig="+to_string(psigma,1)+", kb="+to_string(pkb,1);
    cout<<text<<"  "<<pR0<<"  "<<psigma<<"              testsetsets"<<endl;
    tpt->AddText(text.c_str());
  
    gStyle->SetOptStat("000");
    TCanvas *c1 = new TCanvas("c1","c1",1200, 1200);
    c1->Divide(2,2);
    c1->cd(1);
    c1->SetGrid();
    c1->SetGridx();
    c1->SetGridy();  

    hd->SetMaximum(hd->GetMaximum()*1.8);
    hd->SetTitle("");
    hd->SetXTitle("WET (mm)");
    hd->SetYTitle("Charge (pC)");  
    hd->Draw();
    fphotons->Draw("same");  
    tpt->Draw();


    c1->cd(2);
    LLR0->SetTitle("");
    LLR0->SetXTitle("R0");
    LLR0->SetYTitle("-2 Log Likelihood");
    LLR0->Draw();
    c1->cd(3);
    LLsigma->SetTitle("");
    LLsigma->SetXTitle("sigma");
    LLsigma->SetYTitle("-2 Log Likelihood");
    LLsigma->Draw();
    c1->cd(4);
    LLkb->SetTitle("");
    LLkb->SetXTitle("kb");
    LLkb->SetYTitle("-2 Log Likelihood");
    LLkb->Draw();

  
  
    c1->Update();
    string out="pseudopdf/"+text+".pdf";
    c1->SaveAs(out.c_str());
  }// end of n test




TLegend* leg = new TLegend(0.10, 0.73, 0.7, 0.88);
 leg->SetBorderSize(0);
 leg->SetFillColor(0);
 leg->SetFillStyle(0);
 leg->AddEntry(diforiR0, "OrigFit");
 leg->AddEntry(difminuR0, "FW+Minui");
 leg->AddEntry(diffwR0, "Full FW"); 
  TCanvas *cdiff = new TCanvas("diff","diff",1800, 1200);
  cdiff->Divide(3,2);
  cdiff->cd(1);
  difminuR0->SetMaximum(difminuR0->GetMaximum()*1.5);
  difminuR0->SetTitle("Dif R0 vs true");
  difminuR0->SetXTitle("R0 Fit-True");
  difminuR0->SetYTitle("Fits");
  diforiR0->SetLineColor(1);
  difminuR0->SetLineColor(2);
  diffwR0->SetLineColor(4);
  difminuR0->Draw();
  diforiR0->Draw("same");
  diffwR0->Draw("same");
  leg->Draw("same");
  
  cdiff->cd(2);
  difminusigma->SetMaximum(diforisigma->GetMaximum()*1.5);
  difminusigma->SetTitle("Dif sigma vs true");
  difminusigma->SetXTitle("sigma Fit-True");
  difminusigma->SetYTitle("Fits");
  diforisigma->SetLineColor(1);
  difminusigma->SetLineColor(2);
  diffwsigma->SetLineColor(4);
  difminusigma->Draw();
  diforisigma->Draw("same");
  diffwsigma->Draw("same");
  leg->Draw("same");

  cdiff->cd(4);
  minuoriR0->SetMaximum(minuoriR0->GetMaximum()*1.5);

  minuoriR0->SetTitle("Dif R0 vs orig Fit");
  minuoriR0->SetXTitle("R0 NewFit-OriginalFit");
  minuoriR0->SetYTitle("Fits");
  minuoriR0->SetLineColor(2);
  fworiR0->SetLineColor(4);
  minuoriR0->Draw();
  fworiR0->Draw("same");
  
  cdiff->cd(5);
  minuorisigma->SetMaximum(minuorisigma->GetMaximum()*1.5);
  minuorisigma->SetTitle("Dif sigma orig Fit");
  minuorisigma->SetXTitle("sigma NewFit-OriginalFit");
  minuorisigma->SetYTitle("Fits");
  minuorisigma->SetLineColor(2);
  fworisigma->SetLineColor(4);
  minuorisigma->Draw();
  fworisigma->Draw("same");

  cdiff->cd(3);
  gStyle->SetOptFit(111);
  


  minupullR0->Draw();
  minupullR0->Fit("gaus");

  cdiff->cd(6);
  minupullsigma->Draw();
  minupullsigma->Fit("gaus");


  cdiff->SaveAs("test.pdf");


  return 0;
}