// standard c++ headers                                                                     

#ifndef __CINT__
#include "RooNumConvPdf.h"
#include "TPaveLabel.h"
#include <iostream>
#include <string>
#include <vector>
#include <math.h>
#include "TROOT.h"
#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGlobalFunc.h"
#include "RooChi2Var.h"
#include "RooMCStudy.h"
#endif
#include "RooHist.h"
#include "RooPlot.h"
#include "TColor.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooCBShape.h"
#include "RooConstVar.h"
#include "RooDataHist.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooWorkspace.h"
#include "RooFitResult.h"
#include "RooExtendPdf.h"
#include "RooPlot.h"
#include "TString.h"
#include "TCanvas.h"
#include <sstream>
#include "THStack.h"
#include "TH1.h"
#include "TF1.h"
#include "TTree.h"
#include <TList.h>
#include "TFile.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include <TStyle.h>
#include "TLatex.h"
#include "TGraph.h"
#include "TLimit.h"
#include "TLine.h"
#include "TVector.h"
#include "TMath.h"
#include "RooBreitWigner.h"
#include "RooVoigtian.h"
#include "RooFFTConvPdf.h"
#include "TPaveText.h"
#include "RooChi2MCSModule.h"
using namespace std;
using namespace RooFit;

// Global Variables                                                                                                       
Float_t lumi=14324.0;

