#include "PeakFit.h"

void PeakFit::Fit( TH1F *tdc ){

  TString str_title = "PMT ";str_title += mcp_ch / 4;str_title += " ch# " ;str_title += mcp_ch % 4;
  gPad->SetTitle(str_title);
  tdc->SetTitle(str_title);


  /** Search fit parameter **/
  cout << "Searching fit parameter start!" << endl;

  vector<double> p_x , p_y , p_sigma;
  int cont , pre_count = 0 , pre_max = 0 , pre_min = 100000;
  bool flag_up = 0 , flag_down = 0 , flag_data = 0;
  int p_min_x = 0, p_max_x = 0;

  for(Int_t bin = 1;bin < (range_max - range_min);bin++){

    cont = (int)tdc->GetBinContent(bin);
    
    if( cont > min_th ){

       if( ( cont > pre_min ) && ( cont > pre_count) && ( cont > pre_max ) ){
       //cout << "bin:" << x << " , cont:" << cont << endl;

          pre_max = cont;
          p_max_x = bin;
          //cout << "Over cont , Bin : " << range_min + bin << " , pre_max = " << pre_max << endl;

       }else{

          cout << "cont:" << cont << " & pre:" << pre_count << " , " ;
          if( cont > pre_count ){
             pre_min = cont;
             pre_max = min_th;
             //cout << "Not over cont , Bin : " << range_min + bin << " , cont:" << cont << " , pre_min = " << pre_min << endl;
          }

       }

   }


   //make slope flag
   if( !flag_up ){

     if( (pre_max - pre_min) > up_th ){
        flag_up = 1;
        p_min_x = bin;
        //cout << "Check : Up flag ! " << endl;
     }

   }else{

     if( ( pre_max - cont ) > down_th){
       flag_down = 1;
       //cout << "Check : Down flag ! " << endl;
     }

   }

   pre_count = cont;

   //check flag & set parameter
   if(flag_up && flag_down){

      if(!flag_data){
         p_x.push_back( range_min + p_max_x );
         p_y.push_back( (double)pre_max * 0.8);
         p_sigma.push_back( (double)bin - p_min_x);
         
         cout << "peak_x : " << range_min + p_max_x << " , peak_cont : " << pre_max << " , sigma : " << bin - p_min_x << endl;

         flag_data = 1;
         pre_min = cont;
      }else{

            cout << "Check ," ;
         if( ( cont - pre_min ) > 0 ){
            flag_up = 0;
            flag_down = 0;
            flag_data = 0;

            cout << "All flag are off ! " << endl;
         }
      }
   }
    
  }

  //BG parameter 
  if(bg_fit){
    p_y.push_back(min_th);
    p_x.push_back(tdc->GetMean());
    p_sigma.push_back(tdc->GetRMS());
  }



  //Fit TDC peaks 
  const int p_num =p_x.size();

  TString str_fit = "";
  for(int pn = 0;pn < p_num;pn++){
    str_fit += "gaus(";str_fit += pn * 3;str_fit += ")";
    str_fit += ( pn == ( (int)p_x.size() - 1) ) ?  "" : " + ";
  }
  cout << "Fit Function # : " << p_num << " , " << str_fit << endl;

  TF1 *fit = new TF1("fit",str_fit);

  TString par_name = "";
  for(int pn = 0;pn < p_num;pn++){

     fit->SetParameter( 3 * pn , p_y[pn] );
     cout << "para(" << 3*pn << ") : " << p_y[pn] << " , " ;
     par_name = "Height(";par_name += pn;par_name += ")";
     fit->SetParName(3 * pn , par_name);

     fit->SetParameter( 3 * pn + 1, p_x[pn] );
     cout << "para(" << 3*pn+1 << ") : " << p_x[pn] << " , " ;
     par_name = "Peak(";par_name += pn;par_name += ")";
     fit->SetParName(3 * pn + 1, par_name);

     fit->SetParameter( 3 * pn + 2, p_sigma[pn] );
     cout << "para(" << 3*pn+2 << ") : " << p_sigma[pn] << endl ;
     par_name = "#sigma(";par_name += pn;par_name += ")";
     fit->SetParName(3 * pn + 2, par_name);

  }

  tdc->Fit("fit","revl","goff",range_min,range_max);
  tdc->Draw();


  //Draw all plot
  TF1 *f_fit[p_num];
  double para[3 * p_num]; 
  fit->GetParameters(para);
  //cout << "error" << fit->GetParError(0) << endl;

  cout << endl;
  for(int pn = 0;pn < p_num;pn++){
     TString gaus;
     gaus = "exp(-(x-(";gaus += para[3 * pn + 1];
     gaus += "))**2/(2.*";gaus += para[3 * pn + 2];
     gaus += "**2))*";gaus += para[3 * pn];

     TString str_f = "f";str_f += pn;
     f_fit[pn] = new TF1(str_f,gaus,range_min,range_max);

     int col = (pn < 3) ? pn+2 : pn + 3;
     f_fit[pn]->SetLineColor(col);
     f_fit[pn]->SetLineWidth(1);

     f_fit[pn]->Draw("same");

     int pp;
     if(pn == p_num - 1){
        pp = -1;
     }else{
        pp = pn + 1;
     }

     cout << mcp_ch / 4 << " " << mcp_ch % 4 + 1 << " " << pp << " " << para[3 * pn + 2 ] << " " << fit->GetParError(3 * pn + 2) << " " << para[3 * pn] << " " << fit->GetParError(3 * pn) << " " << para[3 * pn + 1] << " " << fit->GetParError(3 * pn + 1) << " " << endl;

  }

  cout << endl;
}