/*------------------------------------------------------------------------------------
*  zetsubou ver1.3.6 (detuned version)
*  author ASANO Hidemitsu
*  Faculty of Science ,Kyoto university
*  Last Modified 2008/2/6 17:43
*  This library is for analysis of expreriment in P3 2007 graduate research 
-------------------------------------------------------------------------------------*/
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>

#include <TROOT.h>
#include <TStyle.h>
#include <TObject.h>
#include <TApplication.h>
#include <TBranch.h>
#include <TSpectrum.h>
#include <TCanvas.h>
#include <TString.h>
#include <TTree.h>
#include <TPad.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TMath.h>

#define HIST 1
#define SCAT 0

class zetsubou : public TObject
{
  public:
  zetsubou();//Constructor
  ~zetsubou();//Destructor

  Bool_t MakeTree(const Char_t *filename);//for getting data and making tree
  Bool_t MakeCanvas(const Char_t *name,Int_t i);//making a canvas corresponding histname or scatname
  void SetHistStyle();//set hist style
  void SetScatStyle();//set scat style 

  //drawing hist and scat
  TH1D *DrawHist(const Char_t *filename,const Char_t *title,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax,Int_t cutmin,Int_t cutmax);
  TH1D *DrawHist(const Char_t *filename,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax,Int_t cutmin,Int_t cutmax);
  TH1D *DrawHist(const Char_t *filename,const Char_t *title,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax);
  TH1D *DrawHist(const Char_t *filename,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax);
  TH1D *HistCut(Double_t X1,Double_t Y1,Double_t X2,Double_t Y2,Double_t X3,Double_t Y3,Double_t X4,Double_t Y4);
  
  TH2D *DrawScat(const Char_t *filename,const Char_t *title,const Char_t *scatname);
  TH2D *DrawScat(const Char_t *filename,const Char_t *scatname);
  TH2D *DrawScat(const Char_t *filename,const Char_t *title,const Char_t *scatname,Int_t xmin, Int_t xmax,Int_t ymin,Int_t ymax);
  TH2D *DrawScat(const Char_t *filename,const Char_t *scatname,Int_t xmin, Int_t xmax,Int_t ymin,Int_t ymax);  
  //fitting hist
  TF1 *gfit(Double_t xmin,Double_t xmax);
  TF1 *gfit(Double_t xmin,Double_t xmax,Double_t height,Double_t mean,Double_t sigma);
  TF1 *gfit(Double_t xmin,Double_t xmax,Double_t height,Double_t mean,Double_t sigma,Double_t katamuki, Double_t seppen);
  TF1 *gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2);
  TF1 *gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2,Double_t katamuki ,Double_t seppen);
  TF1 *gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2,Double_t height3,Double_t mean3,Double_t sigma3);
  TF1 *gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2,Double_t height3,Double_t mean3,Double_t sigma3,Double_t katamuki ,Double_t seppen);
  TF1 *lgfit(Double_t xmin,Double_t xmax,Double_t width,Double_t mp,Double_t area,Double_t gsigma);
  
   //Searching Peaks Automatically by using TSpectrum class
  TSpectrum *Tspfit(Double_t xmin,Double_t xmax);
  TSpectrum *Tspfit(Double_t xmin,Double_t xmax,Double_t threshold);
  TSpectrum *Tspfit(Double_t xmin,Double_t xmax,Double_t threshold,Double_t resolution);
  
  inline TTree *GetTree() { return tr;}; //return pointer to current tree
 
  private:
  
  std::string currentfname;//input current filename 
  TTree *tr;
  TH1D  *hist;
  TH2D  *scat;
  TCanvas *canvas;
  TSpectrum *Tsp;
    
  TF1 *func;//fit function
  Int_t fitlinecolor;//fitlinecolor   
  
  Int_t Hist_canvas_length;
  Int_t Hist_canvas_width;
  Int_t Scat_canvas_length;
  Int_t Scat_canvas_width;

  Bool_t whichstyle;//HistStyle = 1 ,ScatStyle=0  
  
  ClassDef(zetsubou,1)//(ClassName,ClassVersionID) ClassVersionID should be >=1 

};

#if !defined(__CINT__)
ClassImp(zetsubou);
#endif

