Skip to content
Snippets Groups Projects
WFAnalysis.cpp 21.36 KiB
/** @ingroup ana
 *  @file WFAnalysis.cpp
 *  @brief Implementation of WFAnalysis.
 *
 *  Function definitions for WFAnalysis are provided.
 *  This class is the main  class for the waveform analysis.
 *  It initializes the histograms used for output of processed events.
 *  Also includes methods that accept all waveforms in an event.
 *  This is where the analysis should be done.
 *
 *  @author Sheng Yang, Chad Lantz, Yakov Kulinich
 *  @bug No known bugs.
 */

#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TRandom.h>
#include <TTree.h>
#include <TCanvas.h>
#include <iostream>

#include "WFAnalysis.h"

/** @brief Default Constructor for WFAnalysis.
 */
WFAnalysis::WFAnalysis( ){

}


/** @brief Destructor for WFAnalysis.
 */
WFAnalysis::~WFAnalysis( ){

    delete f;
    delete nlFit;

}


/** @brief Initialization method for WFAnalysis
 *
 *  Can add other things here that you would
 *  perhaps not put into the constructor.
 *  I.e. a TTree, some tools. Etc.
 *
 */
void WFAnalysis::Initialize( ){

     nlFit = new TF1("nl1", "(x<660)*(-0.188891+1.03623*x) + \(x>660 && x<850)*( -51210.2 + 260.781*x - 0.4916*pow(x,2) + 0.00041227*pow(x,3) - (1.29681e-7)*pow(x,4) ) + \(x>850)*( -526857 + 1844.37*x - 2.1493*pow(x,2) + 0.000834971*pow(x,3) )", 0, 1000);

}


/** @brief Historgam Setup method for WFAnalysis
 *
 *  Should instantiate any histograms you wish to output here.
 *
 *  @return none
 */
void WFAnalysis::SetupHistograms( ){

  f = new TF1("f","gaus",-50,50);
  hDiffHisto = new TH1D("diffhist","First Derivative values histogrammed",75,-75,75);
  hPed = new TH1D("ped","ped", 1024, -512, 512);

}

/** @brief Analyze Events method for WFAnalysis
 *  @param vWF 1D vector of all waveforms for N ch
 *
 *  Here is the event-based analysis code.
 *  A const 1D vector of N TH1* is received.
 *  N is the number of channels
 *
 */
void WFAnalysis::AnalyzeEvent( const std::vector< TH1* >& vWFH ){

  // example... you can loop through all histos as follows
  for( unsigned int ch = 0; ch < vWFH.size(); ch++ ){
    TH1* h = vWFH[ch];
    for( unsigned int samp = 0; samp < h->GetNbinsX(); samp++ ){
      //Unused
    }
  }
}


/** @brief Analyze Events method for WFAnalysis
 *  @param vWF const 2D vector for all waveforms in event
 *
 *  Here is the event-based analysis code.
 *  A const 2D vector of NxM is received.
 *  N is the number of channels
 *  M is the number of samples per channel
 *
 *  @return none
 */
void WFAnalysis::AnalyzeEvent( const std::vector< std::vector< float > >& vWF ){

  // example... you can loop through each sample in each channel
  // Can do the same in previous AnalyzeEvent( .. ) function, just
  // looking at vWFH[ch]->GetBinContent( samp + 1 );
  for( unsigned int ch = 0; ch < vWF.size(); ch++ ){
    for( unsigned int samp = 0; samp < vWF[ch].size(); samp++ ){
      //Unused
    }
  }
}


/**
 * @brief Analyze Event method for WF analysis
 * @param vCh A vector of pointers to Channel objects
 *
 *  Receives a const vector of Channels. Loops over all Channels within
 *  the vector for analysis. Inverts channels with negative polarity before processing \n  \n
 *
 *  Processing occurs in the following order: \n
 *  Median filtering - If selected by the user, apply median filtering to the waveform \n
 *  FFT low pass filtering - If selected by the user, apply FFT low pass filter to the waveform \n
 *  Hit detection - Finds a hit within the waveform and gives the time if the hit was present \n
 *  Pedestal finding - Histograms the waveform outside of the hit window to determine pedestal and pedestal rms \n
 *  Pedestal subtraction - Removes the pedestal from the processed waveform \n
 *  Zero suppression - Sets all bins outside of the hit window and below a threshold to zero \n
 *  Integration - Integrates the processed waveform within the hit window accoutinging for bin timing calibration. Provides a value in pC \n
 *  Peak information - Uses the information gathered from the previous functions to save information about the waveform peak and derivative peak to the given channel
 *
 */
