#include <iostream>
#include <iomanip>
#include <string>

//root
#include "TLegend.h"
#include "TMinuit.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TText.h"
#include "TGraph.h"
#include "TLine.h"
#include "TMatrixT.h"
#include "TStopwatch.h"
#include "TDatime.h"

//local
#include "FitBeam.h"
#include "FitPTRW.h"
#include "FitData.h"
#include "FitConf.h"

//minos
#include "Registry/Registry.h"
#include "MCReweight/Zbeam.h"
#include "MCReweight/Zfluk.h"
#include "MCReweight/NuParent.h"
#include "MCReweight/MCReweight.h"
#include "MCReweight/MCEventInfo.h"
//#include "MCReweight/NeugenWeightCalculator.h"

//#include "/afs/fnal/files/home/room3/zarko/.rootalias.C"

using namespace std;
int NPAR;

vector<double> pirat_na49;
vector<double> pz_na49;

void DoTheFit(const char* outname,const char* inputFile);

void SaveBestParamFit(FitPTRW *dobj, double *bpar, double *err, 
		      TMatrixT<double> errmatrix, const char *outname);

void fcn(Int_t &/*npar*/, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t iflag)
{

  FitPTRW *dobj = dynamic_cast<FitPTRW *>(gMinuit->GetObjectFit());
  FitConf *fc = FitConf::GetInstance();
  if (dobj->iterations==0){
    TFile* f=new TFile("data/na49.root");
    if (f) {
      TGraph *pina49=dynamic_cast<TGraph*> (f->Get("pina49"));
      for (int i=0;i<pina49->GetN();i++) {
	pz_na49.push_back(pina49->GetX()[i]);
	pirat_na49.push_back(pina49->GetY()[i]);
      }
    }
    f->Close();
    delete f;
    f=0;
    cout<<"inital values of params: "<<iflag<<endl;
    for(int j=0;j<NPAR;j++){
      cout<<j<<" "<<par[j]<<endl;
    }
  }

  int nZbeam=fc->GetZbeamPar().size();
  int nZfluk=fc->GetZflukPar().size();
  int nDet=fc->GetDetPar().size();

  std::vector<double> zbmVec;
  for (int i=0;i<nZbeam;i++) {zbmVec.push_back(par[i]);}
  std::vector<double> zflVec;
  for (int i=0;i<nZfluk;i++) {zflVec.push_back(par[i+nZbeam]);}
  std::vector<double> detVec;
  for (int i=0;i<nDet;i++) {detVec.push_back(par[i+nZbeam+nZfluk]);}

  dobj->SetZflukPar(zflVec);
  dobj->SetZbeamPar(zbmVec);
  dobj->SetDetPar(detVec);

  dobj->FillRWHist();  

  double totchi2=0.0;
  double inittotchi2=0.0;

  double totndf=0.;
  //loop over mc beam configs
  std::map<FitBeam::FitBeam_t, double> chi2;
  double mychi2=0.;
  int myndf=0;
  for(map<FitBeam::FitBeam_t,  TH1D * >::const_iterator bit = dobj->fDataHist.begin(); 
      bit != dobj->fDataHist.end(); ++bit) { 
   
    mychi2=0.;
    myndf=0;
    const FitBeam::FitBeam_t beam_type = bit -> first;
    chi2[beam_type]=0.;
    //compute chi2
    for(int i=1;i<=dobj->fReWeightHist[beam_type]->GetNbinsX();i++){
      double dc = dobj->fDataHist[beam_type]->GetBinContent(i);
      double mcc = dobj->fReWeightHist[beam_type]->GetBinContent(i);
      double edc = dobj->fDataHist[beam_type]->GetBinError(i);
      double emcc = dobj->fReWeightHist[beam_type]->GetBinError(i);
      if(edc*edc+emcc*emcc==0){
	continue;
      }
      mychi2+=(dc-mcc)*(dc-mcc)/(edc*edc+emcc*emcc);      
      if(dc!=0&&mcc!=0){
	myndf++;
      }
    }
    chi2[beam_type]=mychi2;
    totchi2+=mychi2;
    totndf+=myndf;
    
    if(dobj->iterations==0){
      cout<<FitBeam::NeutrinoTypeAsString(beam_type)<<" "
	  <<FitBeam::BeamTypeAsString(beam_type)
	  <<" initial chi2 (non-weighted MC): "<<endl;
      cout<<mychi2<<"/"<<myndf<<endl;
      dobj->initchi2[beam_type]=mychi2;
      dobj->initndf[beam_type]= myndf;
      dobj->bestchi2[beam_type]=mychi2;
      dobj->bestndf[beam_type]=myndf;  
      inittotchi2+=dobj->initchi2[beam_type];
    }
  }
  mychi2=0;
  myndf=0;
  //add penalty terms
  //penalty term for shift in mean pt
  //first for pi+
  double ptshift = dobj->fZfluk.GetPTshift(8);
  double ptPen=fc->GetPiPlusPTpenalty();
  if (ptPen!=0.) mychi2+=ptshift*ptshift/(ptPen*ptPen);
  
  //then for K+
  ptshift = dobj->fZfluk.GetPTshift(11);
  ptPen=fc->GetPiMinusPTpenalty();
  if (ptPen!=0.) mychi2+=ptshift*ptshift/(ptPen*ptPen);

          
  // check fractional changes in number of pi & K and add penalty terms 
  // for pi+/pi-, k+/k-, k+/pi+, k-/pi-
  // do this for n_xf_slices
  int n_xf_slices=4;
  double binwidth_pi=15.; // for pi+/pi- evaluate change in ratio over 0.<xF<1. 
  double binwidth_k=15.;  // for other ratios data reliable only in 0.<xF<0.5 range

  std::vector<double> pt_slices;
  pt_slices.push_back(0.);
  pt_slices.push_back(0.1);
  pt_slices.push_back(0.3);
  pt_slices.push_back(0.5);
  pt_slices.push_back(1.0);
  int n_pt_slices=int(pt_slices.size())-1;
  for (int i=0;i<n_xf_slices;i++) {
    for (int j=0;j<n_pt_slices;j++) {
      //      double npp1=dobj->fZfluk.GetNFrac(8,double(i)*binwidth_pi,double(i+1)*binwidth_pi,
      //                                        pt_slices[j],pt_slices[j+1]);
      //      double npm1=dobj->fZfluk.GetNFrac(9,double(i)*binwidth_pi,double(i+1)*binwidth_pi,
      //               				pt_slices[j],pt_slices[j+1]);
      double npp2=dobj->fZfluk.GetNFrac(8,double(i)*binwidth_k,double(i+1)*binwidth_k,
					pt_slices[j],pt_slices[j+1]);
      double npm2=dobj->fZfluk.GetNFrac(9,double(i)*binwidth_k,double(i+1)*binwidth_k,
					pt_slices[j],pt_slices[j+1]);
      double nkp=dobj->fZfluk.GetNFrac(11,double(i)*binwidth_k,double(i+1)*binwidth_k,
				       pt_slices[j],pt_slices[j+1]);
      double nkm=dobj->fZfluk.GetNFrac(12,double(i)*binwidth_k,double(i+1)*binwidth_k,
				       pt_slices[j],pt_slices[j+1]);
      
      //      mychi2+=(npp1/npm1-1.)*(npp1/npm1-1.)/(0.2*0.2)/(double(n_xf_slices)*double(n_pt_slices)); // pi+ / pi-
      mychi2+=(nkp/nkm-1.)*(nkp/nkm-1.)/(0.2*0.2)/(double(n_xf_slices)*double(n_pt_slices));     // K+  / K-
      mychi2+=(nkp/npp2-1.)*(nkp/npp2-1.)/(0.2*0.2)/(double(n_xf_slices)*double(n_pt_slices));   // K+  / pi+
      mychi2+=(nkm/npm2-1.)*(nkm/npm2-1.)/(0.2*0.2)/(double(n_xf_slices)*double(n_pt_slices));   // K-  / pi-
    }
  }
  
  // pi+/pi- ratio penalty term (NA49)
  TH1D* piratio=dynamic_cast<TH1D*> ((dobj->fZfluk.GetReweightedPTXF(8))->ProjectionX());
  TH1D* piminus=dynamic_cast<TH1D*> ((dobj->fZfluk.GetReweightedPTXF(9))->ProjectionX());
  piratio->Divide(piminus);
  for (unsigned int i=0;i<pz_na49.size();i++) {
    int bin = piratio->FindBin(pz_na49[i]);
    double rwrat=piratio->GetBinContent(bin);
    mychi2+= (rwrat-pirat_na49[i])*(rwrat-pirat_na49[i])/(0.05*0.05);
  }
  
  for(int j=0;j<nZbeam;j++){//these are for other beam syst.
    double penalty=fc->GetZbeamParPenalty()[j];
    if (penalty!=0.) mychi2+=dobj->GetZbeamPar()[j]*dobj->GetZbeamPar()[j]/(penalty*penalty);
  }
  
  for(int j=0;j<nDet;j++){//these are for detector parameters
    double penalty=fc->GetDetParPenalty()[j];
    if (penalty!=0.) mychi2+=dobj->GetDetPar()[j]*dobj->GetDetPar()[j]/(penalty*penalty);
  }
  //these are for numubar xsec
  for (int j=0;j<4;j++) {
    double a=dobj->GetNumubarXsecPar(0);
    double xp=dobj->GetNumubarXsecPar(1);
    double true_enu=(2*j+1)*xp/8.;
    mychi2+=pow((a+2.*(1.-a)*true_enu/xp+(a-1.)*true_enu*true_enu/(xp*xp)-1.),2)/(0.3*0.3)/4.;
  }
  
  totchi2+=mychi2;
  totndf-=fc->GetNpar();  

  if(dobj->iterations==0) {
    cout<<"Total initial chi2/NDF (non-weighted MC): "
	<<inittotchi2<<" / "<<totndf<<endl;
    dobj->inittotchi2=inittotchi2;
    dobj->besttotchi2=totchi2;
  }
  
  if(totchi2<dobj->besttotchi2) {
    dobj->besttotchi2=totchi2;
    for(map<FitBeam::FitBeam_t,  TH1D * >::const_iterator bit = dobj->fDataHist.begin(); 
	bit != dobj->fDataHist.end(); ++bit) {
      dobj->bestchi2[bit->first]=chi2[bit->first];
    }
  }
  dobj->iterations++;
  if( dobj->iterations%1000 == 0 ) cout<<"Finished iteration "<<dobj->iterations<<endl;
 
  f=totchi2;
  return;
}