zetsubou::zetsubou()//Constructor
{
  currentfname="";
  hist=0;
  scat=0;
  tr=0;
  func=0;
  canvas=0;
  Tsp=0;
  fitlinecolor=3;//1:black 2:red 3:bright green 4:britht blue 5:yellow 6:hot pink 7:aqua 8:green 9:blue
  Hist_canvas_width = 700; 
  Hist_canvas_length = 600;
  Scat_canvas_width = 800;
  Scat_canvas_length = 700;
  SetHistStyle();
  whichstyle=HIST;//Histstyle
}

zetsubou::~zetsubou()//Destructor
{
  delete tr;
  delete hist;
  delete canvas;
  delete Tsp;
}

//setting of style in DrawHist and DrawScat

void zetsubou::SetHistStyle()
{
  Int_t fon=22;
  gROOT->SetStyle("Plain");
  //gStyle->SetOptTitle(0);//delete titile
  gStyle->SetOptTitle(1);
  //gStyle->SetOptStat(0000000);//delete Status
  gStyle->SetOptStat(1111111);
  //gStyle->SetOptFit(0000);//delete Fit Status
  gStyle->SetOptFit(1111);
  gStyle->SetPadLeftMargin(0.1);
  gStyle->SetPadRightMargin(0.31);
  gStyle->SetPadBottomMargin(0.1);
  gStyle->SetTitleXOffset(1.1);
  gStyle->SetPadGridX(1);
  gStyle->SetPadGridY(1);
  gStyle->SetHistLineWidth(3);
  gStyle->SetTitleYOffset(1.1);
  gStyle->SetTitleFont(fon);
  gStyle->SetTitleFont(fon,"Y");
  gStyle->SetLabelFont(fon);
  gStyle->SetLabelFont(fon,"Y");
  gStyle->SetTextFont(fon);
  gStyle->SetStatFont(fon);
  gStyle->SetStatFontSize(0.03);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetStatH(0.1);
  gStyle->SetStatW(0.28);
  gStyle->SetStatY(0.88);
  gStyle->SetLabelSize(0.03,"x");
  gStyle->SetLabelSize(0.03,"y");
  if(canvas)canvas->UseCurrentStyle();//update canvas 
  whichstyle = HIST;
}

void zetsubou::SetScatStyle()
{
  gROOT->SetStyle("Plain");
  //gStyle->SetOptTitle(0);//delete titile
  gStyle->SetOptTitle(1);
  //gStyle->SetOptStat(0000000);//delete Status
  gStyle->SetOptStat("em");
  gStyle->SetPadGridX(2);
  gStyle->SetPadGridY(2);
  gStyle->SetLabelSize(0.03,"x");
  gStyle->SetLabelSize(0.03,"y");
  gStyle->SetHistLineWidth(1);
  gStyle->SetPadLeftMargin(0.1);
  gStyle->SetPadRightMargin(0.1);
  gStyle->SetPadBottomMargin(0.1);
  gStyle->SetTitleXOffset(1);
  gStyle->SetTitleYOffset(1.1);
  gStyle->SetStatH(0.19);
  if(canvas)canvas->UseCurrentStyle(); 
  whichstyle = SCAT;
}

//Making a tree from datafile, Drawing 1-D Histram and 2-D Scatter diagram 
Bool_t zetsubou::MakeCanvas(const Char_t *name,Int_t canvasstyle)
{
  std::string Cnum;//canvas number
  std::string preCName = name;
  Cnum=preCName.substr(1);
  preCName.replace(0,1,"c");//I want canvas name to be "c0","c1",.....
  std::string discri;
  std::string CName;
  Char_t buf[10];
  for(Int_t i=0;i<=10;i++)//in this loop,checking sCName is "c0","c1"...,if not user input canvas name from the console
  {
	std::sprintf(buf,"c%d",i);//buf = "c0","c1"...
	discri = buf;
	if(std::string(name) == discri)//if histname or scatname is "c1","c2","c3", it should not be equall to canvasname
	{
		std::cout << "Please Input Canvas-Name:" ;
		std::cin >> CName; 
		if(!std::cin) std::exit(1);		
		break;
	}
	if(discri == preCName)
	{
		CName = preCName;		
		break;
	}
	if(i == 10)
	{
		std::cout << "Please Input Canvas-Name:" ;
		std::cin >> CName;
		if(!std::cin) std::exit(1);		
	}
  }
  //corresponding to  canvas number,canvas position change.
  Int_t canvas_positionX = 60*std::atoi(Cnum.c_str());
  Int_t canvas_positionY = 30*std::atoi(Cnum.c_str());
  if(canvasstyle == HIST) 
  canvas = new TCanvas(CName.c_str(),CName.c_str(),canvas_positionX,canvas_positionY,Hist_canvas_width,Hist_canvas_length);
  if(canvasstyle == SCAT) 
  canvas = new TCanvas(CName.c_str(),CName.c_str(),canvas_positionX,canvas_positionY,Scat_canvas_width,Scat_canvas_length);
  return true;
}