// Method to vit a Voigtian function for W peak 
void fitVoigforW(TString filename="", TString fitplotname="", TString pullplotname="", TString residplotname="")
{
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(111110110);
  //------------------------------
  // get the file and histo(names to be used)
  //------------------------------

  //TFile *f = new TFile(filename);
  TFile *f = new TFile("/export/share/data2/jveatch/testarea/TopAnalysis/rel18/d5pd/v03/D5PDMacros/JMR.QCDUp100.root");
  if(!f) std::cout<<"No file found"<<std::endl;
  TH1D *histo = (TH1D*) f->Get("data_sig_LeadFJetM2p5GeV_cakt10_fjW_200in"); 
  if(!histo) std::cout<<"No histo found"<<std::endl;
   
  //------------------------------------
  // define x range and frame
  //------------------------------------
  
  RooRealVar x("W Mass","Mass",0.,150.,"GeV/c^{2}") ;
  RooPlot* frame = x.frame(Title("Voigtian Fit for M_{W}")) ;
  TCanvas *c1 = new TCanvas("c1","c1",25,71,400,400);


  //------------------------------------
  // import histo  
  //------------------------------------
  
  RooDataHist *datah;
  datah = new RooDataHist("datah","datah",x,histo);
  datah->plotOn(frame,Name("datah"),
		DataError(RooAbsData::SumW2),MarkerColor(kRed),MarkerSize(0.8),XErrorSize(0.));
  
  //------------------------------------
  // define Voigitian Fitting Function
  //------------------------------------
  // As per truth jets 
//  Floating Parameter    FinalValue +/-  Error
//  --------------------  --------------------------
//                  mean    8.0574e+01 +/-  2.09e-02
//               nsignal    5.0467e+03 +/-  2.44e+01
//                 sigma    7.1535e+00 +/-  5.57e-02
//


  // Construct the voigtian pdf by defining parameters. 
  
//  RooRealVar mean ("mean", "mean", 80.57,75,90);
//  RooRealVar width ("width", "width", 2.0,1.8,2.2);
//  RooRealVar sigma ("sigma", "sigma", 2,0,10);
//  RooVoigtian gauss("gauss", "gauss", x, mean, width, sigma);
//  

  RooRealVar mean ("mean", "mean", 80.57,75,90);
  RooRealVar width ("width", "width", 7.153);
  RooRealVar sigma ("sigma", "sigma", 2,0,10);
  RooVoigtian gauss("gauss", "gauss", x, mean, width, sigma);


  // fix the number of signal events in the histogram 
  RooRealVar nsig("nsignal","nsignal",3000,2000,10000); // MC 3400 and Data 4200 

  //-------------------------------------------
  // Add the pdfs 
  //-----------------------------------------
  RooFitResult* r;
  x.setRange("fitrange",70.,100.);

  RooExtendPdf esig("esig","esig",gauss,nsig,"fitrange");
  RooAddPdf sum("sum","Vogitian fit",RooArgList(esig));


  r = sum.fitTo(*datah,Range("fitrange"),Extended(kTRUE),SumW2Error(kTRUE),Save()) ;

 //------------------------------------
 // plot the fit and Data 
 //------------------------------------  
  //  sum.plotOn(frame,Components("gauss"),LineColor(kGreen),LineStyle(kDashed),Name("mypdf"));
  sum.plotOn(frame,Name("mypdf"));
  sum.paramOn(frame,Layout(0.15,0.45,0.875));
  frame->getAttText()->SetTextSize(0.03);
  
  Int_t floated = static_cast<Int_t>(r->floatParsFinal().getSize());
    
  // Adding some chi2 values 
  RooChi2Var chi2 ("chi2","chi2",sum,*datah);
  Double_t chi2_val = chi2.getVal(); 
  Double_t finalchi2_val=chi2_val/(20.0);
  
  TPaveLabel *t1=new TPaveLabel(0.15,0.53,0.45,0.62, Form("#chi^{2}/dof = %f",finalchi2_val),"brNDC");
  std::cout<<"This is the final chi2 :"<<finalchi2_val<<std::endl;
  
  
  std::cout<<"Let us check another value of chi2 :"<<frame->chiSquare("mypdf","datah",floated)<<std::endl; 

  t1->SetFillColor(0);
  t1->SetTextSize(0.3);
  frame->addObject(t1);
  frame->Draw();
  
  c1->SaveAs(fitplotname);
  //------------------------------------------------
  // Also print out the parameters 
  //-------------------------------------------------
  r->Print();
  r->correlationMatrix().Print() ; 
  
  
  // Construct a histogroam with pullHist 
  RooHist* hpull=frame->pullHist();
  RooHist* hresid=frame->residHist();
  // Create a new frame to draw the residual distribution and add the distribution to the frame
  RooPlot* frame2 = x.frame(Title("Residual Distribution")) ;
  frame2->addPlotable(hresid,"P");
  
  // Create a new frame to draw the pull distribution and add the distribution to the frame
  RooPlot* frame3 = x.frame(Title("Pull Distribution")) ;
  frame3->addPlotable(hpull,"P") ;

  // REsid 
  TCanvas* c = new TCanvas("c","c",25,71,400,400);    
  frame2->Draw() ;
  c->SaveAs(residplotname);

  //Pull 
  TCanvas* c2 = new TCanvas("c2","c2",25,71,400,400) ;
  frame3->Draw() ;
  c2->SaveAs(pullplotname);
  
  //  c->Divide(2) ;
  //  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;
  //  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; 
  //  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; 
  //  c->SaveAs(pullplotname);
  
}// end of voigt fit method for W peak 