void DoTheFit(const char* outname,const char* inputFile)
{
  TStopwatch tstp;
  tstp.Start();
  cout<<"****************************************************"<<endl;
  TDatime *now = new TDatime();
  now->Print();
  delete now;
  cout<<"Starting new run"<<endl;
  cout<<"outname "<<outname<<endl;
  FitPTRW *dobj = new FitPTRW(); 
  FitConf *fc=FitConf::GetInstance();
  fc->ReadInput(inputFile);
  NPAR=fc->GetZbeamPar().size()+fc->GetZflukPar().size()+fc->GetDetPar().size();
  
  bool doinuk=fc->DoInuke();
  if(fc->DoInuke()){
    cout<<"will do inuke"<<endl;
  }
  double cutdpid=fc->GetPIDcut();
  cout<<"cutting on pid>"<<cutdpid<<endl;
  bool cutlowy=fc->DoCutLowY();
  if(cutlowy){
    cout<<"cutting lowy "<<endl;
  }

  if(fc->DoSmooth()) {
    dobj->UseSmooth(true);
    dobj->SetSmoothWidth(fc->GetSmoothWidth());
  }
  
  if (fc->DoNeugen()||fc->DoXSecRW()){
    dobj->UseNeugen(true);
    double maqe=fc->GetMAQE();
    double mares=fc->GetMARES();
    double kno=fc->GetKNO();
    dobj->rwtconfig->UnLockValues();
    dobj->rwtconfig->UnLockKeys();
    dobj->rwtconfig->Clear();
    //  dobj->rwtconfig->Set("neugen:config_name",fc->GetNeugenConfName().c_str());
    //     dobj->rwtconfig->Set("neugen:config_no",fc->GetNeugenConfNumber());
    dobj->rwtconfig->Set("neugen:use_scale_factors",1);
    dobj->rwtconfig->Set("neugen:ma_qe",maqe);
    dobj->rwtconfig->Set("neugen:ma_res",mares);    


    // Currently only using 2 kno parameters
    dobj->rwtconfig->Set("neugen:kno_r112",(0.1+kno*0.1)/(0.1));
    dobj->rwtconfig->Set("neugen:kno_r122",(0.3+kno*0.1)/(0.3));

    //     dobj->rwtconfig->Set("neugen:kno_r112",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r122",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r132",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r142",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r113",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r123",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r133",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r143",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r212",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r222",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r232",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r242",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r213",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r223",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r233",kno);
    //     dobj->rwtconfig->Set("neugen:kno_r243",kno);  
    
    dobj->rwtconfig->LockValues();
    dobj->rwtconfig->LockKeys();
  }
  cout<<"Reweight config:"<<endl;
  dobj->rwtconfig->Print();
  cout<<"*********************"<<endl;

  //Check if MRCC histogram should be loaded for nues
  if (fc->GetNueMRCCFile()!="")
    {
      TFile* ff=new TFile(fc->GetNueMRCCFile().c_str());
      assert (ff->IsOpen() && "Could not find MRCC file");
      
      TH1F* MRCChist=dynamic_cast<TH1F*> (ff->Get(fc->GetNueMRCCHistName().c_str()));
      dobj->SetNueMRCCRatio(MRCChist);
      delete MRCChist;
      MRCChist=0;
      ff->Close();
      delete ff;
      ff=0;
      }
    

  //Fill the data and MC histograms:
  for (Int_t i=0;i<fc->GetNbeams();i++)
    {
    
      FitBeam::FitBeam_t beam_type=FitBeam::IntToEnum(fc->GetBeamType()[i]);
      dobj->SetDataPOTS(beam_type,fc->GetDataPOT()[i]);
      dobj->SetMCPOTS(beam_type,fc->GetMCPOT()[i]);
      dobj->FillDataHist(fc->GetDataPath()[i],fc->GetDataName()[i],beam_type);
      dobj->FillMCVectors(fc->GetMCpath()[i],fc->GetMCname()[i],beam_type,doinuk);
    }
  
  TMinuit *min = new TMinuit(NPAR);
  gMinuit = min;
  min->SetFCN(fcn);
  min->SetObjectFit(dobj);

  Double_t arglist[10];
  Int_t ierflg = 0;
 
  Double_t vstart[NPAR],vmin[NPAR],vmax[NPAR],errstart[NPAR];
		
  Int_t nZbmPar=fc->GetZbeamPar().size();
  Int_t nZflPar=fc->GetZflukPar().size();
  Int_t nDetPar=fc->GetDetPar().size();

  for (Int_t i=0;i<nZbmPar;i++)
    {
      vstart[i]=fc->GetZbeamPar()[i];
      vmin[i]=fc->GetZbeamParMin()[i];
      vmax[i]=fc->GetZbeamParMax()[i];
    }
  for (Int_t i=0;i<nZflPar;i++)
    {
      vstart[i+nZbmPar]=fc->GetZflukPar()[i];
      vmin[i+nZbmPar]=fc->GetZflukParMin()[i];
      vmax[i+nZbmPar]=fc->GetZflukParMax()[i];
    }
  for (Int_t i=0;i<nDetPar;i++)
    {
      vstart[i+nZbmPar+nZflPar]=fc->GetDetPar()[i];
      vmin[i+nZbmPar+nZflPar]=fc->GetDetParMin()[i];
      vmax[i+nZbmPar+nZflPar]=fc->GetDetParMax()[i];
    }

  for(Int_t i=0;i<NPAR;i++){    
    char parname[5];
    sprintf(parname,"p%02i",i);
    errstart[i]=.1;
    gMinuit->DefineParameter(i,parname,vstart[i],errstart[i],vmin[i],vmax[i]);
    
    if (i<nZbmPar){
      if(fc->IsZbeamParFixed()[i]){
	gMinuit->FixParameter(i);
      }
    }
    else if (i<nZflPar+nZbmPar){
      if(fc->IsZflukParFixed()[i-nZbmPar]){
	gMinuit->FixParameter(i);
      }
    }
    else if (i<nDetPar+nZflPar+nZbmPar){
      if(fc->IsDetParFixed()[i-nZbmPar-nZflPar]){
	gMinuit->FixParameter(i);
      }
    }
  }    
  
  // Now ready for minimization step
  const double ERRDEF=1.;
  arglist[0] = 10000;       //max calls
  arglist[1] = 1.;         //tolerance

  gMinuit->SetErrorDef(ERRDEF);

  for (Int_t i=0;i<Int_t(fc->GetMinuitCommand().size());i++)
    {
      cout << "Starting "<<fc->GetMinuitCommand()[i]<<endl;
      //will use default values for tolerance 
      //      gMinuit->mnexcm(fc->GetMinuitCommand()[i].c_str(),arglist,1,ierflg);
      
      if ((fc->GetMinuitCommand()[i].find("SIMPLEX")!=std::string::npos)||
	  (fc->GetMinuitCommand()[i].find("MIGRAD")!=std::string::npos))
	{
	  arglist[1] = 1.;  //for simplex: tolerance=0.1*ERRDEF by default
	  if (fc->GetMinuitCommand()[i].find("MIGRAD")!=std::string::npos) 
	    arglist[1] = 10.;   //for migrad: tolerance=0.001*arglist[1]*ERRDEF by default
	  gMinuit->mnexcm(fc->GetMinuitCommand()[i].c_str(),arglist,2,ierflg);
	}
      if ((fc->GetMinuitCommand()[i].find("IMP")!=std::string::npos)||
	  (fc->GetMinuitCommand()[i].find("MINOS")!=std::string::npos)||
	  (fc->GetMinuitCommand()[i].find("HES")!=std::string::npos))
	{
	  gMinuit->mnexcm(fc->GetMinuitCommand()[i].c_str(),arglist,1,ierflg);
	}
      
      cout <<"Finished "<< fc->GetMinuitCommand()[i] <<"; errflag is " << ierflg <<endl;
    }
  
  double bestfit[NPAR];
  double errs[NPAR];
  for(int i=0;i<NPAR;i++){
    min->GetParameter(i,bestfit[i],errs[i]);
  }

  Int_t nVarPar=min->GetNumFreePars();
  Double_t cmatrix[nVarPar][nVarPar];
  min->mnemat(&cmatrix[0][0],nVarPar);

  TMatrixT<double> errmatrix(nVarPar,nVarPar);
  cout<<"Error matrix:"<<endl;
  for (Int_t i=0;i<nVarPar;i++)
    {
      for (Int_t j=0;j<nVarPar;j++)
        {
	  errmatrix(i,j)=cmatrix[i][j];
          cout << cmatrix[i][j]<<" ";
	}
      cout<<endl; 
    } 
  cout<<setprecision(5);
  cout<<setw(10);
  cout<<endl; 
  for(int i=0;i<NPAR;i++){
    cout<<bestfit[i]<<" "<<errs[i]<<endl;
  }
  cout<<endl;

  cout<<"Start chi2 "<<dobj->inittotchi2<<endl;
  cout<<"Best chi2 "<<dobj->besttotchi2<<endl;
   
  SaveBestParamFit(dobj,bestfit,errs,errmatrix,outname);
  tstp.Stop();
  tstp.Print();
}

