28 #include <boost/math/special_functions.hpp> 29 #include <boost/math/special_functions/gamma.hpp> 33 #include <NovaTimingUtilities/TimingUtilities.h> 48 #include "Utilities/func/MathUtil.h" 84 TVector3
vec(
double zen,
double azi);
85 double GetDot(
double zen1,
double azi1,
double zen2,
double azi2)
const;
90 double GetRate(
double zen,
double azi);
149 fNtuple = tfs->
make<TNtuple>(
"moon_ntuple",
"Moon Shadow Ntuple",
"zen_moon:azi_moon:zen_trk:azi_trk:nhit:len:dot:smart_pre:ang_pre");
169 for(
unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx)
172 if(trk.
Is2D())
continue;
184 if(start.Y() < end.Y())
std::swap(start, end);
188 (start.X() == 0 && end.X() == 0) ||
189 (start.Y() == 0 && end.Y() == 0) ||
190 (start.Z() == 0 && end.Z() == 0))
197 fOutputEvd <<
"(" << start.X() <<
"," << start.Y() <<
"," << start.Z() <<
") ";
198 fOutputEvd <<
"(" << end.X() <<
"," << end.Y() <<
"," << end.Z() <<
")";
208 for (
size_t j=0;
j < hlist.
size(); ++
j) {
228 double xSpread = end.X() - start.X();
229 double zSpread = end.Z() - start.Z();
231 double zSpreadLoc = z - start.Z();
233 double x = start.X() + xSpread * zSpreadLoc/zSpread;
241 double ySpread = end.Y() - start.Y();
242 double zSpread = end.Z() - start.Z();
244 double zSpreadLoc = z - start.Z();
246 double y = start.Y() + ySpread * zSpreadLoc/zSpread;
257 return TVector3(x,y,z);
262 if(e.
event() == 486130)
270 double zen_moon, azi_moon, zen_sun, azi_sun;
282 double smartPrescale =
GetSmartPrescale(zen_moon, azi_moon, zen_sun, azi_sun);
284 double zen_trk, azi_trk, cos_trk_moon, cos_trk_sun;
285 for(
unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx){
287 if(trk.
Is2D())
continue;
293 if(start.Y() < end.Y())
std::swap(start, end);
301 zen_trk =
acos(dir.Y()) * 180 /
M_PI;
302 azi_trk =
atan2(-dir.X(), dir.Z()) * 180 /
M_PI + 332.066111;
303 if(azi_trk < 0) azi_trk += 360;
304 if(azi_trk > 360) azi_trk -= 360;
306 cos_trk_moon =
vec(zen_moon, azi_moon).Dot(
vec(zen_trk,azi_trk));
307 cos_trk_sun =
vec(zen_sun, azi_sun).Dot(
vec(zen_trk,azi_trk));
312 vars[4] = trk.
NCell();
314 vars[6] = cos_trk_moon;
318 vars[7] = smartPrescale;
319 vars[8] = angularPrescale;
353 double dTheta =
acos(
GetDot(zen1, azi1, zen2, azi2)) * 180 /
M_PI;
360 double r = (double)
rand() / (double)RAND_MAX;
380 data.open(filename.c_str());
387 hist->SetBinContent(
i,
j, val);
391 TH2D* htmp = tfs->
make<TH2D>(
"LookupPrescale_stored",
"", 360, 0, 360, 90, 0, 90);
392 htmp = (TH2D*)hist->Clone();
399 double rate1 =
GetRate(zen1, azi1);
400 double rate2 =
GetRate(zen2, azi2);
402 if(rate1 == 0 && rate2 == 0)
409 double rate =
GetRate(zen, azi);
420 double r = (double)
rand() / (double)RAND_MAX;
426 if(zen1 == zen2 && azi1 == azi2)
443 sp.
find_file((
"Eval/PrescaleLookup.txt"), lookupFname);
444 assert(!lookupFname.empty());
SubRunNumber_t subRun() const
int GetX(art::Ptr< rb::CellHit > h, TVector3 start, TVector3 end)
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
int GetY(art::Ptr< rb::CellHit > h, TVector3 start, TVector3 end)
std::map< std::string, double > xmax
int32_t TDC() const
The time of the last baseline sample.
unsigned short Plane() const
void GetMoonPosition_FD(time_t time, double &zen, double &azi)
constexpr std::uint32_t timeHigh() const
bool GetSmartPrescaleDecision(double prescale)
A rb::Prong with full reconstructed trajectory.
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
void WriteEventInfo(const art::Event &event)
fvar< T > round(const fvar< T > &x)
std::string find_file(std::string const &filename) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
const XML_Char const XML_Char * data
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Locate positions of celestial bodies.
virtual double TotalLength() const
Length (cm) of all the track segments.
unsigned short Cell() const
Track finder for cosmic rays.
TH2D * rate_collected_angular
double GetDot(double zen1, double azi1, double zen2, double azi2) const
const double SPEED_OF_LIGHT
fvar< T > exp(const fvar< T > &x)
MoonShadowAna(fhicl::ParameterSet const &p)
art::ServiceHandle< locator::CelestialLocator > loc
EventNumber_t event() const
reference at(size_type n)
EDAnalyzer(Table< Config > const &config)
void GetSunPosition_FD(time_t time, double &zen, double &azi)
TVector3 AnglesToVector(double zen, double azi) const
Vertex location in position and time.
double dot(const std::vector< double > &x, const std::vector< double > &y)
std::string fTrackModuleLabel
std::string FindSmartPrescaleLookup() const
T * make(ARGS...args) const
void GetHisto(std::string filename, int xmax, int ymax, TH2D *hist)
double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
void analyze(art::Event const &e) override
bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
TVector3 vec(double zen, double azi)
assert(nhit_max >=nhit_nbins)
TVector3 GetDetVec(TVector3 vec) const
std::string to_string(ModuleType mt)
TVector3 Stop() const
Position of the final trajectory point.
void WriteTrackInfo(const rb::Track &track, long timeStart)
typedef void(XMLCALL *XML_ElementDeclHandler)(void *userData
Encapsulate the geometry of one entire detector (near, far, ndos)
void WriteHitListInfo(const art::PtrVector< rb::CellHit > hlist, long timeStart, TVector3 start, TVector3 end)
double GetRate(double zen, double azi)
double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)