////////////////////////////////////////////////////////////////////////
/// $Id: SkzpFitter.cxx,v 1.1 2008/02/12 19:22:15 rearmstr Exp $
///
/// (Document me!)
///
/// rearmstr@indiana.edu
////////////////////////////////////////////////////////////////////////
#include <fstream>
#include <iostream>
#include <iomanip>

#include "SkzpFitter.h"
#include "Parsers.h"
#include "GenericFactory.h"
#include "FitData.h"
#include "ZbeamWeight.h"
#include "ZflukWeight.h"
#include "ZdetWeight.h"


#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "MCReweight/NeugenWeightCalculator.h"
#include "MCReweight/Zbeam.h"
#include "MCReweight/Zfluk.h"

#include "TStopwatch.h"
#include "TFile.h"
#include "TMatrix.h"
#include "TMatrixD.h"
#include "TH2D.h"
#include "TDirectoryFile.h"
#include "TGraph.h"
#include "TTree.h"

#include <boost/regex.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string.hpp>

using namespace boost;

JOBMODULE(SkzpFitter, "SkzpFitter",
          "Uses the Skzp Parameterization to weight the MC");
CVSID("");
//......................................................................

// global variable to use in fit function
SkzpFitter *gFitter;
void fit_fcn(Int_t &, double *, double &f, double *par, int);

SkzpFitter::SkzpFitter() 
  :  fInputFile(""), 
     fOutputFile(""),
     fFitter("Minuit"),
     fIterations(0),
     fMaxIterations(10000),
     fTolerance(10.),
     fStrategy(1),
     fERRDEF(1.),					
     fNParFree(0),
     fNParFixed(0),
     fNPar(0),
     fNDF(0),
     fScanPoints(50),
     fUseNeugen(true),
     fNeugenConfigName("MODBYRS"),
     fNeugenVer(4),
     fMaQE(1.),
     fMaRES(1.),
     fKNO(0.),
     fUseOldKNO(false),
     fUseSmoothing(true),
     fSmoothFunc("(1-x**2)**2"),
     fFuncMin(-1.),
     fFuncMax(1.),
     fNormMin(-1.),
     fNormMax(1.),
     fFuncBins(1001),
     fFuncNorm(true),
     fFitParFile("")
{
  ///
  /// (Document me!)
  ///
  gFitter=this;
  mcr = &MCReweight::Instance();
  fSkzpCalc=&SkzpCalculator::Instance();
  smooth=new Smooth();
  nwc_reg = new Registry();
  fPar_reg.SetName("Par");
  fChi2_reg.SetName("Chi2");

}
//......................................................................

SkzpFitter::~SkzpFitter()
{
  ///
  /// (Document me!)
  ///
  
  delete smooth;
  delete nwc_reg;

}

//......................................................................

void SkzpFitter::BeginJob()
{
  ///
  /// (Document me!)
  ///
  bool success;
  // Read the Configuration file;
  success=ReadInputFile(fInputFile);
  if(!success) {
    MSG("SkzpFitter",Msg::kInfo)<<" Error Reading input file "
				<<fInputFile<<endl;
    return;
  }
  success=Setup();
  if(!success) {
    MSG("SkzpFitter",Msg::kInfo)<<" Error in Initialization"<<endl;
    return;
  }

  
}

//......................................................................

void SkzpFitter::EndJob()
{
  ///
  /// (Document me!)
 ///
  Write();
  fOut->Close();
  //TFile *file=

}

//......................................................................