void WFAnalysis::AnalyzeEvent( const std::vector< Channel* > vCh ){
  for( Channel* Ch : vCh ){
    if( !Ch->is_on ) continue;

    //Clear old hits
    Ch->was_hit = false;
    Ch->nHits = 0;
    Ch->hit_window.first.clear();
    Ch->hit_window.second.clear();
    Ch->Peak_center.clear();
    Ch->Diff_Peak_center.clear();
    Ch->Peak_time.clear();
    Ch->Diff_Peak_time.clear();
    Ch->Peak_max.clear();
    Ch->Diff_max.clear();
    Ch->Charge.clear();
    Ch->nPhotons.clear();

    //retrieving information for each channel as a histogram
    TH1D* h = Ch->WF_histo;
    TH1D* hProcessed = Ch->PWF_histo;
    TH1D* hDiff = Ch->FirstDerivative;

    //Invert the signal if it has negative polarity
    if( !Ch->pos_polarity ){
      for(int bin = 0; bin < h->GetNbinsX(); bin++){
        h->SetBinContent( bin+1, -1*Ch->pWF->at(bin) );
        Ch->pWF->at(bin) = -1*Ch->pWF->at(bin) ;
      }
      h->ResetStats();
      h->GetYaxis()->SetRange(0, 0);
    }

    //Apply filters as requested by the user
    if( Ch->doMedianFilter ) MedianFilter( Ch );
    if( Ch->doLPF ) LowPassFilter( Ch, h );

    //If the hit window has been fixed by the user, set was_hit to true and set the
    //hit window to the user provided values
    //Otherwise, find the hit window
    if( m_FixHitWindow ){
      Ch->was_hit = true;
      Ch->nHits = 1;
      Ch->hit_window.first.push_back( m_fixed_hit_window.first );
      Ch->hit_window.second.push_back( m_fixed_hit_window.second );
    }else{
      //Get the first derivative and determine the hit window from it
      //if a valid hit window is not found was_hit is set to false
      GetDifferential( Ch );

      //We set the offset of the DRS4 channels to -250 at first. This lead to clipping
      //of the baseline. If that's the case, estimate the RMS to be 5 based on good runs.
      Ch->FirstDerivativeRMS = ( Ch->offset != -250 ) ? GetRMS( Ch ) : 5.0;
      FindHitWindow( Ch );
    }//end if else fixed hit window


    //Get and subtract the pedestal from the data outside the hit window
    GetPedestal( Ch );

    //If the channel was hit, proceed with processing
    if( Ch->was_hit ){
      SubtractPedestal( Ch );

      //The DRS4 saturates around 800-900 mV. If we have those values, flag it as saturated
      //std::cout << "MAX = " << Ch->PWF_histo->GetMaximum() << std::endl;
      //THIS IS ONLY FOR 2018 TEST BEAM - VERY RARE
      //if( Ch->PWF_histo->GetMaximum() > 890.0 ){
      //  Ch->saturated = true;
      //  std::cout << "====SATURATED=====" << std::endl;
      //}
      //Calibrate out the DRS4 non-linearity if requested by the user
      if( Ch->DRS4NLcomp ) DRS4Cal( Ch );
      //Zero Suppress the processed waveform vector and retrieve the energy related values from it
      ZeroSuppress( Ch );

      for(int hit = 0; hit < Ch->nHits; hit++){
        Ch->Charge.push_back( GetCharge( Ch, Ch->hit_window.first[hit], Ch->hit_window.second[hit] ) );
        Ch->nPhotons.push_back( -Ch->Charge.back()*6.241e6/Ch->gain ); //6.241e6 is electrons per pC

        //Get peak information from the waveform using the hit window
        Ch->Peak_center.push_back( GetMaximumBin( hProcessed, Ch->hit_window.first[hit], Ch->hit_window.second[hit] ) );
        Ch->Peak_max.push_back( Ch->adc_mV*hProcessed->GetBinContent( Ch->Peak_center.back() ) );
        Ch->Peak_time.push_back( Ch->pTimeVec->at( Ch->Peak_center.back() ) );

        //Set the range of the fit function to just around the peak and feed it some starting parameters
        //Then perform the fit
        // int range = 0.05*(Ch->hit_window.second[hit] - Ch->hit_window.first[hit]);
        // Ch->FitFunc->SetRange( Ch->Peak_center.back() - range, Ch->Peak_center.back() + range);
        // Ch->FitFunc->SetParameters(Ch->Peak_max, 3, 0.5, Ch->hit_window.first[hit]);
        // hProcessed->Fit(Ch->FitFunc,"QNR");//Fit (Q)uietly, do (N)ot draw and use the specified (R)ange
      }//end loop over hits
    }//end if channel->was_hit
    hProcessed->ResetStats();
    hProcessed->GetYaxis()->SetRangeUser(0, 0);
    hDiff->ResetStats();
    hDiff->GetYaxis()->SetRange(0, 0);
  }//end channel loop
}//end AnalyzeEvent