void SaveBestParamFit(FitPTRW *dobj, double *bpar, double *err,
		      TMatrixT<double> errmatrix, const char *outname)
{
  FitConf *fc = FitConf::GetInstance();
  TFile *fout = new TFile(outname,"RECREATE");
  fout->cd();
  
  std::map<FitBeam::FitBeam_t, TH2F*> ptvpz_piplus;
  std::map<FitBeam::FitBeam_t, TH2F*> ptvpzstart_piplus;
  std::map<FitBeam::FitBeam_t, TH2F*> ptvpz_piminus;
  std::map<FitBeam::FitBeam_t, TH2F*> ptvpzstart_piminus;
 
  int nZbeam=fc->GetZbeamPar().size();
  int nZfluk=fc->GetZflukPar().size();
  int nDet=fc->GetDetPar().size();
 
  std::vector<double> zbmVec;
  for (int i=0;i<nZbeam;i++) {zbmVec.push_back(bpar[i]);}
  std::vector<double> zflVec;
  for (int i=0;i<nZfluk;i++) {zflVec.push_back(bpar[i+nZbeam]);}
  std::vector<double> detVec;
  for (int i=0;i<nDet;i++) {detVec.push_back(bpar[i+nZbeam+nZfluk]);}

  dobj->SetZbeamPar(zbmVec);
  dobj->SetZflukPar(zflVec);
  dobj->SetDetPar(detVec);
  
  dobj->FillRWHist();
  
  vector<Int_t> *beamlist=new vector<Int_t>;

  for(map<FitBeam::FitBeam_t, TH1D* >::const_iterator beam_it=dobj->fDataHist.begin(); beam_it!=dobj->fDataHist.end();++beam_it){
    beamlist->push_back(FitBeam::AsInt(beam_it->first));
   
    char pn[100];
    sprintf(pn,"ptvpz_piplus_%d",FitBeam::AsInt(beam_it->first));
    ptvpz_piplus[beam_it->first]=new TH2F(pn,pn,60,0,60,100,0,1);

    char pns[100];
    sprintf(pns,"ptvpzstart_piplus_%d",FitBeam::AsInt(beam_it->first));
    ptvpzstart_piplus[beam_it->first]=new TH2F(pns,pns,60,0,60,100,0,1);
   
    char pnm[100];
    sprintf(pnm,"ptvpz_piminus_%d",FitBeam::AsInt(beam_it->first));
    ptvpz_piminus[beam_it->first]=new TH2F(pnm,pnm,60,0,60,100,0,1);

    char pnms[100];
    sprintf(pnms,"ptvpzstart_piminus_%d",FitBeam::AsInt(beam_it->first));
    ptvpzstart_piminus[beam_it->first]=new TH2F(pnms,pnms,60,0,60,100,0,1);

    cout<<"fit hist "<<dobj->fReWeightHist[beam_it->first]->Integral()<<" data "<<dobj->fDataHist[beam_it->first]->Integral()<<endl;

     //fill nu weighted ptpz histograms
     for (vector<FitData>::const_iterator dit=dobj->fMCData[beam_it->first].begin();
	  dit!=dobj->fMCData[beam_it->first].end();++dit)
       {
	 const double pt = dit -> TargetPt();
	 const double pz = dit -> TargetPz();
	 const int tptype = dit -> ParentType();
	 double funcb;
	 funcb = dobj->fZfluk.GetWeight(tptype,pt,pz);

	 if (tptype==8||tptype==211)
	   {
	     ptvpz_piplus[beam_it->first]->Fill(pz,pt,funcb);
	     ptvpzstart_piplus[beam_it->first]->Fill(pz,pt);
	   }
	 else if (tptype==9||tptype==-211)
	   {	   
	     ptvpz_piplus[beam_it->first]->Fill(pz,pt,funcb);
	     ptvpzstart_piplus[beam_it->first]->Fill(pz,pt);
	   }
       }
  }

  cout<<"Creating tree"<<endl;
  TTree *infotree = new TTree("info","info");
  infotree->Branch("beamlist",&beamlist);
  infotree->Branch("bestchi2",&dobj->besttotchi2,"bestchi2/D");
  infotree->Branch("initchi2",&dobj->inittotchi2,"initchi2/D");
  int npar=NPAR;
  infotree->Branch("npar",&npar,"npar/I");
  int nfixpar=NPAR-fc->GetNpar();
  infotree->Branch("nfixpar",&nfixpar,"nfixpar/I");
  
  //store chi2
  for(map<FitBeam::FitBeam_t, TH1D *>::const_iterator bit=dobj->fDataHist.begin();
      bit!=dobj->fDataHist.end();++bit)
    {
      char bn[100];
      sprintf(bn,"bestchi2_%d",FitBeam::AsInt(bit->first));
      char bn2[100];
      sprintf(bn2,"bestchi2_%d/D",FitBeam::AsInt(bit->first));
      infotree->Branch(bn,&dobj->bestchi2[bit->first],bn2);
      
      char bd[100];
      sprintf(bd,"bestndf_%d",FitBeam::AsInt(bit->first));
      char bd2[100];
      sprintf(bd2,"bestndf_%d/I",FitBeam::AsInt(bit->first));
      infotree->Branch(bd,&dobj->bestndf[bit->first],bd2);
      
      char in[100];
      sprintf(in,"initchi2_%d",FitBeam::AsInt(bit->first));
      char in2[100];
      sprintf(in2,"initchi2_%d/D",FitBeam::AsInt(bit->first));
      infotree->Branch(in,&dobj->initchi2[bit->first],in2);
      
      char id[100];
      sprintf(id,"initndf_%d",FitBeam::AsInt(bit->first));
      char id2[100];
      sprintf(id2,"initndf_%d/I",FitBeam::AsInt(bit->first));
      infotree->Branch(id,&dobj->initndf[bit->first],id2);
    }

  //store best fit parameters
  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,&(bpar[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);
  }

  //store number of POTs in MC and data hists.
  map<FitBeam::FitBeam_t, double > datapots;
  map<FitBeam::FitBeam_t, double > mcpots;
  for(map<FitBeam::FitBeam_t, TH1D *>::const_iterator bit=dobj->fDataHist.begin();
      bit!=dobj->fDataHist.end();++bit) 
    {
      char pdn[100];
      sprintf(pdn,"datapots_%d",FitBeam::AsInt(bit->first));
      char pdn2[100];
      sprintf(pdn2,"datapots_%d/D",FitBeam::AsInt(bit->first));
      datapots[bit->first]=dobj->GetDataPOTS(bit->first);
      infotree->Branch(pdn,&datapots[bit->first],pdn2);    
       
      char pmn[100];
      sprintf(pmn,"mcpots_%d",FitBeam::AsInt(bit->first));
      char pmn2[100];
      sprintf(pmn2,"mcpots_%d/D",FitBeam::AsInt(bit->first));
      mcpots[bit->first]=dobj->GetMCPOTS(bit->first);
      infotree->Branch(pmn,&mcpots[bit->first],pmn2); 
    }
  
  infotree->Fill();

  infotree->Write();
  errmatrix.Write("errorMatrix");
  //save the histograms
  cout<<"Writing histograms"<<endl;
  for(map<FitBeam::FitBeam_t, TH1D *>::const_iterator bit=dobj->fDataHist.begin();
      bit!=dobj->fDataHist.end();++bit) 
    {
      dobj->fDataHist[bit->first]->Write();
      dobj->fMCHist[bit->first]->Write();
      dobj->fReWeightHist[bit->first]->Write();
      ptvpz_piplus[bit->first]->Write();
      ptvpzstart_piplus[bit->first]->Write();
      ptvpz_piminus[bit->first]->Write();
      ptvpzstart_piminus[bit->first]->Write();
    }
  vector<int> parentType;
  vector<string> parentName;
  parentType.push_back(8);   parentName.push_back("PiPlus");
  parentType.push_back(9);   parentName.push_back("PiMinus");
  parentType.push_back(11);  parentName.push_back("KPlus");
  parentType.push_back(12);  parentName.push_back("KMinus");
  parentType.push_back(10);  parentName.push_back("K0L");
  const int nPTXF=int(parentType.size());
  for (int i=0;i<nPTXF;i++) 
    { 
      dobj->fZfluk.GetPTXF(parentType[i])->Write(); 
      dobj->fZfluk.GetReweightedPTXF(parentType[i])->Write(); 
      dobj->fZfluk.GetWeightHistogram(parentType[i])->Write();
    } 

  fout->Close();
  
}