void fitGaussforW(TString filename="", TString fitplotname="", TString pullplotname="")
{
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(111110110);
  //------------------------------
  // get the file and histo(names to be used)
  //------------------------------
  TFile *f = new TFile("/export/home/jveatch/testarea/TopAnalysis/rel18/d5pd/v03/D5PDMacros/JMR.2015.04.01/JMR.2015.04.01.NOM.root");
  //  TFile *f = new TFile(filename);
  if(!f) std::cout<<"No file found"<<std::endl;
  TH1D *histo = (TH1D*) f->Get("mc_sig_LeadFJetM2p5GeV_cakt10_fjW_200in"); 
  if(!histo) std::cout<<"No histo found"<<std::endl;
   
  //------------------------------------
  // define x range and frame
  //------------------------------------
  
  RooRealVar x("W Mass","Mass",0.,150.,"GeV/c^{2}") ;
  RooPlot* frame = x.frame(Title("Conv. Gaussian Fit for M_{W}")) ;
  TCanvas *c1 = new TCanvas("c1","c1",25,71,400,400);
  
  //------------------------------------
  // import histo  
  //------------------------------------
  
  RooDataHist *datah;
  datah = new RooDataHist("datah","datah",x,histo);
  datah->plotOn(frame,Name("datah"),
		DataError(RooAbsData::SumW2),MarkerColor(kRed),MarkerSize(0.8),XErrorSize(0.));
  
  //------------------------------------
  // define Voigitian Fitting Function
  //------------------------------------

  // Construct the voigtian pdf by defining parameters. 
  
  RooRealVar mean ("mean", "mean", 80.3);//,70,90);
  RooRealVar sigma ("sigma", "sigma", 4.551);
  
  RooRealVar mean1 ("mean1", "mean1", 80.3,70,90);
  RooRealVar sigma1 ("sigma1", "sigma1", 2.0,0,10);
  RooGaussian gauss("gauss", "gauss", x, mean, sigma);
  RooGaussian gauss1("gauss1", "gauss1", x, mean1, sigma1);
  
  // fix the number of signal events in the histogram 
  RooRealVar nsig("nsignal","nsignal",3000,2000,10000); // MC 3400 and Data 4200 

  //-------------------------------------------
  // Add the pdfs 
  //-----------------------------------------
  RooFitResult* r;
  x.setRange("fitrange",60.,120.);

  RooExtendPdf esig("esig","esig",gauss,nsig,"fitrange");
  RooExtendPdf esig1("esig1","esig1",gauss1,nsig,"fitrange");


  //  RooAddPdf sum("sum","Gaussian fit",RooArgList(esig,esig1));
  //  RooNumConvPdf sum("sum","sum",x,gauss,gauss1);

  RooFFTConvPdf sum("sum","gauss X gauss1",x,esig,esig1);
  r = sum.fitTo(*datah,Range("fitrange"),Extended(kTRUE),SumW2Error(kTRUE),Save()) ;

 //------------------------------------
 // plot the fit and Data 
 //------------------------------------  
  sum.plotOn(frame);
  //  sum.plotOn(frame, Components("bw1"), LineStyle(kDashed));
  
  sum.paramOn(frame,Layout(0.15,0.45,0.875));
  frame->getAttText()->SetTextSize(0.03);
  
  Int_t floated = static_cast<Int_t>(r->floatParsFinal().getSize());
  std::cout<<"These are the # of final Pars :"<<floated<<std::endl;
    
  // Adding some chi2 values 
  RooChi2Var chi2 ("chi2","chi2",sum,*datah);
  Double_t chi2_val = chi2.getVal(); 
  Double_t finalchi2_val=chi2_val/(20.0);
  
  TPaveLabel *t1=new TPaveLabel(0.15,0.53,0.45,0.62, Form("#chi^{2}/dof = %f",finalchi2_val),"brNDC");
  std::cout<<"This is the final chi2 :"<<finalchi2_val<<std::endl;
  t1->SetFillColor(0);
  t1->SetTextSize(0.3);
  frame->addObject(t1);
  frame->Draw();
  
  c1->SaveAs(fitplotname);
  //------------------------------------------------
  // Also print out the parameters 
  //-------------------------------------------------
  r->Print();
  r->correlationMatrix().Print() ; 
  
  
  // Construct a histogroam with pullHist 
  RooHist* hpull=frame->pullHist();
  RooHist* hresid=frame->residHist();
  // Create a new frame to draw the residual distribution and add the distribution to the frame
  RooPlot* frame2 = x.frame(Title("Residual Distribution")) ;
  frame2->addPlotable(hresid,"P");
  
  // Create a new frame to draw the pull distribution and add the distribution to the frame
  RooPlot* frame3 = x.frame(Title("Pull Distribution")) ;
  frame3->addPlotable(hpull,"P") ;
  
  TCanvas* c = new TCanvas("rf109_chi2residpull","rf109_chi2residpull",900,300) ;
  c->Divide(2) ;
  //  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; frame2->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; frame3->Draw() ;
  c->SaveAs(pullplotname);
  
}// end of gaussian fit method for W peak 


