Hough transform for QEvents

hough.cpp

#include <cmath>
#include <TH2F.h>
#include <TVector3.h>
#include <QEvent.h>
#include <QGlobals.h>
#include <QPMTxyz.h>

TH2F* hough(const QEvent* event, TVector3& vfit, float tres_cut, bool qweight) {
  // Calculate Cherenkov angle using FTK energy
  float beta = sqrt(1.0 - (0.511 / (0.511 + e)) * (0.511 / (0.511 + e)));
  float tcher = acos(1.0 / (kWaterEta * beta));

  // Hough transform histogram
  TH2F* hhough = new TH2F("hhough", "", 45, -pi, pi, 45, 0, pi);

  // Loop over PMTs
  for (int ipmt=0; ipmt<event->GetQPMTs()->GetEntries(); ipmt++) {
    QPMT* pmt = event->GetPMT(ipmt);

    // Assuming a straight line photon path (for now)...
    TVector3 vray = gSNO->GetPMTxyz()->GetXYZ(pmt->Getn()) - vvtx;

    // Times residual cut
    if (tres_cut > 0) {
      float ttof = vray.Mag() / (29.9792458 / kWaterEta);  // v in cm/ns
      float tres = pmt->Gett() - ftp->GetTime() - ttof;
      if (std::abs(tres) > tres_cut) {
        continue;
      }
    }

    // Inspired by Beier... insbeiered.
    TVector3 v;
    v.SetMagThetaPhi(1, vray.Theta() + tcher, vray.Phi());
    for (int ixf=0; ixf<kLoopIterations; ixf++) {
      v.Rotate(2 * pi / kLoopIterations, vray);
      hhough->Fill(v.Phi() * sin(v.Theta()), v.Theta());
    }
  }

  // Take the maximum bin in the transform space as the best fit
  int imax = hhough.GetMaximumBin();
  int nx  = hhough.GetXaxis()->GetNbins() + 2;
  int ny  = hhough.GetYaxis()->GetNbins() + 2;
  int ix = imax % nx;
  int iy = ((imax - ix) / nx) % ny;
  vfit.SetMagThetaPhi(1, hhough.GetYaxis()->GetBinCenter(iy),
                         hhough.GetXaxis()->GetBinCenter(ix));

  return hhough;
}

hough.h

/**
 * Modularized Hough transform, takes QEvent*, returns TH2F* of transform
 * space and the best fit direction (by reference).
 *
 * A. Mastbaum <mastbaum@hep.upenn.edu>, 11/2015
 */

#ifndef __HOUGH__
#define __HOUGH__

#include <TVector3.h>

class TH2F;
class QEvent;

// Constants
const float pi = 3.1415927;  // Math pi
const float kWaterEta = 1.34;  // Refractive index
const unsigned kLoopIterations = 100;  // Transform loop iterations


// Perform a Hough transform on the PMT hits.
//
// The caller assumes ownership of the returned pointer.
//
// @param event - Event data
// @param fit - TVector3 by reference, set to the best fit
// @param tres_cut - Symmetric time cut around prompt peak (ns, 0 for none)
// @param qweight - Weight hits by QHS in the transform space?
// @returns 2D histogram of the transform space, in phi/theta.
TH2F* hough(const QEvent* event, TVector3& vfit,
            float tres_cut=15.0, bool qweight=false);

#endif  // __HOUGH__