Bool_t zetsubou::MakeTree(const Char_t *filename)
{
  UShort_t  buf[4];
  std::FILE *fp;
  fp = fopen(filename, "rb");
  if(fp == NULL) 
  {
	std::cout<<"can't open file!"<<std::endl;
	return false;
  } 
  currentfname=filename;
  delete tr;  
  tr = new TTree("tr","tr");
  tr->SetDirectory(0); // not to leave a obj::TTree in memory 
  tr->Branch("E", &buf[1], "E/s");//"s" mean "a 16 bit of unsigned short integer"
  tr->Branch("dE",&buf[2], "dE/s");
  while (1) 
  {
	if (fread(buf,2,4,fp) < 1) break;
	tr->Fill();
  }
  fclose(fp);
  
  return true;
}

TH1D *zetsubou::DrawHist(const Char_t *filename,const Char_t *title,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax,Int_t cutmin,Int_t cutmax)
{
  std::string tmpfname=filename; //temporary file name
  Bool_t isTree = true; //for checking there is a same tree
  if(tmpfname != currentfname) 
  {
	isTree = MakeTree(filename); //making Tree unless same Tree already has been made
  }
  if(!isTree) return 0;

  if(hist)
  {	//string(something):converting something to string type
	if(std::string(hist->GetName()) == std::string(histname))
	{
	delete hist;//if there was a hist that had same name,it must be deleted
	hist=0;
	}
  }
  
  Char_t trdraw[100];
  Char_t trcutL[1000];
  Char_t trcutR[100];
  std::sprintf(trdraw,"%s>>%s",type,histname);//E>>hist or dE>>hist
  MakeCanvas(histname,HIST);
  if(whichstyle  == SCAT) SetHistStyle();  
  hist = new TH1D(histname,title,nbin,xmin,xmax);
  if(type == std::string("E"))
  {
	if(cutmin == 0 && cutmax == 0)
	{
		tr->Draw(trdraw,"E>=0 & dE>=0");
	}
	else
	{
		std::sprintf(trcutL,"dE >= %i &",cutmin);
		std::sprintf(trcutR,"dE < %i "  ,cutmax);
		std::strcat(trcutL,trcutR); //dE >= cutmin &dE < cutmax
		tr->Draw(trdraw,trcutL);
	}
	hist->GetXaxis()->SetTitle("E-counter");
  }
  if(type == std::string("dE"))
  {
	if(cutmin == 0 && cutmax == 0)
	{
		tr->Draw(trdraw,"");
	}
	else
	{
		std::sprintf(trcutL,"E >= %i &",cutmin);
		std::sprintf(trcutR,"E < %i "  ,cutmax);
		std::strcat(trcutL,trcutR); //E >= cutmin &E < cutmax
		tr->Draw(trdraw,trcutL);
	}
	hist->GetXaxis()->SetTitle("#DeltaE-counter");  
  }
  
  return hist;
}

//you can omit Tcut
TH1D *zetsubou::DrawHist(const Char_t *filename,const Char_t *title,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax)
{ 
  DrawHist(filename,title,type,histname,nbin,xmin,xmax,0,0);
  
  return hist;
}

//you can omit title
TH1D *zetsubou::DrawHist(const Char_t *filename,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax,Int_t cutmin,Int_t cutmax)
{ 
  std::string tmpfname=filename;
  if(tmpfname == currentfname) {}//this is necessary but I don't know why  

  DrawHist(filename,filename,type,histname,nbin,xmin,xmax,cutmin,cutmax);
  
  return hist;
}

//you can omit title and Tcut
TH1D *zetsubou::DrawHist(const Char_t *filename,const Char_t *type,const Char_t *histname,Int_t nbin,Double_t xmin,Double_t xmax)
{
  DrawHist(filename,filename,type,histname,nbin,xmin,xmax,0,0);
 
  return hist;
}

