PulseTemplate Class Reference

#include <PulseTemplate.h>

Public Member Functions
	PulseTemplate (int nSigma=3, int refine=100, double convergence=0.01)
	~PulseTemplate ()
void	AddPulse (TPulseIsland *pulse, double shift)
double	Correlation (TPulseIsland *pulse, double &ped, double &amp, double &time)
double	Chi2Fit (TPulseIsland *pulse, double &pedestal, double &amplitude, double &time)
bool	Converged ()
void	Save (TFile *outfile)
void	Print ()
int	GetNumberOfPulses ()
Private Attributes
TH1D *	fTemplate
TFitterMinuit *	fFitter
int	fNPulses
int	fNSigma
int	fRefine
int	fNBins
double	fClockTick
double	fSumX
double	fSumX2
double	fSum2X
double	fConvergence
double	fConvergenceLimit

Detailed Description

Definition at line 11 of file PulseTemplate.h.

Constructor & Destructor Documentation

PulseTemplate::PulseTemplate	(	int	nSigma = `3`,
		int	refine = `100`,
		double	convergence = `0.01`
	)

Definition at line 15 of file PulseTemplate.cpp.

References fFitter.

00015                                                                        : fTemplate(NULL), fFitter(NULL), fNPulses(0),
00016                                                                            fNSigma(nSigma), fRefine(refine), fNBins(0),
00017                                                                            fClockTick(0.), fSumX(0.), fSumX2(0.),
00018                                                                            fSum2X(0.), fConvergence(0.), fConvergenceLimit(convergence) {
00019   HistogramFitFCN* fcn = new HistogramFitFCN();
00020   fFitter = new TFitterMinuit(3); //  Three (3) parameters to modify (amplitude, time, pedestal)
00021   fFitter->SetMinuitFCN(fcn);
00022 }

PulseTemplate::~PulseTemplate ( )

Definition at line 25 of file PulseTemplate.cpp.

References fFitter, and fTemplate.

00025                               {
00026   delete fTemplate;
00027   delete fFitter; // Fitter should own HistogramFitFCN and destroy it according to documentation
00028 }

Member Function Documentation

void PulseTemplate::AddPulse	(	TPulseIsland *	pulse,
		double	shift
	)

Definition at line 30 of file PulseTemplate.cpp.

References fClockTick, fConvergence, fNBins, fNPulses, fNSigma, fRefine, fSum2X, fSumX, fSumX2, fTemplate, TPulseIsland::GetBankName(), TPulseIsland::GetClockTickInNs(), TPulseIsland::GetPedestal(), TPulseIsland::GetSamples(), and TPulseIsland::GetTriggerPolarity().

Referenced by MakeTemplate::ProcessEntry().

