Skip to content

Instantly share code, notes, and snippets.

@mastbaum
Created Dec 4, 2015
Embed
What would you like to do?
Hough transform for QEvents
#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;
}
/**
* 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__
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment