//for wave form analysis
//made by arita 2012.1.4
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>
#include "TFile.h"
#include "TString.h"
#include "TTree.h"
#include "PeakSearch.h"


PeakSearch::PeakSearch(){ //{{{
   n_vect = 0;
   n_accum = 0;
   accum_time = 0;
}; //}}}


PeakSearch::~PeakSearch(){ //{{{
}; //}}}


void PeakSearch::GetWFfromSuperClass( void ){ //{{{

   adc.clear();
   time.clear();

   n_sample = GetWaveform( &time , &adc );

} //}}}


float PeakSearch::GetInnerPointLinear( float x1 , float x2 , float y1 , float y2 , float v_th ){ //{{{
   return  x1 + ( x2 - x1 ) * ( v_th - y1 ) / ( y2 - y1 ) ;
} //}}}


float PeakSearch::GetPeakId(int id ){ //{{{
   return peak_id[id];
} //}}}


float PeakSearch::GetPeakTime(int id ){ //{{{
   return peak_time[id];
} //}}}


float PeakSearch::GetPeakADC(int id ){ //{{{
   return peak_adc[id];
} //}}}


float PeakSearch::GetTimeMean( void ){ //{{{
   return accum_time / (float)n_accum;
} //}}}


int PeakSearch::LeadingEdge( float v_th , int edge ){ //{{{

   GetWFfromSuperClass();

   int j_edge = edge / abs( edge );
   bool vth_over = false;
   int n_peak = 0;

   for( int i = 1;i < n_sample;i++){

      //if( adc[i-1] < (j_edge * v_th ) && adc[i] > ( j_edge * v_th ) && !vth_over ){
      if( ( adc[i-1] * j_edge ) > 0 && ( adc[i] * j_edge ) > 0 ){
         if( fabs( adc[i-1] ) < fabs( v_th ) && fabs( adc[i] ) > fabs( v_th ) && !vth_over ){
            
            if( n_peak >= n_vect ){
               peak_adc.push_back( 0. );
               peak_id.push_back( 0. );
               peak_time.push_back( 0. );
               n_vect++;
            }

            peak_adc[n_peak] = v_th;
            peak_id[n_peak] = GetInnerPointLinear( (float)( i - 1) , (float)i , adc[i-1] , adc[i] , v_th );
            peak_time[n_peak] = GetInnerPointLinear( (float)( time[i-1] ) , (float)time[i] , adc[i-1] , adc[i] , v_th );
            
            n_accum++;
            accum_time += peak_time[n_peak];

            n_peak++;
            vth_over = true;
         }
      }

      //if( adc[i] < ( j_edge * v_th ) && vth_over )vth_over = false;
      if( fabs( adc[i] ) < fabs( v_th ) && vth_over )vth_over = false;
   }


   return n_peak;
} //}}}


int PeakSearch::PeakTop( float v_th , int gradiate_width ){ //{{{

   GetWFfromSuperClass();

   int n_peak = 0;

   if( gradiate_width < 0 ){
      std::cout << "<PeakTop> Error::gradiate_width should be plus value." << std::endl;
      return 0;
   }

   for( int i = gradiate_width;i < (n_sample - gradiate_width - 1);i++){
      if( fabs( adc[i] ) > fabs( v_th ) ){
         bool j_peak = true;
         //std::cout << "<PeakTop> peak search " << std::endl;

         for( int j = (i - gradiate_width);j <= (i + gradiate_width );j++){
            //std::cout << "<PeakTop> center:" << i << " , adc:" << adc[i] << " , com:" << j << " , adc:" << adc[j] << std::endl;
            if( fabs( adc[i] ) < fabs( adc[j] ) )j_peak = false;
         }

         if( j_peak ){
            if( n_peak >= n_vect ){
               peak_adc.push_back( 0. );
               peak_id.push_back( 0. );
               peak_time.push_back( 0. );
               n_vect++;
            }

            peak_id[n_peak] =  i ;
            peak_adc[n_peak] = adc[i];
            peak_time[n_peak] = time[i];

            n_peak++;
         }
      }
   }

   return n_peak;

} //}}}


int PeakSearch::WPeakTop( float v_th , int gradiate_width ){ //{{{
   GetWFfromSuperClass();

   if( gradiate_width < 0 ){
      std::cout << "<WPeakTop> Error::gradiate_width should be plus value." << std::endl;
      return 0;
   }

   //peak top search 
   const int n = PeakTop( v_th , gradiate_width );
   int n_peak = n;

   float peaktime_buf[n];

   //calc Weighted mean of time
   for( int p = 0;p < n_peak;p++){
      //std::cout << "<WPeakTop>peak :" << p << " , peak_id :" << peak_id[p] << " , adc :" << peak_adc[p] << " , time :" << peak_time[p] << std::endl;
      int peak = (int)peak_id[p]; int second = peak - 1; int third = peak + 1;

      peaktime_buf[p] = ( time[peak] * adc[peak]
                   +   time[second] * adc[second]
                   +   time[third] * adc[third] )
                   / ( adc[peak] + adc[second] + adc[third] );

   }

   for( int p = 0;p < n_peak;p++)peak_time[p] = peaktime_buf[p];

   return n_peak;

} //}}}