//for wave form analysis
//made by arita 2011.12.27

#include <iostream> //{{{
#include <fstream>
#include <string>
#include <cmath>
#include <vector>
#include "TApplication.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TCut.h"
#include "TString.h"
#include "TBranch.h"
#include "TGraph.h"
#include "TChain.h"
#include "TTree.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TGraphErrors.h"
#include "TGraph2D.h"
#include "TPaveStats.h"
#include "TText.h"
#include "TNtuple.h"
#include "TRandom.h"
#include <sstream>
#include <stdlib.h>

#include "WaveForm.h"
//}}}


//Constructor
WaveForm::WaveForm(void){ //{{{

   n_accumulate = 0;
   n_sample = 0;

   n_acc_period = 0;
   stck = 0;

   b_acc = false;

   wf_m_type = 1;
   wf_l_col = 1;
   wf_m_size = 1;
   wf_xtitle = "X";
   wf_ytitle = "Y";

   acc_period_sum = 0.; //for approximationa
   for( int i = 0;i < 10000;i++)acc_period_stock[i] = 0.;

   p_nsample = 0;

   n_xpoint[0] = 0;
   n_xpoint[1] = 0;

}; //}}}



//-- Private Function --//


//-- Public Function --//
//General
void WaveForm::SetWaveform( float *wf ){ //{{{

   int n =  sizeof wf / sizeof wf[0];


   if( n != 0 && n != n_sample){
      std::cout << "<Warning in Waveform::SetWaveform,  nuber of sample is changed :" << n << std::endl; 
   }

   if( n_sample < n){
      for( int i = n_sample;i < n;i++){
         adc.push_back( 0. );
      }
      n_sample = n;
   }

   for( int i = 0;i < n_sample;i++){
      adc[i] = wf[i];
   }

} //}}}


void WaveForm::SetWaveform( int id ,  float fl_adc ){ //{{{

   for( int i = n_sample;i <= id;i++){
      adc.push_back( 0. );
      time2f.push_back( 0. );
      adc2f.push_back( 0. );

      n_sample++;
   }

   adc[id] = fl_adc;

   time2f[id] = (float)id;
   adc2f[id] = fl_adc;

} //}}}


void WaveForm::SetWaveform( int id , float fl_time ,  float fl_adc ){ //{{{

   for( int i = n_sample;i <= id;i++){
      adc.push_back( 0. );
      time2f.push_back( 0. );
      adc2f.push_back( 0. );

      n_sample++;
   }

   adc[id] = fl_adc;

   time2f[id] = fl_time;
   adc2f[id] = fl_adc;

} //}}}


int WaveForm::ShiftWaveform( float time_shift , float adc_shift ){ //{{{

   for( int i = 0;i < n_sample;i++){
      time2f[i] += time_shift;
      adc2f[i] += adc_shift;
   }

   return 0;
} //}}}


//Get general values
int WaveForm::GetWaveform( std::vector<float>*rt_time , std::vector<float>*rt_adc ){ //{{{

   for( int i = 0;i < n_sample;i++){
      rt_time->push_back( time2f[i] );
      rt_adc->push_back( adc2f[i] );
   }

   return n_sample;
} //}}}


float WaveForm::GetTime( int id ){ //{{{
   return time2f[ id ];
} //}}}


float WaveForm::GetAdc( int id ){ //{{{
   return adc2f[ id ];
} //}}}


int WaveForm::GetNsample( void ){ //{{{
   return n_sample;
} //}}}


float WaveForm::GetMeanADC( void ){ //{{{

   float sum_adc = 0 ;

   for( int s = 0;s < n_sample;s++){
      sum_adc += adc2f[s];
   }

   return ( sum_adc / (float)n_sample );
} //}}}


float WaveForm::GetRMSADC( void ){ //{{{

   float sum_rms = 0.;
   float adc_mean = GetMeanADC();

   for( int s = 0;s < n_sample;s++){
      sum_rms += pow( adc2f[s] - adc_mean , 2 );
   }

   return ( sum_rms / (float)n_sample );
} //}}}


float WaveForm::Integral( int start_id , int end_id ){ //{{{
   float sum = 0.;

   for( int s = start_id; s <= end_id;s++){
      sum += adc2f[s];
   }

   return sum;
} //}}}


int WaveForm::SortTime( void ){ //{{{

   std::vector<int> id_list;
   std::vector<float> tmp_time;
   std::vector<float> tmp_adc;
   //initialize 
   for( int i = 0;i < n_sample;i++){
      id_list.push_back( i );
      tmp_time.push_back( time2f[ i ] );
      tmp_adc.push_back( adc2f[ i ] );
   }

   //bubble sort
   for( int i = 0;i < n_sample;i++){
      for( int j = n_sample-1;j > i;j--){
         if( time2f[ j ] < time2f[ j - 1 ] ){
            float tmp1 = time2f[j];

            time2f[ j ] = time2f[ j - 1];
            time2f[ j - 1 ] = tmp1;

            id_list[ j ] = j - 1;
            id_list[ j - 1 ] = j;
         }
      }
   }

   //fill adc values
   for( int s = 0;s < n_sample;s++)adc2f[ s ] = tmp_adc[ id_list[ s ] ];

   return 1;
} //}}}


//Plotting
int WaveForm::SetMarkerStyle( int m_style ){ //{{{
   wf_m_type = m_style;
   return 0;
} //}}}


int WaveForm::SetLineColor( int col ){ //{{{

   wf_l_col = col;

   return 0;
} //}}}


int WaveForm::SetMarkerSize( float m_size ){ //{{{

   wf_m_size = m_size;

   return 0;
} //}}}


int WaveForm::SetTitle( TString wf_t ){ //{{{

   wf_title = wf_t;

   return 0;
} //}}}


int WaveForm::SetTitleX( TString wf_t ){ //{{{

   wf_xtitle = wf_t;

   return 0;
} //}}}


int WaveForm::SetTitleY( TString wf_t ){ //{{{

   wf_ytitle = wf_t;

   return 0;
} //}}}


int WaveForm::Draw( TString option , int col ){ //{{{

   TGraph *g_wf = new TGraph();
   for( int s = 0;s < n_sample;s++){
      g_wf->SetPoint( s , time2f[s] , adc2f[s] );
   }

   gPad->SetGrid();
   g_wf->SetMarkerStyle( wf_m_type );
   g_wf->SetMarkerSize( wf_m_size );
   g_wf->SetLineColor( wf_l_col );
   g_wf->SetTitle( wf_title );
   g_wf->GetYaxis()->SetTitle( wf_ytitle );
   g_wf->GetXaxis()->SetTitle( wf_xtitle );
   g_wf->Draw( option );

   return 0;
} //}}}


void WaveForm::Clear( void ){ //{{{

   for( int i = 0;i < n_sample;i++){
      adc[i] = 0.;
   }

   n_accumulate = 0;
   n_sample = 0;

   b_acc = false;
} //}}}


//Fitting
int WaveForm::FitSineWave( void ){ //{{{

   TF1 *f_sin = new TF1( "sine" , "[0] * sin( [1] * x + [2] ) + [3] "); 
   f_sin->SetParameter( 0 , fit_para[0] );
   //f_sin->SetParLimits( 0 , 0.9 , 1.1 );
   f_sin->SetParameter( 1 , fit_para[1] );
   f_sin->SetParameter( 2 , fit_para[2] );
   f_sin->SetParameter( 3 , fit_para[3] );

   TGraphErrors *graph = new TGraphErrors();
   for( int s = 0;s < n_sample;s++){
      graph->SetPoint( s , time2f[s] , adc2f[s] );
      graph->SetPointError( s , 10.e-12 , 0.1 );
   }
   graph->Fit( f_sin , "rel" , "goff" );

   fit_para[0] = f_sin->GetParameter(0);
   fit_para[1] = f_sin->GetParameter(1);
   fit_para[2] = f_sin->GetParameter(2);
   fit_para[3] = f_sin->GetParameter(3);

   fit_para_error[0] = f_sin->GetParError(0);
   fit_para_error[1] = f_sin->GetParError(1);
   fit_para_error[2] = f_sin->GetParError(2);
   fit_para_error[3] = f_sin->GetParError(3);

   return 1;
} //}}}


int WaveForm::SetFitParameter( int id , float para ){ //{{{
   fit_para[id] = para;
   return 1;
} //}}}


float WaveForm::GetFitParameter( int id ){ //{{{
   return fit_para[id];
} //}}}


float WaveForm::GetFitParError( int id ){ //{{{
   return fit_para_error[id];
} //}}}



//cross point analysis
int WaveForm::AccumulateXpoint( float baseline ){ //{{{

   float zero_adc = baseline;

   //initialize Accumulation
   if( !b_acc ){

      for( int i = 0;i < n_sample;i++){
         accum_xpoint.push_back( 0 );
         accum_xpoint_fall.push_back( 0 );
         accum_xpoint_rise.push_back( 0 );
      }

      b_acc = true;
   }

   n_xpoint[0] = 0;
   n_xpoint[1] = 0;

   //scan sample
   for( int s = 0;s < ( n_sample - 1 );s++){
      if( adc[s] > zero_adc && adc[s+1] < zero_adc ){ //falling edge 
         accum_xpoint_fall[s]++;
         accum_xpoint[s]++;

         n_xpoint[0]++;
         n_accumulate++;
      }else if( adc[s] < zero_adc && adc[s+1] > zero_adc ){ //rising edge
         accum_xpoint_rise[s]++;
         accum_xpoint[s]++;

         n_xpoint[1]++;
         n_accumulate++;
      }
   }

   
   return 0;
} //}}}


int WaveForm::GetAccumContent( int sample , int mode ){ //{{{ 
   int acc;

   if( mode == 0){
      acc = accum_xpoint[sample];
   }else if( mode < 0 ){
      acc = accum_xpoint_fall[sample];
   }else if( mode > 0){
      acc = accum_xpoint_rise[sample];
   }

   return acc;
} //}}}


float WaveForm::GetAccumMean( int mode ){ //{{{

   float sum = 0.;
   for( int s = 0;s < n_sample;s++){
      if( mode == 0){
         sum += accum_xpoint[s];
      }else if( mode < 0 ){
         sum += accum_xpoint_fall[s];
      }else if( mode > 0){
         sum += accum_xpoint_rise[s];
      }
   }

   return  sum / n_sample;
} //}}}


float WaveForm::GetAccumRMS( int mode ){ //{{{

   float mean =  GetAccumMean( mode );

   float sum = 0.;
   for( int s = 0;s < n_sample;s++){
      if( mode == 0){
         sum += pow( accum_xpoint[s] - mean , 2 );
      }else if( mode < 0 ){
         sum += pow( accum_xpoint_fall[s] - mean , 2 );
      }else if( mode > 0){
         sum += pow( accum_xpoint_rise[s] - mean , 2 );
      }
   }

   sum /= n_sample;

   return sqrt( sum );

} //}}}


int WaveForm::GetXpoint( int mode ){ //{{{ mode 0:fall , 1:rise
   return n_xpoint[mode];
} //}}}

int WaveForm::AccumulatePeriodX( float baseline ){ //{{{

   float zero_adc = baseline;

   float x_r[1000] = { 0. } , x_f[1000] = { 0. };
   int n_r = 0 , n_f = 0;

   float p_f = 0. , p_r = 0.;

   //scan sample
   for( int s = 0;s < ( n_sample - 1 );s++){
      if( adc[s] > zero_adc && adc[s+1] < zero_adc ){ //falling edge 
         float s1 = (float)s * 2 , s2 = (float)(s+1) * 2;
         x_f[n_f] = ( s1 + ( ( s2 - s1 ) * ( zero_adc - adc[s]) ) / ( adc[s+1] - adc[s] ) );

         if( n_f > 0){
            p_f = x_f[ n_f ] - x_f[ n_f - 1 ];
            acc_period_sum += p_f;

            if( n_acc_period < 10000)acc_period_stock[ stck++ ] = p_f;

            n_acc_period++;
         }


         n_f++;
      }else if( adc[s] < zero_adc && adc[s+1] > zero_adc ){ //rising edge
         float s1 = (float)s * 2 , s2 = (float)(s+1) * 2;
         x_r[n_r] = ( s1 + ( ( s2 - s1 ) * ( zero_adc - adc[s]) ) / ( adc[s+1] - adc[s] ) );

         if( n_r > 0){
            p_r = x_r[ n_r ] - x_r[ n_r - 1 ];
            acc_period_sum += p_r;

            if( n_acc_period < 10000)acc_period_stock[ stck++ ] = p_r;

            n_acc_period++;
         }

         n_r++;
      }
   }

   return 0;
} //}}}


float WaveForm::GetPeriodMean( void ){ //{{{
   return acc_period_sum / (float)n_acc_period;
} //}}}


float WaveForm::GetPeriodError( void ){ //{{{

   float rms = 0.;
   float mean = GetPeriodMean();

   for( int i = 0;i < stck;i++){
      rms += pow( acc_period_stock[ i ] - mean , 2 ) / (float)stck;
   }

   return sqrt( rms / (float)n_acc_period );

} //}}}


int WaveForm::CalcPeriod( float baseline ){ //{{{

   float zero_adc = baseline;

   float x_r[1000] = { 0. } , x_f[1000] = { 0. };
   float xpoint[2][1000] = {{ -1 }}; //[0][] = falling edge , [1][] = rising edge
   int nx_fall = 0 , nx_rise = 0;

   //get cross point over baseline
   for( int s = 0;s < ( n_sample - 1 );s++){
      if( adc2f[s] > zero_adc && adc2f[s+1] < zero_adc ){ //falling edge 
         float s1 = time2f[ s ] , s2 = time2f[ s+1 ];

         xpoint[0][nx_fall++]  = ( s1 + ( ( s2 - s1 ) * ( zero_adc - adc2f[s]) ) / ( adc2f[s+1] - adc2f[s] ) );

      }else if( adc2f[s] < zero_adc && adc2f[s+1] > zero_adc ){ //rising edge
         float s1 = time2f[ s ] , s2 = time2f[ s + 1 ];

         xpoint[1][nx_rise++] = ( s1 + ( ( s2 - s1 ) * ( zero_adc - adc2f[s]) ) / ( adc2f[s+1] - adc2f[s] ) );

      }
   }


   //Fill the period
   int nperiod = 0;
   for(int i = 1;i < nx_fall;i++){
      float p = xpoint[0][i] - xpoint[0][i-1];

      if( period.size() <= nperiod){
         period.push_back( p );
         period_beggining_point.push_back( xpoint[0][i-1] );
         period_end_point.push_back( xpoint[0][i] );
      }else{
         period[nperiod] = p;
         period_beggining_point[nperiod] = xpoint[0][i-1];
         period_end_point[nperiod] = xpoint[0][i];
      }
      nperiod++;
   }

   for(int i = 1;i < nx_rise;i++){
      float p = xpoint[1][i] - xpoint[1][i-1];

      if( period.size() <= nperiod){
         period.push_back( p );
         period_beggining_point.push_back( xpoint[1][i-1] );
         period_end_point.push_back( xpoint[1][i] );
      }else{
         period[nperiod] = p;
         period_beggining_point[nperiod] = xpoint[1][i-1];
         period_end_point[nperiod] = xpoint[1][i];
      }
      nperiod++;
   }


   return nperiod;
} //}}}


float WaveForm::GetPeriodCrossPoint( int id , int begin_or_end ){ //{{{ Cross point for 0:begin 1:end

   if( begin_or_end < 1){
      return period_beggining_point[ id ];
   }else{
      return period_end_point[ id ];
   }

}//}}}


float WaveForm::GetPeriod( int id ){ //{{{
   return period[ id ];
} //}}}


//phase analysis
int WaveForm::SetPhaseADC128( int fwin , float base ){ //{{{
   //if( b_dbk )std::cout << std::endl << "*** Calib_dst23::SetPhaseADC ***" << std::endl;
   //define variable
   float zero_adc = base;

   //if( b_dbk )std::cout << "chk **** zero of adc : " << zero_adc << std::endl;


   //Get odd & even waveform 
   //if( b_dbk )std::cout << "chk *** even  : odd" << zero_adc << std::endl;
   float adc_odd[64] = { 0. } , adc_even[64] = { 0. };
   for( int s = 0;s < 64;s++){
      adc_even[ s ] = adc[2 * s];
      adc_odd[ s ] = adc[2 * s + 1];
   //if( b_dbk )std::cout << "chk *** " << adc_even[ s ] << " : " << adc_odd[s] << std::endl;
   }

   //checking unmatched start sign
   int start_flag = 1; //1 is matched
   if( ( ( adc_even[0] - zero_adc ) * ( adc_odd[0] - zero_adc ) ) < 0){
      start_flag = 0;
   }

   //scan 64 samples
   float xpoint_fall[2][64] = { 0. } , xpoint_rise[2][64] = { 0. }; //[ sample odd/even(0/1) ][ max number of Xpoints ]
   int i_nfall[2] = { 0 } , i_nrise[2] = { 0 }; //[ sample odd/even(0/1) ]
   //int i_lastxp_fall[2][2] = {{ -32 }} , i_lastxp_rise[2][2] = {{ -32 }}; //[ window odd/even(0/1) ][ sample odd/even(0/1) ]

   for( int s = 0;s < 64;s++){
      int win = ( s > 31 ) ? ( fwin + 1) % 2 : fwin ;

      if( start_flag ){
         //check cross point for even samples 
         if( adc_even[s] > zero_adc && adc_even[s+1] < zero_adc ){ //falling edge 

            if( s > 1 && s < 62){ //for skipping about break form

               if( ( adc_even[ s-2 ] >zero_adc && adc_even[ s-1 ] > zero_adc &&  adc_even[ s+2 ] < zero_adc && adc_even[ s+3 ] < zero_adc ) ){ 
                  float s1 = (float)s * 2 , s2 = (float)(s+1) * 2;
                  //xpoint_fall_even[i_nfall_even++] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_even[s]) ) / ( adc_even[s+1] - adc_even[s] );
                  xpoint_fall[0][ i_nfall[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_even[s]) ) / ( adc_even[s+1] - adc_even[s] );
                  //i_lastxp_fall[win][0] = s;

                  //std::cout << " chk *** fall even: xpoint=" << xpoint_fall_even[i_nfall_even-1] << " adc1=" << adc_even[s] << " , adc2=" << adc_even[s+1] <<  std::endl;
               }
            }else{
                  float s1 = (float)s * 2 , s2 = (float)(s+1) * 2;
                  xpoint_fall[0][ i_nfall[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_even[s]) ) / ( adc_even[s+1] - adc_even[s] );
                  //i_lastxp_fall_even = s;
            }

         }else if( adc_even[s] < zero_adc && adc_even[s+1] > zero_adc ){ //rising edge

            if( s > 1 && s < 30){ //for skipping about break form
               if( adc_even[s-2] < zero_adc && adc_even[s-1] < zero_adc &&  adc_even[s+2] > zero_adc && adc_even[s+3] > zero_adc ){ 

                  float s1 = (float)s * 2 , s2 = (float)(s+1) * 2;
                  xpoint_rise[0][ i_nrise[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_even[s]) ) / ( adc_even[s+1] - adc_even[s] );
                  //i_lastxp_rise_even = s;
                  //std::cout << " chk *** rise even: xpoint=" << xpoint_rise_even[i_nrise_even-1] << " adc1=" << adc_even[s] << " , adc2=" << adc_even[s+1] <<  std::endl;
               }
            }else{
                  float s1 = (float)s * 2 , s2 = (float)(s+1) * 2;
                  xpoint_rise[0][ i_nrise[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_even[s]) ) / ( adc_even[s+1] - adc_even[s] );
                  //i_lastxp_rise_even = s;
            }

         }

         //check cross point for odd samples 
         if( adc_odd[s] > zero_adc && adc_odd[s+1] < zero_adc ){ //falling edge 

            if( s > 1 && s < 30){ //for skipping about break form

               if( adc_odd[s-2] >zero_adc && adc_odd[s-1] > zero_adc &&  adc_odd[s+2] < zero_adc && adc_odd[s+3] < zero_adc ){ 
                  float s1 = (float)s * 2 + 1, s2 = (float)(s+1) * 2 + 1;
                  //xpoint_fall_odd[i_nfall_odd++] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_odd[s]) ) / ( adc_odd[s+1] - adc_odd[s] );
                  xpoint_fall[1][ i_nfall[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_odd[s]) ) / ( adc_odd[s+1] - adc_odd[s] );
                  //i_lastxp_fall_odd = s;

                  //std::cout << " chk *** fall odd: xpoint=" << xpoint_fall_odd[i_nfall_odd-1] << " adc1=" << adc_odd[s] << " , adc2=" << adc_odd[s+1] <<  std::endl;
               }
            }else{
                  float s1 = (float)s * 2 + 1, s2 = (float)(s+1) * 2 + 1;
                  xpoint_fall[1][ i_nfall[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_odd[s]) ) / ( adc_odd[s+1] - adc_odd[s] );
                  //i_lastxp_fall_odd = s;
            }

         }else if( adc_odd[s] < zero_adc && adc_odd[s+1] > zero_adc ){ //rising edge

            if( s > 1 && s < 30){ //for skipping about break form
               if( adc_odd[s-2] < zero_adc && adc_odd[s-1] < zero_adc &&  adc_odd[s+2] > zero_adc && adc_odd[s+3] > zero_adc ){ //for skipping about break form
                  float s1 = (float)s * 2 + 1, s2 = (float)(s+1) * 2 + 1;
                  //xpoint_rise_odd[i_nrise_odd++] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_odd[s]) ) / ( adc_odd[s+1] - adc_odd[s] );
                  xpoint_rise[1][ i_nrise[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_odd[s]) ) / ( adc_odd[s+1] - adc_odd[s] );
                  //i_lastxp_rise_odd = s;

                  //std::cout << " chk *** rise odd: xpoint=" << xpoint_rise_odd[i_nrise_odd-1] << " adc1=" << adc_odd[s] << " , adc2=" << adc_odd[s+1] <<  std::endl;
               }
            }else{
                  float s1 = (float)s * 2 + 1, s2 = (float)(s+1) * 2 + 1;
                  xpoint_rise[1][ i_nrise[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - adc_odd[s]) ) / ( adc_odd[s+1] - adc_odd[s] );
                  //i_lastxp_rise_odd = s;
            }
         }

      }else{

         if( ( ( adc_even[s] - zero_adc ) * ( adc_odd[s] - zero_adc ) ) > 0){
            start_flag = 1;
         }

      } //checking start flag
   } //sample scan

   //std::cout << std::endl;

   int evt[2] = { 0 };
   float sum[2] = { 0. };

   int fall_min = ( i_nfall[0] < i_nfall[1] ) ? i_nfall[0] : i_nfall[1] ;
   for( int j = 0;j < fall_min;j++){
      if( xpoint_fall[1][j] < 64 ){
         sum[0] += ( xpoint_fall[1][j] - xpoint_fall[0][j] );
         evt[0]++;
      }else{
         sum[1] += ( xpoint_fall[1][j] - xpoint_fall[0][j] );
         evt[1]++;
      }

      //std::cout << " chk *** time fall:" << ( xpoint_fall_even[j] - xpoint_fall_odd[j] ) << std::endl;
   }


   int rise_min = ( i_nrise[0] < i_nrise[1] ) ? i_nrise[0] : i_nrise[1] ;
   for( int j = 0;j < rise_min;j++){
      if( xpoint_rise[1][j] < 64 ){
         sum[0] += ( xpoint_rise[1][j] - xpoint_rise[0][j] );
         evt[0]++;
      }else{
         sum[1] += ( xpoint_rise[1][j] - xpoint_rise[0][j] );
         evt[1]++;
      }

      //std::cout << " chk *** time rise:" << ( xpoint_rise_even[j] - xpoint_rise_odd[j] ) << std::endl;
   }

   for( int win = 0;win < 2;win++){
      if( evt[win] > 0 ){
         fl_adcsum[ win ] += sum[ win ] / (float)evt[ win ];
         i_cnt[ win ]++;
      }
   }

   //getchar();

   return 0;
} //}}}


float WaveForm::GetPhaseDifference128( int id ){ //{{{ 
   //if( b_dbk )std::cout << std::endl << "*** Calib_dst23::GetPhaseDifference ***" << std::endl;

   //std::cout << "chk *** mean of dt:" <<  fl_adcsum[ id ] / (float)i_cnt[ id ] << " , id:" << id << std::endl;

   return  fl_adcsum[ id ] / (float)i_cnt[id];
} //}}}


int WaveForm::ClearPhaseAdc128( void ){ //{{{

   for( int i = 0;i < 2;i++){
      fl_adcsum[ i ] = 0.;
      i_cnt[ i ] = 0;
   }

   return 0;
} //}}}


void WaveForm::SetPhaseWaveform( float wf1[][2] , float wf2[][2]  ){ //{{{

   int n1 =  sizeof wf1 / sizeof wf1[0][2];
   int n2 =  sizeof wf2 / sizeof wf2[0][2];

   if( n1 != 0 && ( n1 != p_nsample || n2 != p_nsample ) ){
      std::cout << "<Warning in Waveform::SetPhaseWaveform,  nuber of sample is changed :" << n1 << std::endl; 
   }

   if( ( p_nsample < n1 ) || (p_nsample < n2) ){

      p_nsample = ( n1 > n2 ) ? n1 : n2;

      for( int i = p_nsample;i < n1;i++){
         p_time1.push_back( 0. );
         p_time2.push_back( 0. );
         p_adc1.push_back( 0. );
         p_adc2.push_back( 0. );
      }

      for( int i = 0;i < p_nsample;i++){
         p_adc1[i] = wf1[i][0];
         p_time1[i] = wf1[i][1];
         p_adc2[i] = wf2[i][0];
         p_time2[i] = wf2[i][1];
      }

   }

} //}}}


void WaveForm::SetPhaseWaveform( int even_odd , int id , float time , float adc  ){ //set 2 waveform to calc phase difference {{{

   for( int i = p_nsample;i <= id;i++){
      p_adc1.push_back( 0. );
      p_adc2.push_back( 0. );
      p_time1.push_back( 0. );
      p_time2.push_back( 0. );

      p_nsample++;
   }

   if( even_odd % 2 ){
      p_adc2[id] = adc;
      p_time2[id] = time;
   }else{
      p_adc1[id] = adc;
      p_time1[id] = time;
   }


} //}}}


float WaveForm::AccumPhaseDifference( float base ){ //{{{

   //define variable
   float zero_adc = base;
   const int i_num = p_nsample;

   //checking unmatched start sign
   int start_flag = 1; //1 is matched
   if( ( ( p_time1[0] - zero_adc ) * ( p_adc2[0] - zero_adc ) ) < 0){
      start_flag = 0;
   }


   //scan samples
   std::vector< std::vector<float> >  xpoint_fall( 2, std::vector<float>(p_nsample) );
   std::vector< std::vector<float> >  xpoint_rise( 2, std::vector<float>(p_nsample) );
   for( int i = 0;i < 2;i++){
      for( int j = 0;j < p_nsample;j++){
         xpoint_fall[i][j] = 0.;
         xpoint_rise[i][j] = 0.;
      }
   }

   int i_nfall[2] = { 0 } , i_nrise[2] = { 0 }; //[ sample odd/even(0/1) ]

   //for( int s = 0;s < 32;s++){
   for( int s = 0;s < 31;s++){ //ignore last point 

      if( start_flag ){
         //check cross point for even samples 
         if( p_adc1[s] > zero_adc && p_adc1[s+1] < zero_adc ){ //falling edge 

            if( s > 1 && s < 31){ //for skipping about break form

               if( ( p_adc1[ s-2 ] >zero_adc && p_adc1[ s-1 ] > zero_adc ) && (  p_adc1[ s+2 ] < zero_adc && p_adc1[ s+3 ] < zero_adc ) ){ 
                  float s1 = p_time1[ s ] , s2 = p_time1[ s+1 ];
                  xpoint_fall[0][ i_nfall[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc1[s]) ) / ( p_adc1[s+1] - p_adc1[s] );

               }
            }else{
                  float s1 = p_time1[ s ] , s2 = p_time1[ s+1 ];
                  xpoint_fall[0][ i_nfall[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc1[s]) ) / ( p_adc1[s+1] - p_adc1[s] );
            }

         }else if( p_adc1[s] < zero_adc && p_adc1[s+1] > zero_adc ){ //rising edge

            if( s > 1 && s < 30){ //for skipping about break form
               if( p_adc1[s-2] < zero_adc && p_adc1[s-1] < zero_adc &&  p_adc1[s+2] > zero_adc && p_adc1[s+3] > zero_adc ){ 

                  float s1 = p_time1[ s ] , s2 = p_time1[ s+1 ];
                  xpoint_rise[0][ i_nrise[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc1[s]) ) / ( p_adc1[s+1] - p_adc1[s] );
               }
            }else{
                  float s1 = p_time1[ s ] , s2 = p_time1[ s+1 ];
                  xpoint_rise[0][ i_nrise[0]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc1[s]) ) / ( p_adc1[s+1] - p_adc1[s] );
            }

         }

         //check cross point for odd samples 
         if( p_adc2[s] > zero_adc && p_adc2[s+1] < zero_adc ){ //falling edge 

            if( s > 1 && s < 30){ //for skipping about break form

               if( p_adc2[s-2] >zero_adc && p_adc2[s-1] > zero_adc &&  p_adc2[s+2] < zero_adc && p_adc2[s+3] < zero_adc ){ 
                  float s1 = p_time2[ s ] , s2 = p_time2[ s+1 ];
                  xpoint_fall[1][ i_nfall[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc2[s]) ) / ( p_adc2[s+1] - p_adc2[s] );
               }

            }else{
                  float s1 = p_time2[ s ] , s2 = p_time2[ s+1 ];
                  xpoint_fall[1][ i_nfall[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc2[s]) ) / ( p_adc2[s+1] - p_adc2[s] );
            }

         }else if( p_adc2[s] < zero_adc && p_adc2[s+1] > zero_adc ){ //rising edge

            if( s > 1 && s < 30){ //for skipping about break form
               if( p_adc2[s-2] < zero_adc && p_adc2[s-1] < zero_adc &&  p_adc2[s+2] > zero_adc && p_adc2[s+3] > zero_adc ){ //for skipping about break form
                  float s1 = p_time2[ s ] , s2 = p_time2[ s+1 ];
                  xpoint_rise[1][ i_nrise[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc2[s]) ) / ( p_adc2[s+1] - p_adc2[s] );
               }
            }else{
                  float s1 = p_time2[ s ] , s2 = p_time2[ s+1 ];
                  xpoint_rise[1][ i_nrise[1]++ ] = s1 + ( ( s2 - s1 ) * ( zero_adc - p_adc2[s]) ) / ( p_adc2[s+1] - p_adc2[s] );
            }
         }

      }else{

         if( ( ( p_adc1[s] - zero_adc ) * ( p_adc2[s] - zero_adc ) ) > 0){
            start_flag = 1;
         }

      } //checking start flag
   } //sample scan

   //std::cout << std::endl;

   int evt = 0;
   float sum = 0.;

   int fall_min = ( i_nfall[0] < i_nfall[1] ) ? i_nfall[0] : i_nfall[1] ;
   for( int j = 0;j < fall_min;j++){
      sum += ( xpoint_fall[1][j] - xpoint_fall[0][j] );
      evt++;
   }


   int rise_min = ( i_nrise[0] < i_nrise[1] ) ? i_nrise[0] : i_nrise[1] ;
   for( int j = 0;j < rise_min;j++){
      sum += ( xpoint_rise[1][j] - xpoint_rise[0][j] );
      evt++;
   }

   if( evt > 0 ){
      fl_p_adcsum += sum / (float)evt;
      i_p_cnt++;
   }

   return  sum / (float)evt;

} //}}}


float WaveForm::GetPhaseDifference( void ){ //{{{
   return  fl_p_adcsum / (float)i_p_cnt;
} //}}}


int WaveForm::DrawPhaseWaveform( void ){ //{{{

   TGraph *g_wf1 = new TGraph();
   TGraph *g_wf2 = new TGraph();
   for( int s = 0;s < p_nsample;s++){
      g_wf1->SetPoint( s , p_time1[s] , p_adc1[s] );
      g_wf2->SetPoint( s , p_time2[s] , p_adc2[s] );
      std::cout << "s:" << p_time1[s] << " , " << p_adc1[s] << std::endl;
   }

   gPad->SetGrid();
   g_wf1->SetMarkerStyle( wf_m_type );
   g_wf1->SetMarkerSize( wf_m_size );
   g_wf1->SetLineColor( 1 );
   g_wf1->SetTitle( "Comparison 2 waveform" );
   g_wf1->GetYaxis()->SetTitle( wf_ytitle );
   g_wf1->GetXaxis()->SetTitle( wf_xtitle );
   g_wf1->Draw( "apl" );

   g_wf2->SetMarkerStyle( wf_m_type );
   g_wf2->SetMarkerSize( wf_m_size );
   g_wf2->SetLineColor( 2 );
   g_wf2->Draw( "plsame" );

   return 1;
} //}}}


int WaveForm::ClearPhaseWaveform( void ){ //{{{

   for( int i = 0;i < 2;i++){
      fl_adcsum[ i ] = 0.;
      i_cnt[ i ] = 0;
   }

   p_adc1.clear();
   p_adc2.clear();
   p_time1.clear();
   p_time2.clear();

   p_nsample = 0;

   return 1;
} //}}}