void fitBWforW(TString filename="", TString fitplotname="", TString pullplotname="")
{
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(111110110);
  //------------------------------
  // get the file and histo(names to be used)
  //------------------------------
  TFile *f = new TFile("/export/home/jveatch/testarea/TopAnalysis/rel18/d5pd/v03/D5PDMacros/JMR.2015.03.18/JMR.2015.03.18.NOM.root");
  //  TFile *f = new TFile(filename);
  if(!f) std::cout<<"No file found"<<std::endl;
  TH1D *histo = (TH1D*) f->Get("mc_sig_LeadFJetM2p5GeV_takt10_fjW_200in"); 
  if(!histo) std::cout<<"No histo found"<<std::endl;
   
  //------------------------------------
  // define x range and frame
  //------------------------------------
  
  RooRealVar x("W Mass","Mass",0.,150.,"GeV/c^{2}") ;
  RooPlot* frame = x.frame(Title("BW Fit for M_{W}")) ;
  TCanvas *c1 = new TCanvas("c1","c1",25,71,400,400);


  
  //------------------------------------
  // import histo  
  //------------------------------------
  
  RooDataHist *datah;
  datah = new RooDataHist("datah","datah",x,histo);
  datah->plotOn(frame,Name("datah"),
		DataError(RooAbsData::SumW2),MarkerColor(kRed),MarkerSize(0.8),XErrorSize(0.));
  
  //------------------------------------
  // define Voigitian Fitting Function
  //------------------------------------

  // Construct the voigtian pdf by defining parameters. 
  
  RooRealVar mean ("mean", "mean", 80.3,70,100);
  RooRealVar sigma ("sigma", "sigma", 2,0,10);
  RooBreitWigner bw("bw","A Breit-Wigner Distribution",x,mean,sigma);

  // fix the number of signal events in the histogram 
  RooRealVar nsig("nsignal","nsignal",3000,2000,10000); // MC 3400 and Data 4200 

  //-------------------------------------------
  // Add the pdfs 
  //-----------------------------------------
  RooFitResult* r;
  x.setRange("fitrange",70.,90.);

  RooExtendPdf esig("esig","esig",bw,nsig,"fitrange");
  RooAddPdf sum("sum","BW fit",RooArgList(esig));

  r = sum.fitTo(*datah,Range("fitrange"),Extended(kTRUE),SumW2Error(kTRUE),Save()) ;

 //------------------------------------
 // plot the fit and Data 
 //------------------------------------  
  sum.plotOn(frame,Components("bw"),LineColor(kGreen),LineStyle(kDashed));
  sum.paramOn(frame,Layout(0.15,0.45,0.875));
  frame->getAttText()->SetTextSize(0.03);
  
  Int_t floated = static_cast<Int_t>(r->floatParsFinal().getSize());
  std::cout<<"These are the # of final Pars :"<<floated<<std::endl;
  
  
  // Adding some chi2 values 
  RooChi2Var chi2 ("chi2","chi2",sum,*datah);
  Double_t chi2_val = chi2.getVal(); 
Double_t finalchi2_val=chi2_val/(21.0);
  
  TPaveLabel *t1=new TPaveLabel(0.15,0.53,0.45,0.62, Form("#chi^{2}/dof = %f",finalchi2_val),"brNDC");
  std::cout<<"This is the final chi2 :"<<finalchi2_val<<std::endl;
  t1->SetFillColor(0);
  t1->SetTextSize(0.3);
  frame->addObject(t1);
  frame->Draw();
  
  c1->SaveAs(fitplotname);
  //------------------------------------------------
  // Also print out the parameters 
  //-------------------------------------------------
  r->Print();
  r->correlationMatrix().Print() ; 
  
  
  // Construct a histogroam with pullHist 
  RooHist* hpull=frame->pullHist();
  RooHist* hresid=frame->residHist();
  // Create a new frame to draw the residual distribution and add the distribution to the frame
  RooPlot* frame2 = x.frame(Title("Residual Distribution")) ;
  frame2->addPlotable(hresid,"P");
  
  // Create a new frame to draw the pull distribution and add the distribution to the frame
  RooPlot* frame3 = x.frame(Title("Pull Distribution")) ;
  frame3->addPlotable(hpull,"P") ;
  
  TCanvas* c = new TCanvas("rf109_chi2residpull","rf109_chi2residpull",900,300) ;
  c->Divide(2) ;
  //  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; frame2->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; frame3->Draw() ;
  c->SaveAs(pullplotname);
  
}// end of gaussian fit method for W peak 