//for oblique cut 
TH1D *zetsubou::HistCut(Double_t X1,Double_t Y1,Double_t X2,Double_t Y2,Double_t X3,Double_t Y3,Double_t X4,Double_t Y4)
{
  if(!hist) return 0;

  Char_t cutline1[100];//under border
  Char_t cutline2[100];//right border
  Char_t cutline3[100];//top border
  Char_t cutline4[100];//left border
  if(X1 == X2 || X3 == X4) 
  {
	std::cout << "can't cut such a way!" << std::endl;
	return 0;
  }
  std::sprintf(cutline1,"dE >= %g*E +  %g",(Y2-Y1) / (X2-X1) ,Y1 - ((Y2-Y1) / (X2-X1))*X1 );
  std::sprintf(cutline3,"dE <= %g*E + %g",(Y3-Y4) / (X3-X4) ,Y4 - ((Y3-Y4) / (X3-X4))*X4 );
  
  if(X2 == X3){std::sprintf(cutline2,"E <= %g",X3);} 
  if(X2 < X3) {std::sprintf(cutline2,"dE >= %g*E + %g",(Y3-Y2) / (X3-X2) ,Y2 - ((Y3-Y2) / (X3-X2))*X2 );}
  if(X2 > X3) {std::sprintf(cutline2,"dE <= %g*E + %g",(Y3-Y2) / (X3-X2) ,Y2 - ((Y3-Y2) / (X3-X2))*X2 );}
  
  if(X1 == X4){std::sprintf(cutline4,"E >= %g",X1);}
  if(X1 < X4) {std::sprintf(cutline4,"dE <= %g*E + %g",(Y4-Y1) / (X4-X1) ,Y1 - ((Y4-Y1) / (X4-X1))*X1 );}
  if(X1 > X4) {std::sprintf(cutline4,"dE >= %g*E + %g",(Y4-Y1) / (X4-X1) ,Y1 - ((Y4-Y1) / (X4-X1))*X1 );}
  
  std::string  totalcut;
  totalcut = cutline1 +std::string("&&")+ cutline2 +"&&"+ cutline3+ "%%" + cutline4;
  std::cout << "----- border line -----" << std::endl;
  std::cout << "under: " << cutline1 << std::endl;
  std::cout << "right: " << cutline2 << std::endl;
  std::cout << "top: " << cutline3 << std::endl;
  std::cout << "left: "  << cutline4 << std::endl;
  
  Char_t trdraw[100];
  std::sprintf(trdraw,"E >> %s",(hist->GetName()));
  tr->Draw(trdraw,totalcut.c_str());
  return hist;	
}

TH2D *zetsubou::DrawScat(const Char_t *filename,const Char_t *title,const Char_t *scatname)
{
  std::string tmpfname=filename;//temporary file name
  Bool_t isTree = true;//for checking there is same tree
  if(tmpfname != currentfname)
  {
	isTree = MakeTree(filename);
  }				//making Tree unless same Tree already has been made
  if(!isTree) return 0;
  
  if(scat) 
  {
	if(std::string (scat->GetName()) == std::string (scatname) ) 
	{  
	delete scat;
	scat=0;
	}
  }//if there was a scat that had same name,it must be deleted
  MakeCanvas(scatname,SCAT);
  if(whichstyle == HIST) SetScatStyle();
  canvas->Divide(2,2);
  
  canvas->cd(1);
  hist = new TH1D("px","",128,0,4095);
  tr->Draw("E>>px","");
  
  canvas->cd(2);
  hist = new TH1D("py","", 128 ,0,4095);
  tr->Draw("dE>>py","");
  
  canvas->cd(3);
  Char_t trdraw[100]; 
  std::sprintf(trdraw,"dE:E >>%s",scatname);
  scat = new TH2D(scatname,title,1024,0,4095,1024,0,4095);
  tr->Draw(trdraw,"");
  scat->GetXaxis()->SetTitle("E");
  scat->GetYaxis()->SetTitle("#DeltaE");
    
  return scat;
}

//you can omit title
TH2D *zetsubou::DrawScat(const Char_t *filename,const Char_t *scatname)
{
  DrawScat(filename,filename,scatname);

  return scat;
}

