#include #include #include #include "TRandom3.h" #include "JMath/JConstants.hh" #include "JGeometry3D/JVector3D.hh" #include "JGeometry3D/JVersor3D.hh" #include "JGeometry3D/JAxis3D.hh" #include "JGeometry3D/JTrack3D.hh" #include "JGeometry3D/JRotation3D.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to test angle of incidence of a photon from a muon on a PMT. * \author mdejong */ int main(int argc, char**argv) { using namespace std; using namespace JPP; JAxis3D pmt; int numberOfEvents; double precision; UInt_t seed; int debug; try { JParser<> zap("Example program to test angle of incidence of a photon from a muon on a PMT."); zap['n'] = make_field(numberOfEvents) = 1000; zap['P'] = make_field(pmt) = JAxis3D(JVector3D(0,0,1), JVersor3D(0,0,1)); zap['e'] = make_field(precision) = 1.0e-6; zap['S'] = make_field(seed) = 0; zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } gRandom->SetSeed(seed); const JVersor3D U(getSinThetaC(), 0.0, getCosThetaC()); // photon direction in the reference frame of muon track int number_of_errors = 0; int number_of_events = 0; for (int i = 0; i != numberOfEvents; ++i) { STATUS(setw(4) << i << '\r'); DEBUG(endl); const double x = gRandom->Uniform(-1.0, +1.0); const double y = gRandom->Uniform(-1.0, +1.0); const double z = gRandom->Uniform(-1.0, +1.0); const double dx = gRandom->Uniform(-1.0, +1.0); const double dy = gRandom->Uniform(-1.0, +1.0); const double dz = gRandom->Uniform(-1.0, +1.0); JTrack3D ta(JVector3D(x, y, z), JVersor3D(dx, dy, dz), 0.0); JAxis3D tb(pmt); tb.transform(ta); // transform PMT to reference frame of muon track if (ta.getDistance(pmt) > precision) { ++number_of_events; // method JTrack3D::getDot(..) calculates angle of incidence of photon on PMT if (fabs(tb.getDirection().getDot(U) - ta.getDot(pmt)) > precision) { ++number_of_errors; } } DEBUG("distance " << ta.getDistance(pmt) << endl); DEBUG("dot (1) " << tb.getDirection().getDot(U) << endl); DEBUG("dot (2) " << ta.getDot(pmt) << endl); } STATUS(endl); NOTICE("Number of errors: " << number_of_errors << " out of " << number_of_events << " events." << endl); ASSERT(number_of_events > 0); ASSERT(number_of_errors == 0); return 0; }