00030                                                               {
00031   // Function to average in pulse with template
00032   // Input--
00033   // pulse:     Pulse to average in
00034   // shift:     Bin shift (timing offset of peak)
00035   TH1D* old_template;
00036   double norm;
00037   double peak;
00038   double sigma;
00039   double pol;
00040   double ped;
00041   int nsamps;
00042   std::vector<int> samples;
00043   std::vector<double> rectified_samples, reshaped_pulse;
00044 
00045   if (fNPulses > 0)
00046     old_template = (TH1D*)fTemplate->Clone("old_histogram");
00047   pol = (double)pulse->GetTriggerPolarity();
00048   ped = pulse->GetPedestal(0);
00049   samples = pulse->GetSamples();
00050   ped = (double)(samples[0]+samples[1]+samples[2]+samples[3]) / 4.; /*** TEMPERARY PEDESTAL ***/
00051   nsamps = samples.size();
00052   for (int i = 0; i < nsamps; ++i)
00053     rectified_samples.push_back(pol*((double)samples[i]-ped));
00054   // Get peak value for normalization
00055   peak = 0.;
00056   norm = pol * (double)samples[0];
00057   for (int i = 1; i < nsamps; ++i) {
00058     if ((double)samples[i] * pol > norm) {
00059       peak = (double)i;
00060       norm = (double)samples[i];
00061     }
00062   }
00063   norm *= pol;
00064   // Get sigma for scaling
00065   sigma = 0.;
00066   for (int i = 0; i < nsamps; ++i) {
00067     std::cout << i << "\t" << sigma << "," << rectified_samples[i] << "\t";
00068     if (i % 10 == 0)
00069       std::cout << std::endl;
00070     sigma += rectified_samples[i] * std::pow(i - peak, 2.);
00071   }
00072   sigma = std::sqrt(std::abs(sigma)) / (double)nsamps;
00073 
00074   // If this is the first pulse, setup some variables
00075   if (fNPulses == 0) {
00076     fClockTick = pulse->GetClockTickInNs() / (double)fRefine;
00077     fNBins = (int)(2. * (double)fNSigma * sigma) * fRefine;
00078     std::cout << "Making histogram: " << fNBins << " " << fNSigma << " " << sigma << " " << nsamps << " " << ped << std::endl;
00079     TString str = "template_";
00080     str += pulse->GetBankName();
00081     fTemplate = new TH1D(str, "Template", fNBins, -(double)fNSigma, (double)fNSigma);
00082   }
00083 
00084   // Reshape pulse
00085   // To go from course binning somehwere in the pulse
00086   // to fine binning in the template, just use
00087   // a linear f(x)=mx+b
00088   int pulse_index;
00089   double m = 2. * (double)fNSigma * sigma / (double)fNBins;
00090   double b = (peak + shift) - (double)fNSigma * sigma;
00091   for (int i = 0; i < fNBins; ++i) {
00092     pulse_index = (int)((double)i * m + b + 0.5); // The half shift is for bin-centering since int floors the double
00093     reshaped_pulse.push_back((double)samples[pulse_index] / norm);
00094   }
00095 
00096   // If this is the first pulse added to the template, smooth it.
00097   // In C++11, we'll have direct access to the reshaped pulse's
00098   // data array with vector.data(), but now we have to use
00099   // a somewhat circuitous route. TH1::SmoothArray takes an array
00100   // of doubles to smooth as an argument.
00101   if (fNPulses == 0) {
00102     unsigned int l = reshaped_pulse.size();
00103     double *reshaped_pulse_array = new double[l];
00104     for (unsigned int i = 0; i < l; ++i)
00105       reshaped_pulse_array[i] = reshaped_pulse[i];
00106     TH1::SmoothArray(reshaped_pulse.size(), reshaped_pulse_array, 1);
00107     for (unsigned int i = 0; i < l; ++i)
00108       reshaped_pulse.at(i) = reshaped_pulse_array[i];
00109     delete [] reshaped_pulse_array;
00110   }
00111 
00112   // Add in pulse to template
00113   // Histogram indexing is off by one (index 0 is the underflow bin)
00114   // Hence the (i+1), which really corresponds to the ith element
00115   // in the reshaped_pulse array.
00116   double x, x_old;
00117   double e, e_old;
00118   for (int i = 0; i < fNBins; ++i) {
00119     x_old = fTemplate->GetBinContent(i+1);
00120     e_old = fTemplate->GetBinError(i+1);
00121     x = (double)fNPulses * x_old + reshaped_pulse[i];
00122     x /= (double)(fNPulses + 1);
00123     e = (double)(fNPulses - 1) * std::pow(e_old, 2.) + ((double)reshaped_pulse[i] - x_old) * ((double)reshaped_pulse[i] - x);
00124     e = std::sqrt(e / fNPulses);
00125     fTemplate->SetBinContent(i+1, x);
00126     fTemplate->SetBinError(i+1, e);
00127   }
00128 
00129   // Find the change and average it with previous changes.
00130   // That is the dot product of the last template with
00131   // the new one subtracted from 1.
00132   if (fNPulses > 0) {
00133     double mag_old = old_template->Integral(1, fNBins);
00134     double mag_new = fTemplate->Integral(1,fNBins);
00135     double dot_product = 0.;
00136     for (int i = 1; i <= fNBins; ++i)
00137       dot_product += old_template->GetBinContent(i) * fTemplate->GetBinContent(i);
00138     double cos_similarity = dot_product / (mag_old * mag_new);
00139     fConvergence = ((double)(fNPulses - 1) * fConvergence + (1 - cos_similarity))/(double)fNPulses;
00140   }
00141 
00142   // Recalculate values for correlation coefficient
00143   fSumX = 0.;
00144   fSumX2 = 0.;
00145   int c = 0;
00146   for (int i = 1; i <= fNBins; ++i) {
00147     c = fTemplate->GetBinContent(i);
00148     fSumX += c;
00149     fSumX2 += std::pow(c, 2.);
00150   }
00151   fSum2X = std::pow(fSumX, 2.);
00152 
00153   // Keep track of the number of pulses
00154   ++fNPulses;
00155 }

double PulseTemplate::Chi2Fit	(	TPulseIsland *	pulse,
		double &	pedestal,
		double &	amplitude,
		double &	time
	)

Definition at line 167 of file PulseTemplate.cpp.

References fFitter, fRefine, fTemplate, and TPulseIsland::GetSamples().

Referenced by MakeTemplate::ProcessEntry().