TH2D *zetsubou::DrawScat(const Char_t *filename,const Char_t *title,const Char_t *scatname,Int_t xmin, Int_t xmax,Int_t ymin,Int_t ymax)
{
  std::string tmpfname=filename;
  Bool_t isTree = true;
  if(tmpfname != currentfname)
  {
	isTree = MakeTree(filename);
  }				
  if(!isTree) return 0;
  
  MakeCanvas(scatname,SCAT);  
  if(whichstyle == HIST) SetScatStyle(); 
  
  if(scat) 
  {
	if(std::string (scat->GetName()) == std::string (scatname))
	{
	delete scat;
	scat=0;
	}
  } 
  
  Char_t trcut[100]; 
  std::sprintf(trcut,"E>=%i & E<%i & dE>=%i & dE<%i",xmin,xmax,ymin,ymax);
  canvas->Divide(2,2);
  
  canvas->cd(1);
  tr->Draw("E",trcut);
  
  canvas->cd(2);
  tr->Draw("dE",trcut);
  canvas->UseCurrentStyle();
  
  canvas->cd(3);
  Char_t trdraw[100]; 
  std::sprintf(trdraw,"dE:E >>%s",scatname);
  scat = new TH2D(scatname,title,1024,0,4095,1024,0,4095);
  tr->Draw(trdraw,"");
  scat->GetXaxis()->SetTitle("E");
  scat->GetYaxis()->SetTitle("#DeltaE");

  canvas->cd(4);  
  tr->Draw("dE:E",trcut);
  scat->GetXaxis()->SetTitle("E");
  scat->GetYaxis()->SetTitle("#DeltaE");
  return scat;
}

//you can omit title
TH2D *zetsubou::DrawScat(const Char_t *filename,const Char_t *scatname,Int_t xmin, Int_t xmax,Int_t ymin,Int_t ymax)
{
  DrawScat(filename,filename,scatname,xmin,xmax,ymin,ymax);

  return scat;
}

//definision of fit function