void DrawBestFit(const char* fname, const char* save, bool dosave,const char* filetype="eps")
{
  vector<int> parentType;
  vector<string> parentName;
  parentType.push_back(8);   parentName.push_back("PiPlus");
  parentType.push_back(9);   parentName.push_back("PiMinus");
  parentType.push_back(11);  parentName.push_back("KPlus");
  parentType.push_back(12);  parentName.push_back("KMinus");
  parentType.push_back(10);  parentName.push_back("K0L");

  const int nPTXF=int(parentType.size());
  TH2F* hF05ptxf[nPTXF];
  TH2F* hFitptxf[nPTXF];
  TH2F* hWeightptxf[nPTXF];

  TFile *f = new TFile(fname);

  for (int i=0;i<nPTXF;i++) 
    { 
      hF05ptxf[i]=dynamic_cast<TH2F*>   (f->Get(Form("hF05ptxf%s",parentName[i].c_str()))); 
      hFitptxf[i]=dynamic_cast<TH2F*>   (f->Get(Form("hWeightedPTXF%s",parentName[i].c_str()))); 
      hWeightptxf[i]=dynamic_cast<TH2F*>(f->Get(Form("hWeight%s",parentName[i].c_str()))); 
      
      hF05ptxf[i]->SetDirectory(0);
      hFitptxf[i]->SetDirectory(0);
      hWeightptxf[i]->SetDirectory(0);
    } 

  TTree *infotree = (TTree *)f->Get("info");
  int npar,nfixpar;
  vector<Int_t> *beamlist = NULL;
  
  infotree->SetBranchAddress("npar",&npar);
  infotree->SetBranchAddress("nfixpar",&nfixpar);
  infotree->SetBranchAddress("beamlist",&beamlist);
  
  infotree->GetEntry(0);
  const int NV=npar;
  const int NBMS=(*beamlist).size();
  
  cout<< "Making plots for "<<NV-nfixpar<<" parameter fit to "<<NBMS<<" beams"<<endl;
  double besttotchi2, inittotchi2;
  double bestchi2[NBMS];
  double initchi2[NBMS];
  int bestndf[NBMS];
  int initndf[NBMS];
  double par[NV];
  double err[NV];
  double datapots[NBMS];
  double mcpots[NBMS];

  infotree->SetBranchAddress("bestchi2",&besttotchi2);
  infotree->SetBranchAddress("initchi2",&inittotchi2);
  for(int i=0;i<NBMS;i++){
    char bn[100];
    Int_t blc=(*beamlist)[i];
    sprintf(bn,"bestchi2_%d",blc);
    infotree->SetBranchAddress(bn,&bestchi2[i]);

    char bd[100];
    sprintf(bd,"bestndf_%d",blc);
    infotree->SetBranchAddress(bd,&bestndf[i]);

    char in[100];
    sprintf(in,"initchi2_%d",blc);
    infotree->SetBranchAddress(in,&initchi2[i]);

    char id[100];
    sprintf(id,"initndf_%d",blc);
    infotree->SetBranchAddress(id,&initndf[i]);

    char pdn[100];
    sprintf(pdn,"datapots_%d",blc);
    infotree->SetBranchAddress(pdn,&datapots[i]);

    char pmn[100];
    sprintf(pmn,"mcpots_%d",blc);
    infotree->SetBranchAddress(pmn,&mcpots[i]);
  }
  for(int i=0;i<NV;i++){
    char pn[100];
    sprintf(pn,"par_%d",i);
    infotree->SetBranchAddress(pn,&par[i]);
    char en[100];
    sprintf(en,"err_%d",i);
    infotree->SetBranchAddress(en,&err[i]);
  }
  infotree->GetEntry(0);

  int bndf=0;
  for(int i=0;i<NBMS;i++){
    bndf+=bestndf[i];
  }

  for(int i=0;i<NV;i++){
    cout<<"par "<<i<<" "<<par[i]<<" err "<<err[i]<<endl;
  }

  TH1D *start[NBMS];
  TH1D *rwhist[NBMS];
  TH1D *dmcstart[NBMS];
  TH1D *dmcend[NBMS];
  TH2F *ptvpz_piplus[NBMS];
  TH2F *ptvpzstart_piplus[NBMS];
  TH2F *ptvpz_piminus[NBMS];
  TH2F *ptvpzstart_piminus[NBMS];
  TH1F *datahists[NBMS];

  cout<<"Getting histograms"<<endl;
  for(int i=0;i<NBMS;i++){
   
    FitBeam::FitBeam_t beam_type=FitBeam::IntToEnum((*beamlist)[i]);
    datahists[i]=(TH1F *)f->Get(Form("reco_enu_data_%s_%s", FitBeam::NeutrinoTypeAsString(beam_type).c_str(),
				     FitBeam::BeamTypeAsString(beam_type).c_str()));

    start[i]=(TH1D *)f->Get(Form("start_%s_%s", FitBeam::NeutrinoTypeAsString(beam_type).c_str(),
				     FitBeam::BeamTypeAsString(beam_type).c_str()));
    rwhist[i]=(TH1D *)f->Get(Form("reweight_%s_%s", FitBeam::NeutrinoTypeAsString(beam_type).c_str(),
				     FitBeam::BeamTypeAsString(beam_type).c_str()));

    //lets renormalize histograms to have per 0.5GeV
    double norm=0.5;
    for (int xbin=0;xbin<datahists[i]->GetNbinsX()+1;xbin++)
      {
	datahists[i]->SetBinContent(xbin,datahists[i]->GetBinContent(xbin)*norm/datahists[i]->GetBinWidth(xbin));
	start[i]->SetBinContent(xbin,start[i]->GetBinContent(xbin)*norm/start[i]->GetBinWidth(xbin));
	rwhist[i]->SetBinContent(xbin,rwhist[i]->GetBinContent(xbin)*norm/rwhist[i]->GetBinWidth(xbin));
      }
    ptvpz_piplus[i]=(TH2F *)f->Get(Form("ptvpz_piplus_%i",FitBeam::AsInt(beam_type)));
    ptvpzstart_piplus[i]=(TH2F *)f->Get(Form("ptvpzstart_piplus_%i",FitBeam::AsInt(beam_type)));

    ptvpz_piminus[i]=(TH2F *)f->Get(Form("ptvpz_piminus_%i",FitBeam::AsInt(beam_type)));
    ptvpzstart_piminus[i]=(TH2F *)f->Get(Form("ptvpzstart_piminus_%i",FitBeam::AsInt(beam_type)));

    char dmcs[100];
    sprintf(dmcs,"datamcratstart_%d",FitBeam::AsInt(beam_type));
    dmcstart[i]=dynamic_cast<TH1D*> (datahists[i]->Clone(dmcs));
    char dmce[100];
    sprintf(dmce,"datamcratend_%d",FitBeam::AsInt(beam_type));
    dmcend[i]=dynamic_cast<TH1D*> (datahists[i]->Clone(dmce));

    dmcstart[i]->Reset();
    dmcend[i]->Reset();

    rwhist[i]->SetLineColor(2);
    rwhist[i]->SetMarkerColor(2);

    datahists[i]->SetLineWidth(1);
    rwhist[i]->SetLineWidth(1);
    start[i]->SetLineWidth(1);

    datahists[i]->SetMarkerSize(.8);
    rwhist[i]->SetMarkerSize(.6);
    start[i]->SetMarkerSize(.8);

    datahists[i]->SetMarkerStyle(21);
    rwhist[i]->SetMarkerStyle(20);
    start[i]->SetMarkerStyle(27);
    start[i]->SetLineColor(4);

    datahists[i]->SetTitle(Form("%s_%s;Reco Energy (GeV);CC Events/10^{18}POT/kt/%2.1fGeV",
				FitBeam::NeutrinoTypeAsString(FitBeam::IntToEnum((*beamlist)[i])).c_str(),
				FitBeam::BeamTypeAsString(FitBeam::IntToEnum((*beamlist)[i])).c_str(),
				norm));
    dmcstart[i]->SetTitle(Form(";Reco Energy (GeV);Data/MC",
			       FitBeam::NeutrinoTypeAsString(FitBeam::IntToEnum((*beamlist)[i])).c_str(),
			       FitBeam::BeamTypeAsString(FitBeam::IntToEnum((*beamlist)[i])).c_str()));
   
    ptvpz_piplus[i]->SetLineColor(2);
    ptvpz_piplus[i]->SetMarkerColor(2);
    ptvpzstart_piplus[i]->SetLineColor(4);
    ptvpzstart_piplus[i]->SetMarkerColor(4);
   
    ptvpz_piminus[i]->SetLineColor(2);
    ptvpz_piminus[i]->SetMarkerColor(2);
    ptvpzstart_piminus[i]->SetLineColor(4);
    ptvpzstart_piminus[i]->SetMarkerColor(4);
    
    cout<<"Scaling "<<FitBeam::NeutrinoTypeAsString(beam_type)<<" "
	<< FitBeam::BeamTypeAsString(beam_type) 
	<< " to 1e18 pot * 1/"<<datapots[i]<<endl;

    datahists[i]->Scale(1.e1/datapots[i]);
    start[i]->Scale(1.e1/datapots[i]);
    rwhist[i]->Scale(1.e1/datapots[i]);
    
    datahists[i]->SetStats(0);
    start[i]->SetStats(0);
    rwhist[i]->SetStats(0);

    dmcstart[i]->Divide(datahists[i],start[i]);
    dmcend[i]->Divide(datahists[i],rwhist[i]);

    dmcstart[i]->SetLineColor(4);
    dmcend[i]->SetLineColor(2);

    dmcstart[i]->SetLineWidth(1);
    dmcend[i]->SetLineWidth(1);

    dmcstart[i]->SetStats(0);
    dmcend[i]->SetStats(0);

  
  }
  cout <<endl;
  double tc=0.;
  int tndf=0;
  cout.width(18);
  cout << "Beam type";
  cout.width(16);
  cout<<"Initial chi2";   
  cout.width(8);
  cout<<"NDF";
  cout.width(15);
  cout<<"Final chi2";
  cout.width(8);
  cout<<"NDF"<<endl;
  for(int nt=0;nt<NBMS;nt++){
    //compute chi2
    double mychi2=0.;
    int myndf=0;
    for(int i=1;i<=start[nt]->GetNbinsX();i++){
      double dc = datahists[nt]->GetBinContent(i);
      double mcc = start[nt]->GetBinContent(i);
      double edc = datahists[nt]->GetBinError(i);
      double emcc = start[nt]->GetBinError(i);
      if(edc*edc+emcc*emcc==0){
	continue;
      }
      mychi2+=pow((dc-mcc),2)/(edc*edc+emcc*emcc);      
      if(dc!=0&&mcc!=0){
	myndf++;
      }
    }
   
    cout.width(8);
    cout<< FitBeam::NeutrinoTypeAsString(FitBeam::IntToEnum((*beamlist)[nt]))<<" ";
    cout.width(10);
    cout<< FitBeam::BeamTypeAsString(FitBeam::IntToEnum((*beamlist)[nt]));
    cout.width(13);
    cout<<initchi2[nt];
    cout.width(10);
    cout<<myndf;
    cout.width(13);
    cout<<bestchi2[nt]<<endl;
    tc+=mychi2;
    tndf+=myndf;
  }
  cout.width(19);
  cout<<"Total";
  cout.width(13);
  cout<<tc;
  cout.width(10);
  cout<<tndf;
  cout.width(13);
  cout<<besttotchi2;
  cout.width(10);
  cout<<tndf-NV+nfixpar<<endl;

  TLegend *leg2 = new TLegend(0.55,0.7,0.93,0.88);
  leg2->AddEntry(dmcstart[0],"Fluka 2005","l");
  leg2->AddEntry(dmcend[0],"Tuned Had. Prod.","l");
  leg2->SetFillStyle(0);
  leg2->SetBorderSize(0);
  TLine *l = new TLine(0.5,1,120,1);
  l->SetLineWidth(1);
  for (int i=0;i<NBMS;++i){
    TLegend *leg1 = new TLegend(0.5,0.62,0.88,0.9);
    leg1->AddEntry(datahists[0],"Data","p");
    leg1->AddEntry(start[0],"Fluka 2005","l");
    leg1->AddEntry(dmcstart[0],Form("#chi^{2}=%i",int(initchi2[i])),"");
    leg1->AddEntry(rwhist[0],"Tuned Had. Prod.","l");
    leg1->AddEntry(dmcend[0],Form("#chi^{2}=%i",int(bestchi2[i])),"");
    leg1->SetFillStyle(0);
    leg1->SetBorderSize(0);
    
    TCanvas *can = new TCanvas(Form("can_%i",i),Form("can_%i",i),1000,900);
    can->Divide(1,2);
    can->cd(1);
    can->GetPad(1)->SetLogy();
    datahists[i]->DrawCopy("p");
    datahists[i]->GetXaxis()->CenterTitle();
    datahists[i]->GetYaxis()->CenterTitle();
    start[i]->Draw("sames hist");
    rwhist[i]->Draw("sames hist");
    leg1->Draw();
        
    can->cd(2);
    dmcstart[i]->Draw("hist");
    dmcstart[i]->GetXaxis()->CenterTitle();
    dmcstart[i]->GetYaxis()->CenterTitle();
    dmcstart[i]->GetYaxis()->SetRangeUser(0.6,1.7);
    dmcend[i]->Draw("same hist");
    l->Draw();
    leg2->Draw();
    
    TCanvas *can2 = new TCanvas(Form("can2_%i",i),Form("can2_%i",i),600,600);
    can2->SetLogy(0);
    datahists[i]->GetXaxis()->SetRangeUser(0.,39.);
    datahists[i]->Draw("p");
    start[i]->Draw("sames hist");
    rwhist[i]->Draw("sames hist");
    leg1->Draw();

    if (dosave)
      {
	can->Print(Form("enu_%s_%s.%s",datahists[i]->GetTitle(),save,filetype));
	can2->Print(Form("enu_%s_%s_zoom.%s",datahists[i]->GetTitle(),save,filetype));
      }
  }
  
  TH1D* hF05ptxf_px[nPTXF];
  TH1D* hF05ptxf_py[nPTXF];
  TH1D* hFitptxf_px[nPTXF]; 
  TH1D* hFitptxf_py[nPTXF];

  for (int i=0;i<nPTXF;i++)
    {
      hF05ptxf_px[i]=dynamic_cast<TH1D*> (hF05ptxf[i]->ProjectionX()->Clone());
      hFitptxf_px[i]=dynamic_cast<TH1D*> (hFitptxf[i]->ProjectionX()->Clone());
      hF05ptxf_py[i]=dynamic_cast<TH1D*> (hF05ptxf[i]->ProjectionY()->Clone());
      hFitptxf_py[i]=dynamic_cast<TH1D*> (hFitptxf[i]->ProjectionY()->Clone());

      hF05ptxf_px[i]->SetTitle(";p_{z} (GeV);##/10^{7}POT/bin");
      hF05ptxf_py[i]->SetTitle(";p_{T} (GeV);##/10^{7}POT/bin");
      hFitptxf_px[i]->SetTitle(";p_{z} (GeV);##/10^{7}POT/bin"); 
      hFitptxf_py[i]->SetTitle(";p_{T} (GeV);##/10^{7}POT/bin");
      
      hF05ptxf_px[i]->SetLineColor(4);
      hF05ptxf_py[i]->SetLineColor(4);
      hFitptxf_px[i]->SetLineColor(2);
      hFitptxf_py[i]->SetLineColor(2);
 
      hF05ptxf_px[i]->SetStats(0); 
      hF05ptxf_py[i]->SetStats(0);  
      hFitptxf_px[i]->SetStats(0);  
      hFitptxf_py[i]->SetStats(0); 

      hF05ptxf[i]->SetStats(0);  
      hF05ptxf[i]->SetStats(0);   
      hFitptxf[i]->SetStats(0);   
      hFitptxf[i]->SetStats(0);  

      if (hFitptxf_py[i]->GetMaximum()>hF05ptxf_py[i]->GetMaximum())
	{
	  hF05ptxf_py[i]->SetMaximum(hFitptxf_py[i]->GetMaximum()*1.1);
	}

      TCanvas *canf=new TCanvas(Form("%s_ptpz",parentName[i].c_str()),
				Form("%s_ptpz",parentName[i].c_str()));
      
      canf->Divide(2,2);
      canf->cd(1);
      gPad->SetLogz();
      hF05ptxf[i]->Draw("col");
      canf->cd(2);
      gPad->SetLogz();
      hFitptxf[i]->Draw("col");
      canf->cd(3);
      gPad->SetLogy();
      hF05ptxf_px[i]->Draw();
      hFitptxf_px[i]->Draw("sames");
      canf->cd(4);
      hF05ptxf_py[i]->Draw(); 
      hFitptxf_py[i]->Draw("sames"); 
      
      TLegend *leg3 = new TLegend(0.55,0.6,0.93,0.85); 
      leg3->AddEntry(dmcstart[0],"Fluka 2005","l"); 
      leg3->AddEntry(dmcstart[0],Form("<p_{T}>=%i", 
				      int(hF05ptxf[i]->ProjectionY()->GetMean()*1000)),""); 
      leg3->AddEntry(dmcend[0],"Tuned Had. Prod.","l"); 
      leg3->AddEntry(dmcend[0],Form("<p_{T}>=%i", 
				    int(hFitptxf[i]->ProjectionY()->GetMean()*1000)),""); 
      leg3->SetFillStyle(0); 
      leg3->SetBorderSize(0); 
      leg3->Draw(); 

      if (dosave) 
	canf->Print(Form("%s_ptpz_dist_%s.%s",parentName[i].c_str(),save,filetype)); 

      TCanvas *cw = new TCanvas(Form("%s_weights",parentName[i].c_str()),
				Form("%s_weights",parentName[i].c_str()),1000,800);
      cw->cd();
      gPad->SetLeftMargin(0.12);
      gPad->SetRightMargin(0.12);
      hWeightptxf[i]->GetXaxis()->SetNoExponent();
      hWeightptxf[i]->SetStats(0);
      hWeightptxf[i]->Draw("colz");

      if (dosave) 
	cw->Print(Form("%s_weights_%s.%s",parentName[i].c_str(),save,filetype));
    
      cout <<parentName[i].c_str()<<":"<<endl; 
      cout<<"  start MEAN PT: "<<hF05ptxf[i]->ProjectionY()->GetMean()<<endl; 
      cout<<"  fit MEAN PT: "<<hFitptxf[i]->ProjectionY()->GetMean()<<endl; 

    }

}