JobCResult SkzpFitter::Ana(const MomNavigator* )
{
  ///
  /// (Document me!)
  ///


  std::vector<std::string>::iterator com_iter;
  arglist[0]=fMaxIterations;
  arglist[1]=fTolerance;
  
  int errflag=0;

  for(com_iter=fCommands.begin();
      com_iter!=fCommands.end(); ++com_iter) {
    
    std::string cmd=*com_iter;

    MSG("SkzpFitter",Msg::kInfo)<<" Starting minimization routine " <<cmd<<endl;
    
    if( cmd=="MINOS" || cmd=="IMPROVE" || cmd=="HESSE") {
      
      errflag=fitter->ExecuteCommand(cmd.c_str(),arglist,1);
      MSG("SkzpFitter",Msg::kInfo)<<cmd<<" finished with flag "<<errflag<<endl;
    }
    else if( cmd=="MIGRAD" || cmd=="SIMPLEX") {
      
      errflag=fitter->ExecuteCommand(cmd.c_str(),arglist,2);
      MSG("SkzpFitter",Msg::kInfo)<<cmd<<" finished with flag "<<errflag<<endl;
    }
    else if( cmd.find("SCAN")!=std::string::npos) { 
      Scan(cmd);
     }
    
    else {
      MSG("SkzpFitter",Msg::kWarning) << cmd << " is not a valid command"<<endl;
    }
  }
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

JobCResult SkzpFitter::Get(MomNavigator* )
{
  ///
  /// (Document me!)
  ///
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

JobCResult SkzpFitter::Put(const MomNavigator* )
{
  ///
  /// (Document me!)
  ///
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

const Registry& SkzpFitter::DefaultConfig() const
{
  ///
  /// Supply the default configuration for the module
  ///
  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.Set("InputFile",       "fit.inp");
  r.Set("OutputFile",       "fit.root");

  return r;
}

//......................................................................

void SkzpFitter::Config(const Registry& r)
{
  ///
  /// Configure the module given the Registry r
  ///
  fInputFile = r.GetCharString("InputFile");
  fOutputFile = r.GetCharString("OutputFile");

}

//......................................................................

void SkzpFitter::Help()
{
  ///
  /// (Document me!)
  ///
}

//......................................................................

void SkzpFitter::HandleCommand(JobCommand* )
{
  ///
  /// (Document me!)
  ///
}

//......................................................................

bool SkzpFitter::ReadInputFile(const std::string file)
{
  ///
  /// (Document me!)
  ///

  std::ifstream stream;
  stream.open(file.c_str());
  if(!stream.is_open()) {
    MSG("SkzpFitter",Msg::kError)<<"Can't Open Input file "<<file<<endl;
    return false;
  }
  
  MSG("SkzpFitter",Msg::kInfo)<<"Reading beam data from "<<file<<endl;
  
  std::string sline;
  
  while (getline(stream,sline)) {
   
      if (sline.find("*",0)==0 || 
	  sline.size()==0 ||
	  sline.find("//",0)==0) continue;  
      
      std::vector<std::string> parse;
      split(parse, sline, is_any_of(",{} "), token_compress_on);

      // if the last entry or first is blank delete it
      if(parse.size()>1) {
	if(parse[parse.size()-1]=="") parse.pop_back();
	if(parse[0]=="") parse.erase(parse.begin());
      }
      
      
      if (parse[0]=="NEUGEN") {
	fUseNeugen=true;
	if (parse[1]=="false") {
	  fUseNeugen=false;
	  fNeugenConfigName=parse[2];
	  fNeugenVer=atoi(parse[3].c_str());
	}
 	else {
 	  fNeugenConfigName=parse[1];
 	  fNeugenVer=atoi(parse[2].c_str());
 	}
      }
      else if (parse[0]=="XSECRW") {

	if(parse[1]=="true") {
	  {
	    fUseNeugen=true;
	    fMaQE=atof(parse[2].c_str());
	    fMaRES=atof(parse[3].c_str());
	    fKNO=atof(parse[4].c_str());
	  }
	  continue;
	}
      }
      else if (parse[0]=="USENEWKNO") {
	fUseOldKNO=false;
      }
      else if (parse[0]=="MINUIT") {

	// delete the first entry which is MINUIT;
	parse.erase(parse.begin());
	fCommands=parse;
      }
      else if (parse[0]=="STRATEGY") {
	fStrategy=atoi(parse[1].c_str());
      }
      else if (parse[0]=="FITTER") {
	fFitter=parse[1];
      }
      else if (parse[0]=="MAX_ITERS") {
	fMaxIterations=atoi(parse[1].c_str());
      }
      else if (parse[0]=="TOL") {
	fTolerance=atof(parse[1].c_str());
      }
      else if (parse[0]=="SMOOTHING") {
	Parsers::parse_val(parse[1],fUseSmoothing);
      }
      else if (parse[0]=="SMOOTH_FUNC") {
	fSmoothFunc=parse[1];
      }
      else if (parse[0]=="FUNC_MIN") {
	fFuncMin=atof(parse[1].c_str());
      }
      else if (parse[0]=="FUNC_MAX") {
	fFuncMax=atof(parse[1].c_str());
      }
      else if (parse[0]=="NORM_MIN") {
	fNormMin=atof(parse[1].c_str());
      } 
      else if (parse[0]=="NORM_MAX") {
	fNormMax=atof(parse[1].c_str());
      } 
      else if (parse[0]=="FUNC_BINS") {
	fFuncBins=atoi(parse[1].c_str());
      } 
      else if (parse[0]=="FITPARFILE") {
	fFitParFile=parse[1];
      } 
      else if (parse[0]=="CHI2FUNCS") {

	// Check for a multiline definition
	sline=boost::trim_right_copy(sline);
	if(sline.substr(sline.size()-1)=="{") {
	  string val("{");
	  do {
	    getline(stream,sline);
	    sline=trim_copy(sline);
	    val+=sline;
	  } while (sline!="}");
	  
	  //	  split(fChi2FuncNames, val, is_any_of(",{} "), token_compress_on);
	  //if(fChi2FuncNames[parse.size()-1]=="" || 
	  // fChi2FuncNames[parse.size()-1]==" ") fChi2FuncNames.pop_back();
	  Parsers::parse_val(val,fChi2FuncNames);
	  
	  
	}
	else {
	  parse.erase(parse.begin());
	  fChi2FuncNames=parse;
	}
	
      }
      else if (parse[0]=="WEIGHTS") {
	parse.erase(parse.begin());
	fWeightNames=parse;
      }
      // if the first word of the line is not one of the
      // above assume you are setting the properties of
      // your weight or chi2func classes
      // these lines should be of the form of
      // Name:Property=Value,Value2
      // can be multiple on one line
      else {
	string obj=sline.substr(0,sline.find(":"));
	string properties=sline.substr(sline.find(":")+1);

	properties=boost::trim_right_copy(properties);
 	if(properties.substr(properties.size()-1)=="{") {
	  
	  do {
	    getline(stream,sline);
	    properties+=sline;
	    sline=trim_copy(sline);
	  } while (sline!="}");
	  MSG("SkzpFitter",Msg::kVerbose)<<"Adding "<<properties<<" to "<<obj<<endl;
	  fProperties[obj].push_back(properties);
	}
	else {
	  std::vector<std::string> props;
	  split(props, properties, is_any_of(", "), token_compress_on);
	  std::vector<std::string>::iterator iter=props.begin();
	  std::vector<std::string>::iterator iter_end=props.end();
	  for( ; iter!=iter_end ; ++iter) {
	    MSG("SkzpFitter",Msg::kVerbose)<<"Adding "<<(*iter)<<" to "<<obj<<endl;
	    fProperties[obj].push_back(*iter);
	    
	  }
	}
	
      }
      
  }
  
  stream.close();
  return true;
  
}

//......................................................................

bool SkzpFitter::Setup() 
{

  bool success=true;
  // setup generator reweighting
  if (fUseNeugen) {

    NeugenWeightCalculator *nwc=new NeugenWeightCalculator();

    
    nwc_reg->UnLockValues();
    nwc_reg->UnLockKeys();
    nwc_reg->Clear();

    nwc_reg->Set("neugen:config_name",fNeugenConfigName.c_str());
    nwc_reg->Set("neugen:config_no",fNeugenVer);

    // Set the standard to have the configuration
    //nwc->SetStandardConfig(nwc_reg);
    

    nwc_reg->Set("neugen:use_scale_factors",1);
    nwc_reg->Set("neugen:ma_qe",fMaQE);
    nwc_reg->Set("neugen:ma_res",fMaRES); 
    if(fUseOldKNO) {
      nwc_reg->Set("neugen:kno_r112",(0.1+fKNO*0.1)/(0.1));
      nwc_reg->Set("neugen:kno_r122",(0.3+fKNO*0.1)/(0.3));
    }
    else {
      
    // This is specific to MODBYRS v4 where we weight them 
    // by their 1 sigma error as given in docdb-2989
    // we need to find the ratio of the one sigma errors to the
    // nominal value so they must be hard-wired in.
    
    double kno1_val=0.1; //for ri12 and ri42
    double kno2_val=0.3; //for ri22 and ri32
    double kno3_val=1.0; //for rij3 
    
    double kno1_err=0.1; //for rij2 
    double kno2_err=0.2; //for rij3

    nwc_reg->Set("neugen:kno_r112",
			 (kno1_val+fKNO*kno1_err)/(kno1_val));
    nwc_reg->Set("neugen:kno_r142",
			 (kno1_val+fKNO*kno1_err)/(kno1_val));
    nwc_reg->Set("neugen:kno_r122",
			 (kno2_val+fKNO*kno1_err)/(kno2_val));
    nwc_reg->Set("neugen:kno_r132",
			 (kno2_val+fKNO*kno1_err)/(kno2_val));
    nwc_reg->Set("neugen:kno_r113",
			 (kno3_val+fKNO*kno2_err)/(kno3_val));
    nwc_reg->Set("neugen:kno_r123",
			 (kno3_val+fKNO*kno2_err)/(kno3_val));
    nwc_reg->Set("neugen:kno_r133",
			 (kno3_val+fKNO*kno2_err)/(kno3_val));
    nwc_reg->Set("neugen:kno_r143",
			 (kno3_val+fKNO*kno2_err)/(kno3_val));
    }
    nwc_reg->LockValues();
    nwc_reg->LockKeys();

    nwc->SetReweightConfig(nwc_reg);
    mcr->AddWeightCalculator(nwc); 
    fSkzpCalc->SetMCReweight(mcr);
    nwc_reg->Print();
    cout<<"***********************************"<<endl;
  }

  // Setup smoothing
  if(fUseSmoothing) {
    smooth->SetPoints(fFuncBins);
    smooth->SetFunc(fSmoothFunc);
    smooth->SetMin(fFuncMin);
    smooth->SetMax(fFuncMax);
    smooth->SetNormMin(fNormMin);
    smooth->SetNormMax(fNormMax);
    smooth->SetNormalize(fFuncNorm);
    smooth->Initialize();
    fSkzpCalc->SetMCSmooth(smooth);
  }


  
  // Create Chi2 funcs
  // the names should be of the form of Type/Name
  // where Type refers to the concrete class of the
  // object and Name is its name
  std::vector<std::string>::iterator st_iter=fChi2FuncNames.begin();
  std::vector<std::string>::iterator st_end=fChi2FuncNames.end();

  for( ; st_iter!=st_end ; ++st_iter) {

   string type=(*st_iter).substr(0,(*st_iter).find("/"));
   string name=(*st_iter).substr((*st_iter).find("/")+1);

   Chi2FuncPtr chi2 = GenericFactory<Chi2Func>::Instance().Create(type);
   chi2->SetName(name);

   // Set Properties common to all
   PropIter prop_iter=fProperties.find("AllChi2s");
   if(prop_iter!=fProperties.end()) {
     
     std::vector<std::string>::iterator pr_iter=(*prop_iter).second.begin();
     std::vector<std::string>::iterator pr_end=(*prop_iter).second.end();
     for(  ; pr_iter!=pr_end; ++pr_iter) {chi2->SetProperty(*pr_iter);}
   }

   // Set own Properties
   prop_iter=fProperties.find(name);
   if(prop_iter!=fProperties.end()) {
     
     std::vector<std::string>::iterator pr_iter=prop_iter->second.begin();
     std::vector<std::string>::iterator pr_end=prop_iter->second.end();
     for(  ; pr_iter!=pr_end; ++pr_iter) chi2->SetProperty(*pr_iter);
   }

   // Initialize the func
   if(!chi2->Init()) return false;
   fChi2Funcs.push_back(chi2);
  }


  // Create weights and add to the SkzpCalculator
  st_iter=fWeightNames.begin();
  st_end=fWeightNames.end();
  
  for( ; st_iter!=st_end ; ++st_iter) {
    string type=(*st_iter).substr(0,(*st_iter).find("/"));
    string name=(*st_iter).substr((*st_iter).find("/")+1);
    
    if(name=="") name=type;

    SkzpWeightPtr weight = GenericFactory<SkzpWeight>::Instance().Create(type);
    weight->SetName(name);
    
    // Set Properties common to all
    PropIter prop_iter=fProperties.find("AllWeights");
    if(prop_iter!=fProperties.end()) {
      std::vector<std::string>::iterator pr_iter=(*prop_iter).second.begin();
      std::vector<std::string>::iterator pr_end=(*prop_iter).second.end();
      for(  ; pr_iter!=pr_end; ++pr_iter) weight->SetProperty(*pr_iter);
    }
    
   // Set own Properties
   prop_iter=fProperties.find(name);
   if(prop_iter!=fProperties.end()) {
     std::vector<std::string>::iterator pr_iter=prop_iter->second.begin();
     std::vector<std::string>::iterator pr_end=prop_iter->second.end();
     for(  ; pr_iter!=pr_end; ++pr_iter) weight->SetProperty(*pr_iter);
   }

   success=weight->Init();
   if(!success) return false;

   fSkzpCalc->AddSkzpWeight(weight);
  }


  // Setup minuit
  TVirtualFitter::SetDefaultFitter(fFitter.c_str());
  TVirtualFitter::SetMaxIterations(fMaxIterations);
  fitter = TVirtualFitter::Fitter(0, 50);
  arglist[0]=fStrategy;
  fitter->ExecuteCommand("SET STR",arglist,1);
  fitter->PrintResults(3,0);
  fitter->SetMaxIterations(fMaxIterations);
  fitter->SetErrorDef(fERRDEF);
  fitter->SetFCN(fit_fcn);


  // Define Parameters

  // Get the Paramters from the skzpcalculator
  fSkzpCalc->GetParameters(fParameters);

  
  fNPar=fParameters.size();

  Registry *reg=0;
  TFile *file=0;
  double npar_loaded=0;
  if(fFitParFile!="") {

   file=new TFile(fFitParFile.c_str());
    reg=(Registry*)file->Get("Par");

    if(!reg) {
      cout<<"Could not find Registry Par in file "<<fFitParFile
	  <<". Not loading any values"<<endl;
      fFitParFile="";
    }
    reg->Get("NPar",npar_loaded);

    if(npar_loaded==0) {
      cout<<"Could not find number of parameters in Registry Par."
	  <<" No values will load"<<endl;
      fFitParFile="";
    }
  }

  for(int i=0; i<fNPar ; ++i) {

    if(fFitParFile!="") {
      int index=-1;
      for(int j=0; j<npar_loaded ; ++j) {
	const char* get_name=Form("par_%i_name",j);
	const char* got_name;

	reg->Get(get_name,got_name);

	string tmp(got_name);
	if(tmp==fParameters[i].name) { 
	  index=j;
	  break;
	}
      }
      
      if(index!=-1) {
	
	double par_val,par_error;
	reg->Get(Form("par_%i_val",index),par_val);
	reg->Get(Form("par_%i_error",index),par_error);
	fParameters[i].init_val=par_val;
	if(par_error!=0) fParameters[i].err=par_error;

      }
      else {
	cout<<"Could not find parameter "<<fParameters[i].name<<" in registry."
	    <<" Using values in file"<<endl;
      }
    }
    fitter->SetParameter(i,fParameters[i].name.c_str(),
			 fParameters[i].init_val,fParameters[i].err,
			 fParameters[i].min,fParameters[i].max);
    if(fParameters[i].fixed) { 
      fitter->FixParameter(i);
      fNParFixed++;
    }
  }
  fNParFree=fNPar-fNParFixed;
  
  
  fOut=new TFile(fOutputFile.c_str(),"RECREATE");
  fOut->cd();

  
  return true;

}

//......................................................................

void SkzpFitter::Write()
{
  fOut->cd();

  // a list of the free parameter names
  std::vector<std::string> free_names;
  // Get the final parameter values and errors from Minuit
  Registry fPar_val;
  fPar_val.SetName("vals");
  std::vector<double> par,err;
  for(int i=0; i<fNPar ; ++i) {
    fParameters[i].val=fitter->GetParameter(i);
    fParameters[i].err=fitter->GetParError(i);
    par.push_back(fParameters[i].val);
    err.push_back(fParameters[i].err);

    char * name=Form("par_%i_name",i);
    fPar_reg.Set(name,fParameters[i].name.c_str());
    name=Form("par_%i_val",i);
    fPar_reg.Set(name,fParameters[i].val);
    fPar_val.Set(fParameters[i].name.c_str(),fParameters[i].val);
    name=Form("par_%i_error",i);
    fPar_reg.Set(name,fParameters[i].err);
    name=Form("par_%i_fixed",i);
    fPar_reg.Set(name,fParameters[i].fixed);
    name=Form("par_%i_min",i);
    fPar_reg.Set(name,fParameters[i].min);
    name=Form("par_%i_max",i);
    fPar_reg.Set(name,fParameters[i].min);
    name=Form("par_%i_init_val",i);
    fPar_reg.Set(name,fParameters[i].init_val);

    if(!fParameters[i].fixed) free_names.push_back(fParameters[i].name);
  }
  fPar_reg.Set("NPar",fNPar);
  fPar_val.Write();
  fPar_reg.Write();



  // Write tree for backward compatibility

  TTree *infotree = new TTree("info","info");
  int npar=fNPar;
  infotree->Branch("npar",&npar,"npar/I");
  int nfixpar=fNParFixed;

  infotree->Branch("nfixpar",&nfixpar,"nfixpar/I");
  
  
  for(int i=0;i<npar;i++){
    char pn[100];
    sprintf(pn,"par_%d",i);
    char pn2[100];
    sprintf(pn2,"par_%d/D",i);
    infotree->Branch(pn,&(par[i]),pn2);
    
    char en[100];
    sprintf(en,"err_%d",i);
    char en2[100];
    sprintf(en2,"err_%d/D",i);
    infotree->Branch(en,&(err[i]),en2);
  }


  infotree->Fill();
  infotree->Write();


  fSkzpCalc->SetParameters(par);
  // Find and store the best fit values of the chi2 funcs
  std::vector<Chi2FuncPtr>::iterator iter=fChi2Funcs.begin();
  std::vector<Chi2FuncPtr>::iterator end=fChi2Funcs.end();


  
  TDirectoryFile *chi2_dir=new TDirectoryFile("chi2","");
  chi2_dir->cd();
  double totchi2=0;
  int    totndf=0;
  for( ; iter != end ; ++iter) {
    
    double chi2=(*iter)->GetChi2();
    int ndf=(*iter)->GetNDF();
    string name=(*iter)->GetName()+":chi2";
    (*iter)->Write();

    fChi2_reg.Set(name.c_str(),chi2);
    name=(*iter)->GetName()+":ndf";
    fChi2_reg.Set(name.c_str(),ndf);

    totchi2+=chi2;
    totndf+=ndf;
  }

  fOut->cd();
  TDirectoryFile *weights_dir=new TDirectoryFile("weights","");
  weights_dir->cd();
  fSkzpCalc->Write();


  double pt=fSkzpCalc->GetSkzpPT();
  fChi2_reg.Set("weight:pt",pt);
  totchi2+=pt;
  totndf-=fNParFree;

  double param_pt=0.;
  // Add parameter penalty terms
  for(int i=0; i<fNPar ; ++i) {
    if(fParameters[i].pen>1e-4) {
      param_pt+= (fParameters[i].val*fParameters[i].val)/
	(fParameters[i].pen*fParameters[i].pen);
    }
  }

  fChi2_reg.Set("param:pt",param_pt);
  totchi2+=param_pt;

  fOut->cd();
  string name="total:chi2";
  fChi2_reg.Set(name.c_str(),totchi2);
  name="total:ndf";
  fChi2_reg.Set(name.c_str(),totndf);
  fChi2_reg.Write();

  if(fNParFree>0 ) {
    //Get the covariance matrix
   TMatrixT<double> err_matrix(fNParFree,fNParFree);
   TMatrixT<double> corr_matrix(fNParFree,fNParFree);
   for (Int_t i=0;i<fNParFree;i++)
     for (Int_t j=0;j<fNParFree;j++) {
      
      err_matrix(i,j)=fitter->GetCovarianceMatrixElement(i,j);
      corr_matrix(i,j)=fitter->GetCovarianceMatrixElement(i,j)/
	sqrt(fitter->GetCovarianceMatrixElement(i,i)*
	     fitter->GetCovarianceMatrixElement(j,j));
     }

   err_matrix.Write("errorMatrix");
   corr_matrix.Write("CorrMatrix");



  TH2D *corr=new TH2D("Corr","Correlation Matrix",fNParFree,0,fNParFree,
		      fNParFree,0,fNParFree);
  corr->SetLabelFont(62);
  // Set Bin Labels
  for(int i=0;i<fNParFree;i++) {
    corr->GetXaxis()->SetBinLabel(i+1,free_names[i].c_str());
    corr->GetYaxis()->SetBinLabel(i+1,free_names[i].c_str());
  }

  corr->GetXaxis()->LabelsOption("v");

  for(int i=0;i<fNParFree;i++) 
    for(int j=0;j<fNParFree;j++) {
      corr->SetBinContent(i+1,j+1,corr_matrix(i,j));
    } 
  
  corr->Write();

  TDirectoryFile *corr_dir=new TDirectoryFile("corr","");
  corr_dir->cd();
  for(int i=0;i<fNParFree;i++) {
    
    string hname=free_names[i];

    // Replace + and - in names
    if(hname.find("+")!=std::string::npos) {
      hname.replace(hname.find("+"),1,"plus",4);
    }
    if(hname.find("-")!=std::string::npos) {
      hname.replace(hname.find("-"),1,"minus",5);
    }

    const char *name=Form("corr_%s",hname.c_str());
    const char *title=Form("Correlations for %s",free_names[i].c_str());

    TH1D *h=new TH1D(name,title,fNParFree,0,fNParFree);

    for(int j=0;j<fNParFree;j++) {
      h->GetXaxis()->SetBinLabel(j+1,free_names[j].c_str());
      h->SetBinContent(j+1,corr->GetBinContent(i+1,j+1));
    }
    h->GetXaxis()->LabelsOption("v");
    h->Write();
 }
  }
fOut->cd();
  

  
  
}

void SkzpFitter::SetParameters(vector<double> par)
{
  MSG("SkzpFitter",Msg::kVerbose)<<" Parameter Values"<<endl;;
  for(int i=0;i<fNPar ; ++i) {
    if(!fParameters[i].fixed) {
      MSG("SkzpFitter",Msg::kVerbose)<<"    "<<fParameters[i].name<<" :"<<par[i]<<endl;
    }
    fParameters[i].val=par[i];
  }

  fSkzpCalc->SetParameters(par);
}

//......................................................................

double SkzpFitter::GetChi2() 
{
  double chi2=0;
  fNDF=0;
  std::vector<Chi2FuncPtr>::iterator iter=fChi2Funcs.begin();
  std::vector<Chi2FuncPtr>::iterator end=fChi2Funcs.end();

  for( ; iter != end ; ++iter) {
    
    double mychi2=(*iter)->GetChi2();
    int   myNDF=(*iter)->GetNDF();
    if(fIterations==0) { 
      MSG("SkzpFitter",Msg::kInfo) << "Initial chi2 for "<<(*iter)->GetName()
				       <<" "<<mychi2<<" / "<<myNDF<<endl;
    }
    chi2+=mychi2;
    fNDF+=myNDF;
  }
  
  
  chi2+=fSkzpCalc->GetSkzpPT();

  // Add parameter penalty terms
  for(int i=0; i<fNPar ; ++i) {
    if(fParameters[i].pen>1e-6) {
      chi2+= (fParameters[i].val*fParameters[i].val)/
	(fParameters[i].pen*fParameters[i].pen);
    }
  }
  


  return chi2;
}


void  SkzpFitter::Scan(const std::string &cmd)
{

  // 
  for(int i=0; i<fNPar; ++i) {
    fParameters[i].val=fitter->GetParameter(i);
    fParameters[i].err=fitter->GetParError(i);
  }
    

  // Get the properties for this scan
  vector<std::string> params;   // list of parameters by name
  int n_points=fScanPoints;     // number of points
  bool save_values=false;       // set parameter value to lowest chi2
  bool write=true;              // write out chi2 graphs
  double step=-1;               // overide step size


  //  Are there any properties
  string scan_cmd=cmd.substr(0,cmd.find("/"));
  string name=cmd.substr(cmd.find("/")+1);
  
  if(fProperties.find(name)!=fProperties.end()) {
    vector<string>::const_iterator iter=fProperties[name].begin();
    vector<string>::const_iterator end=fProperties[name].end();

    for( ; iter!=end ; ++iter) {
      
      // Split into property and value
      
      string property=(*iter).substr(0,(*iter).find("="));
      string val=(*iter).substr((*iter).find("=")+1);
      
      if(property=="Parameters")  Parsers::parse_val(val,params);
      else if(property=="Save")   Parsers::parse_val(val,save_values);
      else if(property=="Write")  Parsers::parse_val(val,write);
      else if(property=="Points") Parsers::parse_val(val,n_points);
      else if(property=="Step")   Parsers::parse_val(val,step);
    }
  }
  
  vector<int> param_list; // list of parameter by index

  // if it is emtpy scan over all parameters that are not fixed
  if(params.empty()) {
    for(int i=0; i<fNPar ; ++i) {
      if(!fParameters[i].fixed) param_list.push_back(i);
    }
  }
  else {

    vector<string>::const_iterator p_iter=params.begin();
    vector<string>::const_iterator p_end=params.end();
    
    for( ; p_iter!=p_end ; ++p_iter) {
      bool found=false;
      
      for(int i=0; i<fNPar ; ++i) {
	
	if((*p_iter)==fParameters[i].name) {
	  param_list.push_back(i);
	  found=true;
	}
      }
      if(!found) {
	MSG("SkzpFitter",Msg::kError) << "Cannot find parameter " << *p_iter 
				      << " it will not be scanned" << endl;
      }
    }
  }

  if(write) {
    fOut->cd();
    if(!fScanDir) {
      fScanDir=new TDirectoryFile("scan","");
    }
  }
  

  vector<int>::const_iterator iter=param_list.begin();
  vector<int>::const_iterator end=param_list.end();
  for( ; iter!=end ; ++iter) {

    double step_size;
    if(step==-1) step_size=fParameters[*iter].err;
    else step_size=step;

    // find max and min points
    double max=fParameters[*iter].val+n_points/2.*step_size;
    double min=fParameters[*iter].val-n_points/2.*step_size;

    Scan(name,*iter,n_points,min,max,save_values,write);
    
  }

}

void  SkzpFitter::Scan(const std::string &id,const int index,const int &npoints,
		       const double &start,const double &end,const bool &save,
		       const bool &write)
{

  double chi2_min=1e9;
  double val_min=1e9;
  TGraph *gr=new TGraph;;
  MSG("SkzpFitter",Msg::kInfo) << " Scanning "<<fParameters[index].name<<" from "
			       <<start<<" to "<<end<<" in "<<npoints<<" steps"<<endl;
  double default_value=fParameters[index].val;
  vector<double> values;
  for(int i=0;i<fNPar; ++i) {
    values.push_back(fParameters[i].val);
  }

  double step=(end-start)/npoints;

  for(int i=0; i<=npoints; ++i) {

    values[index]=start+i*step;


    SetParameters(values);

    double chi2=GetChi2();
    MSG("SkzpFitter",Msg::kDebug) 
      << "   Current value "<<values[index]<< "   chi2: "<<chi2<<endl;
    if(chi2<chi2_min) {
      chi2_min=chi2;
      val_min=values[index];
    }

    gr->SetPoint(i,values[index],chi2);
    IncrementIteration();
  }

  if(save) {
    fitter->SetParameter(index,fParameters[index].name.c_str(),
			 val_min,fParameters[index].err,
			 fParameters[index].min,fParameters[index].max);
  }else fParameters[index].val=default_value;
  

  if(write) {
    if(!fScanDir) {
      fOut->cd();
      fScanDir=new TDirectoryFile("scan","");
    }
    fScanDir->cd();

    string gr_name=id+"_"+fParameters[index].name;

    // replace the '+' and '-' to word values
    if(gr_name.find("+")!=std::string::npos) {
      gr_name.replace(gr_name.find("+"),1,"plus",4);
    }
    if(gr_name.find("-")!=std::string::npos) {
      gr_name.replace(gr_name.find("-"),1,"minus",5);
    }
    gr->SetName(gr_name.c_str());
    gr->Write();
  }

}


    

    
  
////////////////////////////////////////////////////////////////////////
void fit_fcn(Int_t &, double *, double &f, double *par, int)
{

  TStopwatch timer;

  std::vector<double> vec;
  for( int i=0 ; i<gFitter->GetNPar() ; ++i) vec.push_back(par[i]);

  
  gFitter->SetParameters(vec);
  
  // You can only get ndf after you call getchi2
  f=gFitter->GetChi2();
  MSG("SkzpFitter",Msg::kDebug) << "Iteration " <<gFitter->GetIteration()<<" chi2:"<<f<<endl;
  int ndf=gFitter->GetNDF();
  ndf-=gFitter->GetNParFree();
  
  gFitter->IncrementIteration();
  //  timer.Print("u");
}