//expression of single gaussian 
Double_t singlegaussian(Double_t *x, Double_t *par)
{
  Double_t fitval = par[0]/(sqrt(2.*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
  return fitval;
}

//expression of single gaussian + linear
Double_t singlegaussian_linear(Double_t *x, Double_t *par)
{
 Double_t fitval = par[0]/(sqrt(2.*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2])
                       + par[3]*x[0]+par[4];
 return fitval;
}

//expresion of double gaussian 
Double_t doublegaussian(Double_t *x, Double_t *par)
{
   Double_t fitval = par[0]/(sqrt(2.*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2])
                         + par[3]/(sqrt(2.*TMath::Pi())*par[5])*TMath::Exp(-(x[0]-par[4])*(x[0]-par[4])/2./par[5]/par[5]);
  return fitval;
}

//expresion of double gaussian +linear
Double_t doublegaussian_linear(Double_t *x, Double_t *par)
{
  Double_t fitval = par[0]/(sqrt(2.*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2])
                        + par[3]/(sqrt(2.*TMath::Pi())*par[5])*TMath::Exp(-(x[0]-par[4])*(x[0]-par[4])/2./par[5]/par[5])
                        + par[6]*x[0]+par[7];
  return fitval;
}

//expresion of triple gaussian
Double_t triplegaussian(Double_t *x, Double_t *par)
{
  Double_t fitval =par[0]/(sqrt(2.*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2])
                        +par[3]/(sqrt(2.*TMath::Pi())*par[5])*TMath::Exp(-(x[0]-par[4])*(x[0]-par[4])/2./par[5]/par[5])
                        +par[6]/(sqrt(2.*TMath::Pi())*par[8])*TMath::Exp(-(x[0]-par[7])*(x[0]-par[7])/2./par[8]/par[8]);
  return fitval;
}

//expression of triple gaussian + linear
Double_t triplegaussian_linear(Double_t *x, Double_t *par)
{
  Double_t fitval =par[0]/(sqrt(2.*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2])
                        +par[3]/(sqrt(2.*TMath::Pi())*par[5])*TMath::Exp(-(x[0]-par[4])*(x[0]-par[4])/2./par[5]/par[5])
                        +par[6]/(sqrt(2.*TMath::Pi())*par[8])*TMath::Exp(-(x[0]-par[7])*(x[0]-par[7])/2./par[8]/par[8])
                        +par[9]*x[0]+par[10];
  return fitval;
}

//Convoluted Landau and Gaussian Fitting Function
Double_t langaufun(Double_t *x, Double_t *par)
{
  // Numeric constants
  Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
  Double_t mpshift  = -0.22278298;       // Landau maximum location
  
  // Control constants
  Double_t np = 100.0;      // number of convolution steps
  Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
  
  // Variables
  Double_t xx;
  Double_t mpc;
  Double_t fland;
  Double_t sum = 0.0;
  Double_t xlow,xupp;
  Double_t step;
  Double_t i;

  // MP shift correction
  mpc = par[1] - mpshift * par[0]; 

  // Range of convolution integral
  xlow = x[0] - sc * par[3];
  xupp = x[0] + sc * par[3];

  step = (xupp-xlow) / np;

  // Convolution integral of Landau and Gaussian by sum
  for(i=1.0; i<=np/2; i++) 
  {
  xx = xlow + (i-.5) * step;
  fland = TMath::Landau(xx,mpc,par[0]) / par[0];
  sum += fland * TMath::Gaus(x[0],xx,par[3]);

  xx = xupp - (i-.5) * step;
  fland = TMath::Landau(xx,mpc,par[0]) / par[0];
  sum += fland * TMath::Gaus(x[0],xx,par[3]);
  }

  return (par[2] * step * sum * invsq2pi / par[3]);
}

//single gaussian 
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax)
{
  if(!hist) return 0; 
  func = new TF1("func",singlegaussian,xmin,xmax,3);//fitfunctionname,expression,area,number of parameters
  func->SetLineColor(fitlinecolor);
  func->SetParNames("const","mean","sigma");
  Double_t ms,mmax;
  ms= hist->GetMaximumBin()*hist->GetBinWidth(2);//position of maximumbin in the canvas
  mmax=ms + canvas->GetUxmin(); //absolute position 
  func->SetParameters(100,mmax,10);//set initial parameter
  func->SetParLimits(2,0,xmax-xmin+10);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event=" << func->GetParameter("const")  / hist->GetBinWidth(2) << std::endl;
  std::cout << "event error=" << func->GetParError(0) / hist->GetBinWidth(2) <<std::endl;
  std::cout << "mean=" << func->GetParameter("mean") << std::endl;
  std::cout << "mean error=" << func->GetParError(1) << std::endl;
 
  return func;
}

//sigle gaussian 
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax,Double_t height,Double_t mean,Double_t sigma)
{
  if(!hist) return 0; 
  func = new TF1("func",singlegaussian,xmin,xmax,3);
  func->SetLineColor(fitlinecolor);
  func->SetParameters(height,mean,sigma);
  func->SetParNames("const","mean","sigma");
  func->SetParLimits(1,xmin,xmax);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event=" << func->GetParameter("const") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event error=" << func->GetParError(0) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean=" << func->GetParameter("mean") << std::endl;
  std::cout << "mean error=" << func->GetParError(1) << std::endl;  
  return func;
}

//single gaussian + linear
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax,Double_t height,Double_t mean,Double_t sigma,Double_t katamuki ,Double_t seppen)
{
  if(!hist) return 0;   
  func = new TF1("func",singlegaussian_linear,xmin,xmax,5);
  func->SetLineColor(fitlinecolor);
  func->SetParameters(height,mean,sigma,katamuki,seppen);
  func->SetParNames("const","mean","sigma","slope","intercept");
  func->SetParLimits(1,xmin,xmax);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event=" << func->GetParameter("const") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event error=" << func->GetParError(0) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean=" << func->GetParameter("mean") << std::endl;
  std::cout << "mean error=" << func->GetParError(1) << std::endl;  
  return func;
}

//double gaussian
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2)
{
  if(!hist) return 0;   
  TF1 *func = new TF1("func",doublegaussian,xmin,xmax,6);
  func->SetLineColor(fitlinecolor);
  func->SetParameters(height1,mean1,sigma1,height2,mean2,sigma2);
  func->SetParNames("const1","mean1","sigma1","const2","mean2","sigma2");
  func->SetParLimits(1,xmin,xmax);
  func->SetParLimits(4,xmin,xmax);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event1=" << func->GetParameter("const1") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event1 error=" << func->GetParError(0) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean1=" << func->GetParameter("mean1") << std::endl;
  std::cout << "mean1 error=" << func->GetParError(1) << std::endl;
  std::cout << "event2=" << func->GetParameter("const2") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event2 error=" << func->GetParError(3) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean2=" << func->GetParameter("mean2") << std::endl;
  std::cout << "mean2 error=" << func->GetParError(4) << std::endl;
  
  return func;
}