void MakeTables(const char* fname, const char* cnfname="")
{
  TFile *f = new TFile(fname);
  vector<int> parentType;
  vector<string> parentName;
  parentType.push_back(8);   parentName.push_back("PiPlus");
  parentType.push_back(9);   parentName.push_back("PiMinus");
  parentType.push_back(11);  parentName.push_back("KPlus");
  parentType.push_back(12);  parentName.push_back("KMinus");
  parentType.push_back(10);  parentName.push_back("K0L");

  std::map<std::string, TH2F* > hF05ptxf;
  std::map<std::string, TH2F* > hFitptxf;
  for (std::vector<string>::iterator itParent=parentName.begin();
       itParent!=parentName.end();itParent++)
    {
      hF05ptxf[*(itParent)]=dynamic_cast<TH2F*> (f->Get(Form("hF05ptxf%s", (*(itParent)).c_str()))); 
      hFitptxf[*(itParent)]=dynamic_cast<TH2F*> (f->Get(Form("hWeightedPTXF%s", (*(itParent)).c_str()))); 
      
      hF05ptxf[*(itParent)]->SetDirectory(0);
      hFitptxf[*(itParent)]->SetDirectory(0);
    } 

  TTree *infotree = (TTree *)f->Get("info");
  int npar;
  vector<Int_t> *beamlist = NULL;
  
  infotree->SetBranchAddress("npar",&npar);
  infotree->SetBranchAddress("beamlist",&beamlist);
  infotree->GetEntry(0);
  const int NV=npar;
  const int NBMS=(*beamlist).size();
  
  double besttotchi2, inittotchi2;
  double bestchi2[NBMS];
  double initchi2[NBMS];
  int bestndf[NBMS];
  int initndf[NBMS];
  double par[NV];
  double err[NV];
  double datapots[NBMS];
  double mcpots[NBMS];

  infotree->SetBranchAddress("bestchi2",&besttotchi2);
  infotree->SetBranchAddress("initchi2",&inittotchi2);
  
  for(int i=0;i<NBMS;i++){
    char bn[100];
    Int_t blc=(*beamlist)[i];
    sprintf(bn,"bestchi2_%d",blc);
    infotree->SetBranchAddress(bn,&bestchi2[i]);

    char bd[100];
    sprintf(bd,"bestndf_%d",blc);
    infotree->SetBranchAddress(bd,&bestndf[i]);

    char in[100];
    sprintf(in,"initchi2_%d",blc);
    infotree->SetBranchAddress(in,&initchi2[i]);

    char id[100];
    sprintf(id,"initndf_%d",blc);
    infotree->SetBranchAddress(id,&initndf[i]);

    char pdn[100];
    sprintf(pdn,"datapots_%d",blc);
    infotree->SetBranchAddress(pdn,&datapots[i]);

    char pmn[100];
    sprintf(pmn,"mcpots_%d",blc);
    infotree->SetBranchAddress(pmn,&mcpots[i]);
  }
  for(int i=0;i<NV;i++){
    char pn[100];
    sprintf(pn,"par_%d",i);
    infotree->SetBranchAddress(pn,&par[i]);
    char en[100];
    sprintf(en,"err_%d",i);
    infotree->SetBranchAddress(en,&err[i]);
  }
  infotree->GetEntry(0);

  // find total NDF
  int bndf=0;
  for(int i=0;i<NBMS;i++){
    bndf+=bestndf[i];
  }
  //subtract number of fitted parameters
  bndf-=NV;
  //check if some of the parameters were fixed and add that back to NDF
  for(int i=0;i<NV;i++){
    if (err[i]==0.)
      bndf++;
  }
  if (cnfname=="") cnfname=fname;
  cout<<"***************TOTAL CHI2**********************"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>Cnf";
  for(int i=0;i<NBMS;i++)
    {
      FitBeam::FitBeam_t beam_type=FitBeam::IntToEnum((*beamlist)[i]);
      cout<<"<th>"<<FitBeam::NeutrinoTypeAsString(beam_type)<<" "
	  <<FitBeam::BeamTypeAsString(beam_type);
    } 
  cout<<"<th>Total<th>NDF<th>Reduced</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname;
  for(int i=0;i<NBMS;i++)
    {
      cout<<"<td>"<<bestchi2[i];
    }
  cout <<"<td>"<<besttotchi2<<"<td>"<<bndf<<"<td>"<<besttotchi2/bndf<<"</tr>"<<endl;
  cout<<"</table>"<<endl;
  
  cout<<"*****************nonbeam params*****************"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>Cnf<th>miscal (down)<th>err<th>shwoffset<th>err<th>nc_numu<th>err<th>nc_numubar<th>err<th>xsec1<th>err<th>xsec2<th>err</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname<<"<td>";
  // for(int j=21;j<25;j++){
    cout<<par[29]<<"<td>"<<err[29]<<"<td>";
    cout<<par[31]<<"<td>"<<err[31]<<"<td>";
    cout<<par[32]<<"<td>"<<err[32]<<"<td>";
    cout<<par[33]<<"<td>"<<err[33]<<"<td>";
    cout<<par[35]<<"<td>"<<err[35]<<"<td>";
    cout<<par[36]<<"<td>"<<err[36]<<"<td>";
    //}
  cout<<"</tr>"<<endl;
  cout<<"</table>"<<endl;

  cout<<"*****************other params*****************"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>Cnf<th>horn offset<th>err<th>baffle<th>err<th>pot<th>err<th>horn I miscal<th>err<th>horn I dist<th>err<th>tgtz (le10z170)<th>err<th>tgtz (le10z185)<th>err<th>tgtz (le10z200)<th>err<th>tgtz (le100)<th>err<th>tgtz (le150)<th>err<th>tgtz (le250)<th>err<th>tgtz (le10 post)<th>err<th>tgtz (le250 post)<th>err</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname<<"<td>";
  for(int j=0;j<13;j++){
    cout<<par[j]<<"<td>"<<err[j]<<"<td>";
  }
  cout<<"</tr>"<<endl;
  cout<<"</table>"<<endl;

  cout<<"*****************hadron prod. params*****************"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>Pi+<th>par0<th>dpar0<th>par1<th>dpar1<th>par2<th>dpar2<th>par3<th>dpar3<th>par4<th>dpar4<th>par5<th>dpar5</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname<<"<td>";
  for(int j=13;j<19;j++){
    cout<<par[j]<<"<td>"<<err[j]<<"<td>";
  }
  cout<<"</tr>"<<endl;
  cout<<"</table>"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>K+<th>par0<th>dpar0<th>par1<th>dpar1<th>par2<th>dpar2<th>par3<th>dpar3<th>par4<th>dpar4<th>par5<th>dpar5</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname<<"<td>";
  for(int j=19;j<25;j++){
    cout<<par[j]<<"<td>"<<err[j]<<"<td>";
  }
  cout<<"</tr>"<<endl;
  cout<<"</table>"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>Pi-<th>par0<th>dpar0<th>par1<th>dpar1</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname<<"<td>";
  for(int j=25;j<27;j++){
    cout<<par[j]<<"<td>"<<err[j]<<"<td>";
  }
  cout<<"</tr>"<<endl;
  cout<<"</table>"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>K-<th>par0<th>dpar0<th>par1<th>dpar1</tr>"<<endl;
  cout<<"<tr><th>"<<cnfname<<"<td>";
  for(int j=27;j<29;j++){
    cout<<par[j]<<"<td>"<<err[j]<<"<td>";
  }
  cout<<"</tr>"<<endl;
  cout<<"</table>"<<endl;

  cout<<"*****************mean pt*****************"<<endl;
  cout<<"<table align=center, border=1, cellspace=10, cellpadding=3>"<<endl;
  cout<<"<tr><th>Cnf<th>Fluka05<th>fit mean pt</t>"<<endl;
  
  for (std::vector<string>::iterator itParent=parentName.begin();
       itParent!=parentName.end();itParent++)
    {
      cout<<Form("<tr><th>%s %s<td>%d<td>%d</tr>",
		 cnfname,
		 (*(itParent)).c_str(),
		 int(hF05ptxf[*(itParent)]->ProjectionY()->GetMean()*1000),
		 int(hFitptxf[*(itParent)]->ProjectionY()->GetMean()*1000))
	  <<endl;
    }
  cout<<"</table>"<<endl;
}