-
Chad Lantz authoredChad Lantz authored
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();
}