/** @brief MedianFilter method for WFAnalysis
 *  @param Ch channel to perform Median filtering on
 *
 *  Uses the median of the raw waveform in a sliding
 *  window to determine the value of the bin at the
 *  center of the window
 *
 */
void WFAnalysis::MedianFilter( Channel* Ch ){
  int nBins = Ch->WF.size();
  std::vector< float > output( nBins, 0 );
  std::vector< float > window;

  //Fill the front and back of the output vector since we can't filter
  //those bins
  for(int bin = 0; bin < Ch->medFiltHalfWindow; bin++){
    output[bin] = Ch->WF[bin];
    output[nBins -1 - bin] = Ch->WF[nBins -1 - bin];
  }
  //Go through each bin except for the first and last halfWindow bins
  //Set the filtered bin value to the median of values of the bins in
  //the window around the current bin
  for(int bin = Ch->medFiltHalfWindow; bin < nBins - Ch->medFiltHalfWindow; bin++){
    for(int i = 0; i < 2*Ch->medFiltHalfWindow; i++){
      window.push_back( Ch->WF[ bin - Ch->medFiltHalfWindow + i ] );
    }
    output[bin] = TMath::Median(2*Ch->medFiltHalfWindow, &window[0]);
    window.clear();
  }
  //Replace the raw waveform with the filtered one
  for(int bin = 0; bin < nBins; bin++){
    Ch->WF[bin] = output[bin];
  }
}


/** @brief GetDifferential method for WFAnalysis
 *  @param Ch channel who's waveform will be differentiated
 *
 * The derivative is calculated as the difference of the summation N points before and after the ith point
 *
 *  @f[
 *       \delta_i(N)= \sum_{k=1}^{N} (s_i + k)  -  \sum_{k=1}^{N} (s_i - k)
 *  @f]
 * N is given by Ch::Threshold and is set by Detector::SetWFAthresholds()
 */
void WFAnalysis::GetDifferential( Channel* Ch ){


  //Pad the front of the vector with 0s
  for(int bin = 0; bin < Ch->diffSmoothing; bin ++) Ch->vDiff.push_back( 0.0 );


  //decalare temp variable to store the sum before and after the i data point
  double sum_before = 0;
  double sum_after = 0;

  // Generate 2 sliding window, sum_before and sum_after
  // Calculate the sum before and after the start data point
  for (int i = 1; i <= Ch->diffSmoothing; i++  ){
    sum_before += Ch->WF.at( Ch->diffSmoothing - i );
    sum_after  += Ch->WF.at( Ch->diffSmoothing + i );
  }

  // Loop over histogram
  for( unsigned int bin = Ch->diffSmoothing; bin < Ch->WF_histo->GetNbinsX() - Ch->diffSmoothing - 1 ; bin++ ){
        //set the bin to the calculated derivative value
        Ch->vDiff.push_back( sum_after - sum_before);
        Ch->FirstDerivative->SetBinContent(bin,(sum_after - sum_before));

        // Move the sum_before window forward
        sum_before += Ch->WF.at( bin );
        sum_before -= Ch->WF.at( bin - Ch->diffSmoothing );

        // Move the sum_after window forward
        sum_after -= Ch->WF.at( bin + 1);
        sum_after += Ch->WF.at( bin + Ch->diffSmoothing + 1);

    }//end derivative loop

    //Pad the back of the vector with 0s
    for(int bin = 0; bin <= Ch->diffSmoothing; bin ++) Ch->vDiff.push_back( 0.0 );
}