// Crystal Ball fitting function for Top peak 
void fit2CrystalBall(TString filename="", TString fitplotname="fitTop_qcdup100.eps", TString pullplotname="pullTop_qcdup100.eps", TString residplotname="residTop_qcdup100.eps")
{
  
  gROOT->SetStyle("Plain");
  gStyle->SetOptFit(1111);
  // get the file and histo you wish to fit 
  TFile *f = new TFile("/export/share/data2/jveatch/testarea/TopAnalysis/rel18/d5pd/v03/D5PDMacros/JMR.QCDCorr.root");
  //  TFile *f = new TFile(filename);
  if(!f) std::cout<<"No file found "<<std::endl;
  TH1D *histo = (TH1D*) f->Get("mc_sig_LeadFJetM2p5GeV_cakt10_fjTop_300in"); 

  if(!histo) std::cout<<"No histo found"<<std::endl;

 //------------------------------------
 // define x range and frame
 //------------------------------------

 RooRealVar x("Top Mass","Mass",100,250.,"GeV/c^{2}") ;
 RooPlot* frame = x.frame(Title("Crystall Ball fit for M_{t}")) ;
 TCanvas *c1 = new TCanvas("c1","c1",25,71,400,400);

 //------------------------------------
 // import histo  
 //------------------------------------

 RooDataHist *datah;
 datah = new RooDataHist("datah","datah",x,histo);
 datah->plotOn(frame,Name("datah"),
 	       DataError(RooAbsData::SumW2),MarkerColor(kRed),MarkerSize(0.8),XErrorSize(0.));
 
 //------------------------------------
 // define Crystal Ball
 //------------------------------------

 RooRealVar cbmean("cb mean","cb_mean",173.2,165,175);
 RooRealVar cbwidth("cb width","cb_width",10,0,15);
 RooRealVar cbn("cb n","cb_n",150,100,180);
 RooRealVar cbalpha("cb alpha","cb_alpha",5,0,10.);
 RooCBShape cb("cb","cb",x,cbmean,cbwidth,cbalpha,cbn) ;
 
 // cbalpha.setVal(10);   
 // cbalpha.setConstant(kTRUE);
 cbn.setVal(5); 
 cbn.setConstant(kTRUE);
 
 RooRealVar nsig("N.top","nsignal",2500,2000,3000);
 RooExtendPdf esig("esig","esig",cb,nsig);

 // Final pdf that I use: 
 RooAddPdf sum("cbfit","crystal ball fit",RooArgList(esig));

 RooFitResult* r;
 x.setRange("fitrange",140.,200.);
 r = sum.fitTo(*datah,Range("fitrange"),Extended(kTRUE),SumW2Error(kTRUE),Save()) ;
 sum.paramOn(frame,Layout(0.1,0.4,0.875));
 frame->getAttText()->SetTextSize(0.03);

 Int_t floated = static_cast<Int_t>(r->floatParsFinal().getSize());

 sum.chi2FitTo(*datah);
 // Adding some chi2 values 
 RooChi2Var chi2 ("chi2", "chi2", sum, *datah);
 Double_t chi2_val = chi2.getVal(); 
 Double_t finalchi2_val=chi2_val/(20.0);
 std::cout<<" This is the final value of chi2/DOF :"<<finalchi2_val<<std::endl;

 // Put the Chi2 label on the plot 


 //------------------------------------
 // plot
 //------------------------------------  
 sum.plotOn(frame,Name("mypdf")) ;
 frame->Draw() ;
 c1->SaveAs(fitplotname);


 r->Print();
 r->correlationMatrix().Print() ; 
 double a=frame->chiSquare("mypdf","datah",floated);
 std::cout<<"Let us check another value of chi2 :"<<a<<std::endl; 

 TPaveLabel *t1=new TPaveLabel(0.15,0.53,0.35,0.62, Form("#chi^{2}/dof = %f",a),"brNDC");
 t1->SetFillColor(0);
 t1->SetTextSize(0.3);
 frame->addObject(t1);

 
 // Construct a histogroam with pullHist 
 RooHist* hpull=frame->pullHist();
 RooHist* hresid=frame->residHist();

 // Create a new frame to draw the residual distribution and add the distribution to the frame
 RooPlot* frame2 = x.frame(Title("Residual Distribution")) ;
 frame2->addPlotable(hresid,"P");
 
 // Create a new frame to draw the pull distribution and add the distribution to the frame
 RooPlot* frame3 = x.frame(Title("Pull Distribution")) ;
 frame3->addPlotable(hpull,"P") ;
 

 TCanvas* c = new TCanvas("pull","pull",25,71,400,400);
 frame3->Draw() ;c->SaveAs(pullplotname);

 TCanvas* c2 = new TCanvas("resid","resid",25,71,400,400);
 frame2->Draw() ;c2->SaveAs(residplotname);


 // ************ Create manager for MC Study ******************
 RooMCStudy* mcstudy =new RooMCStudy(sum,x,Binned(true),Verbose(false),Extended(),Silence() );

 // Add a chi2 module
 RooChi2MCSModule chi2mod ;
 mcstudy->addModule(chi2mod) ;

 // Generate and fit results 
 mcstudy->generateAndFit(1000,2400);

 
 // Fill histograms with distributions chi2 and prob(chi2,ndf) that are calculated by RooChiMCSModule
 TH1* hist_chi2 = mcstudy->fitParDataSet().createHistogram("chi2") ; 
 TH1* hist_chi2red = mcstudy->fitParDataSet().createHistogram("chi2red") ; 
 TH1* hist_prob = mcstudy->fitParDataSet().createHistogram("prob") ;

 TCanvas* c3 = new TCanvas("mcstudy_addons","rf802_mcstudy_addons",800,400) ;
 c3->Divide(2);
 c3->cd(1) ; gPad->SetLeftMargin(0.15) ; hist_chi2->GetYaxis()->SetTitleOffset(1.4) ; hist_chi2->Draw() ;
 c3->cd(2) ; gPad->SetLeftMargin(0.15) ; hist_chi2red->GetYaxis()->SetTitleOffset(1.4) ; hist_chi2red->Draw() ;
 c3->Draw();


 TCanvas* c4 = new TCanvas("mcstudy_addons_2","rf802_mcstudy_addons_2",800,400) ;
 RooPlot* mframe= mcstudy->plotParam(cbmean,Bins(40));
 mframe->Draw();
 c4->Draw();
}//end of method for CrytalBall Fitting for top 


//main method  used to call the fitting functions 
void fitMassPeak(TString filename="",TString fitplotname="", TString pullplotname="", TString residplotname="")
{
  // setting gstyle options.......                                                                                                 
  gStyle->SetPalette(1);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetOptStat(1111111111);
  gStyle->SetOptFit(111);
  gROOT->SetStyle("Plain");
  
  std::cout<<"Going to fit ():"<<std::endl;
  fit2CrystalBall(filename,fitplotname,pullplotname,residplotname);
  //fitVoigforW(filename,fitplotname,pullplotname,residplotname);
  //fitGaussforW(filename,fitplotname,pullplotname);
  //  fitBWforW(filename,fitplotname,pullplotname);
}// end of main method