00167                                                                                          {
00168   // First make a histogram out of the pulse with the same bin width as the template,
00169   // Then pass to fitter.
00170   // Returns ped in units of ADC counts, amp in... scale units(?), and time in units
00171   // of bins since the beginning of pulse
00172   std::vector<int> samples = pulse->GetSamples();
00173   int nel = samples.size();
00174   TH1D* rebinned_pulse = new TH1D("pulse2fit", "pulse2fit", nel * fRefine, -0.5, nel - 0.5);
00175   for (int i = 0; i < nel; ++i)
00176     for (int j = 1; j <= fRefine; ++j)
00177       rebinned_pulse->SetBinContent(i * fRefine + j, samples.at(i));
00178   time *= fRefine;
00179 
00180   // Prepare for minimizations
00181   fFitter->Clear();
00182   HistogramFitFCN* fcn = (HistogramFitFCN*)fFitter->GetMinuitFCN();
00183   fcn->SetH1(fTemplate);
00184   fcn->SetH2(rebinned_pulse);
00185   fFitter->SetParameter(0, "Pedestal", ped, 0.1, 0, 0);
00186   fFitter->SetParameter(1, "Amplitude", amp, 0.1, 0, 0);
00187   fFitter->SetParameter(2, "Time", time, 1., 0, 0); // Timing should have step size no smaller than binning,
00188                                                     // *IF* the fourth argument is step size this is okay,
00189                                                     // or later implement some interpolation method, note
00190                                                     // *DERIVATIVES* at bounderies of interpolation may cause
00191                                                     // problems since MIGRAD (the default method) relies on
00192                                                     // these heavily.
00193   fFitter->CreateMinimizer();
00194 
00195   // Minimize and notify if there was a problem
00196   int status = fFitter->Minimize();
00197   if (status != 0)
00198     std::cout << "ERROR: Problem with fit (" << status << ")!" << std::endl;
00199 
00200   // Get parameters and return chi2 of fit
00201   ped  = fFitter->GetParameter(0);
00202   amp = fFitter->GetParameter(1);
00203   time = fFitter->GetParameter(2) / (double)fRefine;
00204   fcn = (HistogramFitFCN*)fFitter->GetMinuitFCN();
00205   std::vector<double> params;
00206   params.push_back(ped);
00207   params.push_back(amp);
00208   params.push_back(time);
00209 
00210   return (*fcn)(params);
00211 }

bool PulseTemplate::Converged ( )

Definition at line 213 of file PulseTemplate.cpp.

Referenced by MakeTemplate::ProcessEntry().

00213                               {
00214   return (fConvergence <= fConvergenceLimit && fNPulses >= 2);
00215 }

double PulseTemplate::Correlation	(	TPulseIsland *	pulse,
		double &	ped,
		double &	amp,
		double &	time
	)

Definition at line 157 of file PulseTemplate.cpp.

00157                                                                                              {
00158   // Fits by finding the maximum correlation coefficient of template with pulse
00159   // Input--
00160   // pulse:    Pulse to fit
00161   // Output--
00162   // ped:      Pedestal of fit
00163   // amp:      Amplitude of fit
00164   return 0.;
00165 }

int PulseTemplate::GetNumberOfPulses ( )

Definition at line 245 of file PulseTemplate.cpp.

References fNPulses.

Referenced by MakeTemplate::ProcessEntry().

00245                                      {
00246   return fNPulses;
00247 }

void PulseTemplate::Print ( )

Definition at line 228 of file PulseTemplate.cpp.

References fClockTick, fConvergence, fConvergenceLimit, fNBins, fNPulses, fNSigma, fRefine, fSum2X, fSumX, and fSumX2.

Referenced by MakeTemplate::ProcessEntry().

00228                           {
00229   std::cout << "Pulse Template Info:"
00230             << "\t" << "NP" << "\t" << "NS"
00231             << "\t" << "Ref" << "\t" << "NB"
00232             << "\t" << "T" << "\t" << "SX"
00233             << "\t" << "SX2" << "\t" << "S2X"
00234             << "\t" << "Conv" << "\t" << "ConLim"
00235             << std::endl;
00236   std::cout << "                    "
00237             << "\t" << fNPulses << "\t" << fNSigma
00238             << "\t" << fRefine << "\t" << fNBins
00239             << "\t" << fClockTick << "\t" << fSumX
00240             << "\t" << fSumX2 << "\t" << fSum2X
00241             << "\t" << fConvergence << "\t" << fConvergenceLimit
00242             << std::endl;
00243 }

void PulseTemplate::Save ( TFile * outfile )

Definition at line 217 of file PulseTemplate.cpp.

References fTemplate.

Referenced by MakeTemplate::ProcessEntry().

00217                                      {
00218   TDirectory* cwd = gDirectory;
00219   if (ofile == NULL || !ofile->IsOpen()) {
00220     std::cout << "ERROR: Cannot save template to file (file NULL or not open)!" << std::endl;
00221     return;
00222   }
00223   ofile->cd();
00224   fTemplate->Write();
00225   cwd->cd();
00226 }