/** @brief GetRMS method for WFAnalysis
 *  @param Input Channel
 *  @return Width of gaussian fit
 *
 *  Histgrams the second derivative for a given Channel and gets the baseline RMS
 *  from a gaussian fit which ignores tails created by the peaks.
 *
 */
double WFAnalysis::GetRMS( Channel* Ch ){

    Double_t xmin,xmax;
    Ch->FirstDerivative->GetMinimumAndMaximum(xmin,xmax);
    if(xmax == 0) return 0;

    hDiffHisto->Reset();
    hDiffHisto->ResetStats();
    hDiffHisto->SetBins( (xmax - xmin)/36, xmin/3, xmax/3);

    //Loop over the histogram excluding the window used for differentiating to fill hRMS
    Int_t nbins = Ch->FirstDerivative->GetNbinsX();
    for(int bin = Ch->diffSmoothing; bin < nbins - Ch->diffSmoothing; bin++){
        hDiffHisto->Fill( Ch->vDiff.at( bin ) );
    }

    //Set the range of the gaussian fit to something reasonable, then perform the fit
    f->SetRange( hDiffHisto->GetMean() - 1.5*hDiffHisto->GetRMS(), hDiffHisto->GetMean() + 1.5*hDiffHisto->GetRMS() );
    f->SetParameters( hDiffHisto->GetMaximum() , hDiffHisto->GetMean(), hDiffHisto->GetRMS()/2.0);
    hDiffHisto->Fit("f","qR");

    //Return parameter 2. "gaus" is [0]*exp(-0.5*((x-[1])/[2])**2)
    //And reset the histogram
    return f->GetParameter(2);
}


/** @brief Defines the hit window for a given channel
  * @param ch Channel to be processed
  *
  * This function first checks if the FirstDerivative has any bins above threshold. If not, there will be no hit recorded.
  * In this case, we set datum to 0.0 in the case of float/doubles, and -1 in the case of ints.
  * If the derivative has bins above threshold the function steps through the waveform first derivative bin by bin looking for the first bin above threshold.
  * Once it finds the first bin above threshold it looks for the following landmarks in order:
  * FirstDerivative maximum, FirstDerivative zero crossing, FirstDerivative minimum, FirstDerivative zero crossing.
  * Which give us the corresponding datum:
  * Hit window start, Peak center, N/A, Hit window end.
  */