//double gaussian + linear
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2,Double_t katamuki ,Double_t seppen)
{
  if(!hist) return 0;  
  func = new TF1("func",doublegaussian_linear,xmin,xmax,8);
  func->SetLineColor(fitlinecolor);
  func->SetParameters(height1,mean1,sigma1,height2,mean2,sigma2,katamuki,seppen);
  func->SetParNames("const1","mean1","sigma1","const2","mean2","sigma2","slope","intercept");
  func->SetParLimits(1,xmin,xmax);
  func->SetParLimits(4,xmin,xmax);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event1=" << func->GetParameter("const1") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event1 error=" << func->GetParError(0) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean1="<< func->GetParameter("mean1") << std::endl;
  std::cout << "mean1 error=" << func->GetParError(1) << std::endl;
  std::cout << "event2=" << func->GetParameter("const2") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event2 error=" << func->GetParError(3) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean2=" << func->GetParameter("mean2") << std::endl;
  std::cout << "mean2 error=" << func->GetParError(4)<< std::endl;
  return func;
}

//triple gaussian
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2,Double_t height3,Double_t mean3,Double_t sigma3)
{
  if(!hist) return 0;  
  func = new TF1("func",triplegaussian,xmin,xmax,9);
  func->SetLineColor(fitlinecolor);
  func->SetParameters(height1,mean1,sigma1,height2,mean2,sigma2,height3,mean3,sigma3);
  func->SetParNames("const1","mean1","sigma1","const2","mean2","sigma2","const3","mean3","sigma3");
  func->SetParLimits(1,xmin,xmax);
  func->SetParLimits(4,xmin,xmax);
  func->SetParLimits(7,xmin,xmax);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event1=" << func->GetParameter("const1") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event1 error=" << func->GetParError(0) / hist->GetBinWidth(2)<< std::endl;
  std::cout << "mean1=" << func->GetParameter("mean1") << std::endl;
  std::cout << "mean1 error=" << func->GetParError(1) << std::endl;
  std::cout << "event2=" << func->GetParameter("const2") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event2 error=" << func->GetParError(3) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean2=" << func->GetParameter("mean2") << std::endl;
  std::cout << "mean2 error=" << func->GetParError(4) << std::endl;
  std::cout << "event3=" << func->GetParameter("const3") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event3 error=" << func->GetParError(6) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean3=" << func->GetParameter("mean3") << std::endl;
  std::cout << "mean3 error=" << func->GetParError(7) << std::endl;
  return func;
}

//triple gaussian + linear
TF1 *zetsubou::gfit(Double_t xmin,Double_t xmax,Double_t height1,Double_t mean1,Double_t sigma1,Double_t height2, Double_t mean2, Double_t sigma2,Double_t height3,Double_t mean3,Double_t sigma3,Double_t katamuki ,Double_t seppen)
{
  if(!hist) return 0;
  func = new TF1("func",triplegaussian_linear,xmin,xmax,11);
  func->SetLineColor(fitlinecolor);
  func->SetParameters(height1,mean1,sigma1,height2,mean2,sigma2,height3,mean3,sigma3,katamuki,seppen);
  func->SetParNames("const1","mean1","sigma1","const2","mean2","sigma2","const3","mean3","sigma3","slope","intercept");
  func->SetParLimits(1,xmin,xmax);
  func->SetParLimits(4,xmin,xmax);
  func->SetParLimits(7,xmin,xmax);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" << func->GetNDF() << std::endl;
  std::cout << "event1=" << func->GetParameter("const1") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event1 error=" << func->GetParError(0) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean1=" << func->GetParameter("mean1") << std::endl;
  std::cout << "mean1 error=" << func->GetParError(1)<< std::endl;
  std::cout << "event2=" << func->GetParameter("const2") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event2 error=" << func->GetParError(3) / hist->GetBinWidth(2)<< std::endl;
  std::cout << "mean2=" << func->GetParameter("mean2") << std::endl;
  std::cout << "mean2 error=" << func->GetParError(4) << std::endl;
  std::cout << "event3=" << func->GetParameter("const3") / hist->GetBinWidth(2) << std::endl;
  std::cout << "event3 error=" << func->GetParError(6) / hist->GetBinWidth(2) << std::endl;
  std::cout << "mean3=" << func->GetParameter("mean3") << std::endl;
  std::cout << "mean3 error=" << func->GetParError(7) << std::endl;
  return func;
}

//Convoluted Landau and Gaussian Fitting 
TF1 *zetsubou::lgfit(Double_t xmin,Double_t xmax,Double_t width,Double_t mp,Double_t area,Double_t gsigma)
{
  if(!hist) return 0;
  std::cout << "Now fitting......." << std::endl;
  func = new TF1("func",langaufun,xmin,xmax,4);
  func->SetParNames("Width","MP","Area","GSigma");
  func->SetParameters(width,mp,area,gsigma);
  func->SetLineColor(fitlinecolor);
  hist->Fit("func","","",xmin,xmax);
  hist->Draw();
  std::cout << "Chisquare/NDF=" << func->GetChisquare() / func->GetNDF() << "=" << func->GetChisquare() << "/" <<func->GetNDF() << std::endl;
  return func;
}

TSpectrum *zetsubou::Tspfit(Double_t xmin,Double_t xmax)
{
  if(!hist) return 0;
  std::cout << "Now fitting..." << std::endl;
  Tsp = new TSpectrum(3,1); //(maximum number of peaks, resolution of the neighboring peaks)
  Tsp->Search(hist,2,"",0.10); //(hist,sigma of searched peaks,option,threshold)
  Int_t np =0; 
  Float_t *xpeaks = Tsp->GetPositionX();
  Float_t x[3];
  for(Int_t i=0;i <=2;i++)
  {
	if ( xpeaks[i] > xmin && xpeaks[i] < xmax)
	{
	x[np] = xpeaks[i];
	np++;
	}
  }		
  if(np == 1) gfit(xmin,xmax,100,x[0],10);
  if(np == 2) gfit(xmin,xmax,100,x[0],10,100,x[1],10);
  if(np == 3) gfit(xmin,xmax,100,x[0],10,100,x[1],10,100,x[2],10);
  return Tsp;
}

TSpectrum  *zetsubou::Tspfit(Double_t xmin,Double_t xmax,Double_t threshold)
{
  if(!hist) return 0;
  std::cout << "Now fitting..." << std::endl;
  Tsp = new TSpectrum(3,1); //(maximum number of peaks, resolution of the neighboring peaks)
  Tsp->Search(hist,2,"",threshold); //(hist,sigma of searched peaks,option,threshold)
  Int_t np =0;
  Float_t *xpeaks = Tsp->GetPositionX();
  Float_t x[3];
  for(Int_t i=0;i <=2;i++)
  {
	if ( xpeaks[i] > xmin && xpeaks[i] < xmax)
	{
	x[np] = xpeaks[i];
	np++;
	}
  }		
  if(np == 1) gfit(xmin,xmax,100,x[0],10);
  if(np == 2) gfit(xmin,xmax,100,x[0],10,100,x[1],10);
  if(np == 3) gfit(xmin,xmax,100,x[0],10,100,x[1],10,100,x[2],10);
  return Tsp;
}

TSpectrum *zetsubou::Tspfit(Double_t xmin,Double_t xmax,Double_t threshold,Double_t resolution)
{
  if(!hist) return 0;
  std::cout << "Now fitting..." << std::endl;
  Tsp = new TSpectrum(3,resolution); //(maximum number of peaks, resolution of the neighboring peaks)
  Tsp->Search(hist,2,"",threshold); //(hist,sigma of searched peaks,option,threshold)
  Int_t np =0;
  Float_t *xpeaks = Tsp->GetPositionX();
  Float_t x[3];
  for(Int_t i=0;i <=2;i++)
  {
	if ( xpeaks[i] > xmin && xpeaks[i] < xmax)
	{
	x[np] = xpeaks[i];
	np++;
	}
  }		
  if(np == 1) gfit(xmin,xmax,100,x[0],10);
  if(np == 2) gfit(xmin,xmax,100,x[0],10,100,x[1],10);
  if(np == 3) gfit(xmin,xmax,100,x[0],10,100,x[1],10,100,x[2],10);
  return Tsp;
}