void WFAnalysis::FindHitWindow( Channel* ch ){

  double threshold = ch->Threshold*ch->FirstDerivativeRMS;

  //If there are no bins above threshold, we won't find a hit
  if( ch->FirstDerivative->GetMaximum() < threshold ||  ch->FirstDerivative->GetMinimum() > -1*threshold ) return;

  hPed->Reset();
  int min = ch->WF_histo->GetMinimum();
  int max = ch->WF_histo->GetMaximum();
  hPed->SetBins( (max - min)/2, min, max); //Set the range of the histogram tailored to each waveform


  int nBins = ch->WF.size();
  for(int bin = 0; bin < nBins; bin++ ){
    hPed->Fill( ch->WF[bin] );
  }

  double mean = hPed->GetMean();
  double rms = hPed->GetRMS();

  hPed->Reset();


  float currentVal, prevVal;
  int startBin=0, maxBin=0, zeroBin=0, minBin=0, endBin=0;
  //Search for the first bin above threshold
  for(int bin = 1; bin < nBins; bin++){
    startLooking:
    startBin = maxBin = zeroBin = minBin = endBin = 0;
    if( ch->vDiff.at(bin) < threshold )continue;

    startBin = bin;

    //From the start of the hit window, look for the maximum of the derivative
    prevVal = ch->vDiff.at(bin);
    for( bin; bin < nBins; bin++){
      currentVal = ch->vDiff.at(bin);

      //If we have passed the maximum, the previous value was the max
      if( currentVal < prevVal ){
        maxBin = bin - 1;
        max = prevVal;
        break;
      }else{
        prevVal = currentVal;
      }//end if else currentVal < prevVal
    }//end looking for maxBin

    //From the maxBin, start looking for the zero crossing
    for(bin; bin < nBins; bin++){
      // if( ch->FirstDerivative->GetBinContent(bin) < 0.0 && ch->WF.at(bin) > 1.5*rms ){
      if( ch->vDiff.at(bin) < 0.0 ){
        zeroBin = bin - 1;
        break;
      }
    }//end looking for zero crossing

    //From the zero bin, start looking for the minimum bin
    //If the minimum bin is below threshold (technically above) go back to the beginning
    prevVal = ch->vDiff.at( zeroBin );
    for(bin; bin < nBins; bin++){
      currentVal = ch->vDiff.at( bin );
      if( currentVal > prevVal ){

        // if( prevVal > -1*threshold )continue;
        if( prevVal > -1*threshold )goto startLooking;

        minBin = bin - 1;
        min = prevVal;
        break;
      }else{
        prevVal = currentVal;
      }//end if currentVal > prevVal
    }//end looking for minimum

    if(minBin == 0){
      startBin = maxBin = zeroBin = 0;
      continue;
    }

    //Find the next time the derivative goes to zero,
    //this is the end of the hit
    for(bin; bin < nBins; bin++){
      if( ch->vDiff.at(bin) > 0.0 ){
        endBin = bin - 1;
        break;
      }
    }

    //If any of the values remain at default, start the loop over
    //This also exits the loop if we made it to the last bin
    if( startBin == 0 || maxBin == 0 || zeroBin == 0 || minBin == 0 || endBin == 0 ) continue;

    //If we've gotten here, we should have a hit
    //For now just assign values and break
    //TODO: Add multiple hit storage
    ch->was_hit = true;
    ch->nHits++;
    ch->hit_window.first.push_back( startBin );
    ch->hit_window.second.push_back( endBin );
    ch->Diff_Peak_center.push_back( maxBin );
    ch->Diff_Peak_time.push_back( ch->pTimeVec->at(maxBin) );
    ch->Diff_max.push_back( ch->vDiff.at(maxBin) );

  }//end loop over all bins

}//end FindHitWindowNew


/** @brief Finds the pedestal and pedestal rms of the raw waveform
 *  @param ch Channel to be processed
 *
 *  Histograms data from the raw waveform excluding the hit window.
 *  Saves the mean and rms of that histogram as PedMean and PedRMS.
 */
void WFAnalysis::GetPedestal( Channel* ch ){
    int nBins = ch->WF_histo->GetNbinsX();
    int min = ch->WF_histo->GetMinimum();
    int max = ch->WF_histo->GetMaximum();
    hPed->SetBins( (max - min)/2, min, max); //Set the range of the histogram tailored to each waveform

    //If the channel was hit, we only calculate the pedestal until the first hit
    //Otherwise use every bin
    int end = ( ch->was_hit ) ? ch->hit_window.first.front() : nBins;
    for(int bin = 0; bin < end; bin++){
          hPed->Fill( ch->pWF->at(bin) );
    }

    ch->PedMean = hPed->GetMean();
    ch->PedRMS  = hPed->GetRMS();

}



/** @brief Subtracts channel pedestal from the processed waveform histogram
*  @param ch Channel to be processed
*
*  Saves a pedestal subtracted version of the waveform to PWF_histo
*/
void WFAnalysis::SubtractPedestal( Channel* ch ){

  int nBins = ch->WF.size();
  for(int bin = 0; bin < nBins; bin++){
    ch->PWF_histo->SetBinContent( bin, ch->pWF->at(bin) - ch->PedMean );
  }

}


/** @brief Supress values below PedRMS
 *
 *  Changes all values outside of the hit window smaller than PedRMS after
 *  pedestal subtraction to zero.
 */
void WFAnalysis::ZeroSuppress( Channel* ch ){
    int nBins = ch->PWF_histo->GetNbinsX();
    for(int bin = 0; bin < nBins; bin++){
        double content = ch->PWF_histo->GetBinContent(bin);

        //Skip the hit windows
        for(int hit = 0; hit < ch->nHits; hit++){
          if( bin == ch->hit_window.first[hit] ){
            bin = ch->hit_window.second[hit] + 1;
          }
        }

        if( content <= 1.5*ch->PedRMS ){ ch->PWF_histo->SetBinContent( bin, 0 ); }
    }
}

/** @brief Applies a low pass filter to the processed waveform
 *  @param ch Channel to be processed
 *
 *  Applies FFT to the raw waveform, removes high frequencies
 *  and reconstructs the signal in PWF_histo.
 */
void WFAnalysis::LowPassFilter( Channel* ch, TH1D* hIn ){
   if(!hIn){ hIn = ch->WF_histo; }
   Int_t n = hIn->GetNbinsX();
   double re[n], im[n], reOut[n], imOut[n];

   TH1 *hm =0;
   hm = hIn->FFT(hm, "MAG");

   // Apply the filter
   TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
   fft->GetPointsComplex(re,im);
   for(int i = 0; i< ch->LPFfreq; i++){
       reOut[i] = re[i];
       imOut[i] = im[i];
   }
   for(int i = ch->LPFfreq; i < n; i++){
       reOut[i] = 0.0;
       imOut[i] = 0.0;
   }

   //Make a new TVirtualFFT with an inverse transform. Set the filtered points and apply the transform
   TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &n, "C2R M K");
   fft_back->SetPointsComplex(reOut,imOut);
   fft_back->Transform();
   ch->PWF_histo = (TH1D*)TH1::TransformHisto(fft_back,ch->PWF_histo,"MAG");

   //The transform scales the signal by sqrt(n) each time, so rescale by 1/n
   ch->PWF_histo->Scale(1.0/n);

   delete hm;
   delete fft;
   delete fft_back;

}


/** @brief Provides calibration for DRS4 voltage response
 *  @param ch Channel to be processed
 *
 *  Uses a TF1 with the DRS4 response to convert from the recorded value in mV
 *  to the actual value of the waveform in mV and reconstructs the signal in PWF_histo.
 *s
 */
void WFAnalysis::DRS4Cal( Channel* ch ){

    for(int bin = 0; bin < ch->PWF_histo->GetNbinsX(); bin++){
        ch->PWF_histo->SetBinContent(bin, nlFit->Eval( ch->PWF_histo->GetBinContent(bin) ) );
    }
}

/** @brief Uses time vector and input resistance to find the actual charge detected by the DRS4
 *  @param ch Channel to be processed
 *
 *  Integrates the processed waveform over the hit window. Uses calibrated timing information
 *  from the DRS4 to determine charge in pC
 */
double WFAnalysis::GetCharge( Channel* ch, int start, int end){
  double charge = 0.0;
  for(int bin = start; bin < end; bin++){
    //Charge for a given time bin = dt*I = (t(bin+1)-t(bin))*V/R
    //Units are Charge(pC), time(ns), Voltage(mV), resistance(ohm), current(Amp = C/s)
    charge += (ch->pTimeVec->at(bin+1) - ch->pTimeVec->at(bin)) * ch->adc_mV * ch->PWF_histo->GetBinContent(bin)/Rin;
  }//end hit window
  return charge;
}//end GetCharge


/** @brief Get the maximum bin in a given range
 *  @return maxBin
 *  @param h - histogram under inspection
 *  @param _start - First bin of the desired range
 *  @param _end - Last bin in the desired range
 *
 *  Finds the maximum bin for the given histogram in the given range (in bins)
 *
 */
int WFAnalysis::GetMaximumBin( TH1* h, int _start, int _end ){
  float max = h->GetBinContent(_start);
  int maxBin = _start;

  for(int bin = _start+1; bin < _end; bin++){
    //If this bin isn't larger than the maximum identified, continue to the next
    if( h->GetBinContent(bin) < max) continue;

    max = h->GetBinContent(bin);
    maxBin = bin;
  }
  return maxBin;
}

/** @brief Finalize method for WFAnalysis
 *  @return none
 *
 *  Write histograms, TTree if it exists.
 *
 */
void WFAnalysis::Finalize( ){

  // If these exist...
  // m_tree->Write();
  // m_hist->Write();
}