MoonShadow_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MoonShadow
3 // Module Type: filter
4 // \file: MoonShadow_module.cc
5 // \brief Select tracks originating from the direction of the moon or sun
6 // \author Christopher Backhouse - bckhouse@caltech.edu
7 ////////////////////////////////////////////////////////////////////////
8 
14 
19 
20 // To convert TDC to NOvA time
21 #include <NovaTimingUtilities/TimingUtilities.h>
22 // To convert vectors to lenght
24 
25 // This will ultimately come from an external UPS product
27 
28 //
29 #include <TH2D.h>
30 #include <iostream>
31 #include <stdlib.h>
32 #include <fstream>
33 
34 namespace novaddt
35 {
36  class MoonShadow: public art::EDFilter
37  {
38  public:
39  explicit MoonShadow(const fhicl::ParameterSet& p);
40  virtual ~MoonShadow();
41 
42  virtual bool filter(art::Event & e);
43 
44  virtual void beginJob();
45  virtual void endJob();
46 
47  // Why is it protected
48  protected:
51 
52  // Angle manipulation
53  double inline RadToDeg(double rad) const;
54  double inline DegToRad(double deg) const;
55 
56  // Vector manipulation
57  TVector3 AnglesToVector(double zen, double azi) const;
58  double GetDot(double zen1, double azi1, double zen2, double azi2) const;
59  TVector3 FarDetToAstron(TVector3& fdVector) const;
60  TVector3 AstronToFarDet(TVector3& astrVector) const;
61  double GetPhysLength(TVector3 vec) const;
62  double GetEffectiveArea();
63  void InitLookupTable();
64  double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2);
65  bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2);
66  double GetRate(double zen, double azi);
67  double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2);
68  double GetSmartPrescale(double zen, double azi);
69  bool GetSmartPrescaleDecision(double prescale);
70  void GetHisto(std::string filename, int xmax, int ymax, TH2D* hist);
71 
72  TVector3 GetPhysVec(TVector3 vec) const;
73 
74 
75  // Time
76  time_t GetTime(art::Event& e);
77 
78  int fPrescale;
79 
83  double fDotCut;
84  double fLengthCut;
85  int fBinCut;
86 
88 
90 
91  // Lookup table for introducing "smart prescale"
93  int bin_size;
94 
95  // Limit on what rate is safe for detector hardware
96  double max_rate;
97  // We will cut more tracks on the border of 10 degree window than in the middle.
98  double angular_stdev;
99  // file name of lookup table
101 
102  };
103 
104 
106  : fPrescale(pset.get<int>("prescale")),
107  fTrack3DLabel(pset.get<std::string>("Track3DLabel")),
108  fInputHitLabel(pset.get<std::string>("InputHitLabel")),
109  fInputHitInstance(pset.get<std::string>("InputHitInstance")),
110  fDotCut(cos(5*M_PI/180) /*cos(pset.get<double>("AngleCut")*M_PI/180)*/),
111  fLengthCut(pset.get<double>("LengthCut")),
112  fBinCut(pset.get<int>("BinCut")),
114  fNEvents(0), fNTracks(0), fNSunTrigs(0), fNMoonTrigs(0),
115  max_rate(pset.get<double>("AvgRate")),
116  angular_stdev(pset.get<double>("AngularSuppression")),
118  {
119  produces< std::vector<TriggerDecision>>();
120  InitLookupTable();
121  srand(time(NULL));
122 
123  bin_size = 1;
124  }
125 
126  double MoonShadow::GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)
127  {
128  double dTheta = acos(GetDot(zen1, azi1, zen2, azi2)) * 180 / M_PI;
129  return exp(-dTheta*dTheta/(angular_stdev*angular_stdev));
130  }
131  bool MoonShadow::GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
132  {
133  double prescale = GetAngularPrescale(zen1, azi1, zen2, azi2);
134 
135  double r = (double)rand() / (double)RAND_MAX;
136 
137  return r < prescale;
138  }
139 
141 {
142  std::ifstream data;
143  data.open(filename.c_str());
144  for(int j=0; j<ymax;j++)
145  {
146  for(int i=0; i<xmax; i++)
147  {
148  double val;
149  data >> val;
150  hist->SetBinContent(i,j, val);
151  }
152  }
153 
154  data.close();
155 }
156 
157 
159  {
160  // Init lookup table and make 2d hist with expected rate
161  // Get histogram with full rates in 2degx2deg bins
162  int azi_max = 360;
163  int zen_max = 90;
164 
165  rate_collected_angular = new TH2D("LookupPrescale", "", azi_max, 0, 360, zen_max, 0, 90);
166 
167  GetHisto(fLookupFile, azi_max, zen_max, rate_collected_angular);
168  }
169  // Prescale factor taken for sum of rates in 2 positions
170  double MoonShadow::GetRate(double zen, double azi)
171  {
172  if(zen > 90)
173  return 0;
174  // Smart prescale is ratio of maximal rate to expected rate divided by angular prescale factor for given position
175  int bin_azi = round(azi/bin_size);
176  int bin_zen = round(zen/bin_size);
177 
178  //return 0;
179  return rate_collected_angular->GetBinContent(bin_azi, bin_zen);
180  }
181  double MoonShadow::GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
182  {
183  double rate1 = GetRate(zen1, azi1);
184  double rate2 = GetRate(zen2, azi2);
185 
186  if(rate1 == 0 && rate2 == 0)
187  return 0;
188 
189  return max_rate / (rate1 + rate2);
190  }
191  double MoonShadow::GetSmartPrescale(double zen, double azi)
192  {
193  double rate = GetRate(zen, azi);
194 
195  if(rate == 0)
196  return 0;
197 
198  return max_rate / rate;
199  }
201  {
202  // Is it good enough for small prescale factors ~1/100?
203  double r = (double)rand() / (double)RAND_MAX;
204  return r < prescale;
205  }
206 
208  {
209  }
210 
212  {
213  }
214 
215  // Conversion radians <-> degrees
216  double inline MoonShadow::DegToRad(double deg) const
217  {
218  return deg * M_PI / 180;
219  }
220  double inline MoonShadow::RadToDeg(double rad) const
221  {
222  return rad * 180 / M_PI;
223  }
224 
225  // ???
227  {
228  cet::search_path sp("FW_SEARCH_PATH");
229  std::string ephemerisFname;
230  sp.find_file(("MoonShadowTrigger/novas/JPLEPH"), ephemerisFname);
231  assert(!ephemerisFname.empty());
232  return ephemerisFname;
233  }
234  // find lookup table file
236  {
237  cet::search_path sp("FW_SEARCH_PATH");
238  std::string lookupFname;
239  sp.find_file(("MoonShadowTrigger/smart_prescale/PrescaleLookup.txt"), lookupFname);
240  assert(!lookupFname.empty());
241  return lookupFname;
242  }
243 
244 
245  // Creates unit vector pointing to a point on unit sphere with theta=zen, phi = azi
246  TVector3 MoonShadow::AnglesToVector(double zen, double azi) const
247  {
248  zen *= M_PI / 180;
249  azi *= M_PI / 180;
250  // Doesn't matter which axis is which
251  return TVector3(sin(zen)*cos(azi), sin(zen)*sin(azi), cos(zen));
252  }
253  // Gets dot product of unit vectors with given angles
254  double MoonShadow::GetDot(double zen1, double azi1, double zen2, double azi2) const
255  {
256  if(zen1 == zen2 && azi1 == azi2)
257  return 1;
258 
259  zen1 *= M_PI / 180;
260  zen2 *= M_PI / 180;
261  azi1 *= M_PI / 180;
262  azi2 *= M_PI / 180;
263  double dot = sin(zen1)*cos(azi1)*sin(zen2)*cos(azi2) + sin(zen1)*sin(azi1)*sin(zen2)*sin(azi2) + cos(zen1)*cos(zen2);
264 
265  return dot;
266  }
267 
268  // Transfoms vector in FarDet coordinate system to vector in astronomical coordinate system
269  TVector3 MoonShadow::FarDetToAstron(TVector3& fdVector) const
270  {
271  throw "Not implemented";
272 
273  }
274  // Transfoms vector in astronomical coordinate system to vector in FarDet coordinate system
275  TVector3 MoonShadow::AstronToFarDet(TVector3& fdVector) const
276  {
277  throw "Not implemented";
278 
279  }
280  // Get length of the vector in cm
281  double MoonShadow::GetPhysLength(TVector3 vec) const
282  {
283  double x = vec.X() * novaddt::smt::Constants::CELL_WIDTH;
284  double y = vec.Y() * novaddt::smt::Constants::CELL_WIDTH;
285  double z = vec.Z() * novaddt::smt::Constants::PLANE_WIDTH;
286 
287  return sqrt(x*x + y*y + z*z);
288  }
289 
290 
291  // Gets absolute time of the event. This function shoudl be used to determine time up to uncertainty bigger than hit time
292  // difference
294  {
296 
298 
299  uint64_t secs = header->timeStart/novadaq::timeutils::NOVA_TIME_FACTOR;
300 
301  uint64_t nova_time = novadaq::timeutils::NOVA_EPOCH + secs;
302 
303  return time_t(nova_time);
304  }
305 
306  TVector3 MoonShadow::GetPhysVec(TVector3 vec) const
307 {
308  double x = vec.X() * novaddt::smt::Constants::CELL_WIDTH;
309  double y = vec.Y() * novaddt::smt::Constants::CELL_WIDTH;
310  double z = vec.Z() * novaddt::smt::Constants::PLANE_WIDTH;
311 
312  return TVector3(x,y,z);
313 }
314 
315 
317  {
318  ++fNEvents;
319 
320  std::unique_ptr<std::vector<TriggerDecision>> td(new std::vector<TriggerDecision>);
321 
322  // Grab the tracks from the event
324  evt.getByLabel(fTrack3DLabel, tracks);
325 
326  // We need the hitlist just so we can get the times for the trigger
327  art::FindOneP<HitList> hls(tracks, evt, fTrack3DLabel);
328 
329  //time_t t = GetTime(hls);
330  time_t t = GetTime(evt);
331  // std::cout << "Time_t: " << asctime(gmtime(&t)) << std::endl;
332 
333  double zen_moon, azi_moon;
334  fLocate.GetMoonPosition(t, zen_moon, azi_moon);
335  double zen_sun, azi_sun;
336  fLocate.GetSunPosition(t, zen_sun, azi_sun);
337 
338  int trackIndex = 0;
339 
340  double smartPrescale = GetSmartPrescale(zen_moon, azi_moon, zen_sun, azi_sun);
341  double smartPrescaleSun = GetSmartPrescale(zen_sun, azi_sun);
342  double smartPrescaleMoon = GetSmartPrescale(zen_moon, azi_moon);
343  double smartPrescaleTot = smartPrescaleSun + smartPrescaleMoon;
344  smartPrescaleSun = smartPrescale * (smartPrescaleSun/smartPrescaleTot);
345  smartPrescaleMoon = smartPrescale * (smartPrescaleMoon/smartPrescaleTot);
346 
347  for(const Track3D& track: *tracks){
348  const HitList& hlist = *hls.at(trackIndex);
349  ++trackIndex;
350 
351  TVector3 start = track.Start();
352  TVector3 end = track.End();
353  if(start.Y() < end.Y()) std::swap(start, end);
354 
355  // WORK
356  //Check whether 3D track vector makes sense
357  if(
358  (start.X() == 0 && end.X() == 0) ||
359  (start.Y() == 0 && end.Y() == 0) ||
360  (start.Z() == 0 && end.Z() == 0))
361  {
362  continue;
363  }
364 
365 
366 
367  // Apply length cuts
368  if(GetPhysLength(start-end) < fLengthCut || (int)hlist.size() < fBinCut)
369  {
370  continue;
371  }
372 
373 
374  ++fNTracks;
375 
376 
377 
378 
379 
380  // Point back along the track towards the moon
381  const TVector3 dir = GetPhysVec((start-end)).Unit();
382 
383  const double zen_trk = acos(dir.Y()) * 180 / M_PI;
384  // docdb 5485 / numix 17 says FD is oriented -27o51'26'' from true North
385  // azi_trk = atan2(-dir.X(), dir.Z()) * 180 / M_PI - 27.857222;
386  // Or: 332o03'58.071769" (clockwise from North). Email from Virgil Bocean
387  // to Alec Habig.
388  double azi_trk = atan2(-dir.X(), dir.Z()) * 180 / M_PI + 332.142778;
389  if(azi_trk < 0) azi_trk += 360;
390  if(azi_trk > 360) azi_trk -= 360;
391 
392  const double dotMoon = GetDot(zen_moon, azi_moon, zen_trk, azi_trk);
393  const double dotSun = GetDot(zen_sun, azi_sun, zen_trk, azi_trk);
394 
395  // Inline &&'s call all these functions, so this is made to avoid extra calling
396  bool trigMoon = (dotMoon > fDotCut);
397  if(trigMoon) { trigMoon = trigMoon && GetSmartPrescaleDecision(smartPrescale); }
398  if(trigMoon) { trigMoon = trigMoon && GetAngularPrescaleDecision(zen_moon, azi_moon, zen_trk, azi_trk); }
399 
400  bool trigSun = (dotSun > fDotCut);
401  if(trigSun) { trigSun = trigSun && GetSmartPrescaleDecision(smartPrescale); }
402  if(trigSun) { trigSun = trigSun && GetAngularPrescaleDecision(zen_sun, azi_sun, zen_trk, azi_trk); }
403 
404 
405  if(trigMoon || trigSun){
406  auto mint = hlist.front().TDC().val;
407  auto maxt = hlist.back().TDC().val;
408  for(const DAQHit& it: hlist){
409  mint = std::min(mint, it.TDC().val);
410  maxt = std::max(maxt, it.TDC().val);
411  }
412 
413  // Add 12.5us padding to each end of the window
414  mint -= 800;
415  maxt += 800;
416 
417  if(trigMoon){
418  ++fNMoonTrigs;
419  td->emplace_back(mint, maxt-mint,
421  fPrescale);
422  }
423  if(trigSun){
424  ++fNSunTrigs;
425  td->emplace_back(mint, maxt-mint,
427  fPrescale);
428  }
429  }
430 
431  } // end of tracks
432 
433 
434 
435  const bool ret = !td->empty();
436  evt.put(std::move(td));
437  return ret;
438  }
439 
441  {
442  delete rate_collected_angular;
443  std::cout << "=== novaddt::MoonShadow endJob " << std::endl;
444  std::cout << "\tNumber of events: " << fNEvents << std::endl;
445  std::cout << "\tNumber of tracks: " << fNTracks << std::endl;
446  std::cout << "\tNumber of sun triggers: " << fNSunTrigs << std::endl;
447  std::cout << "\tNumber of moon triggers: " << fNMoonTrigs << std::endl;
448  }
449 
450 
452 
453 } // namespace
T max(const caf::Proxy< T > &a, T b)
virtual bool filter(art::Event &e)
static constexpr Double_t deg
Definition: Munits.h:165
const uint32_t NOVA_EPOCH
std::map< std::string, double > xmax
double RadToDeg(double rad) const
set< int >::iterator it
double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
double GetRate(double zen, double azi)
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string FindEphemerisFile() const
Definition: event.h:19
double GetDot(double zen1, double azi1, double zen2, double azi2) const
const double CELL_WIDTH
Definition: Constants.h:12
double GetEffectiveArea()
string filename
Definition: shutoffs.py:106
TVector3 AstronToFarDet(TVector3 &astrVector) const
const uint64_t NOVA_TIME_FACTOR
DEFINE_ART_MODULE(TestTMapFile)
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
T acos(T number)
Definition: d0nt_math.hpp:54
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
std::string find_file(std::string const &filename) const
#define M_PI
Definition: SbMath.h:34
TVector3 AnglesToVector(double zen, double azi) const
bool GetSmartPrescaleDecision(double prescale)
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Double_t ymax
Definition: plot.C:25
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void GetMoonPosition(time_t t, double &zen, double &azi)
double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)
shadow::NOVASLocate fLocate
std::string FindSmartPrescaleLookup() const
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
time_t GetTime(art::Event &e)
const double j
Definition: BetheBloch.cxx:29
Eigen::VectorXd vec
static constexpr Double_t rad
Definition: Munits.h:162
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
void GetSunPosition(time_t t, double &zen, double &azi)
TVector3 GetPhysVec(TVector3 vec) const
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
std::string fInputHitInstance
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
MoonShadow(const fhicl::ParameterSet &p)
bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
double DegToRad(double deg) const
void GetHisto(std::string filename, int xmax, int ymax, TH2D *hist)
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
TVector3 FarDetToAstron(TVector3 &fdVector) const
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
double GetPhysLength(TVector3 vec) const
const double PLANE_WIDTH
Definition: Constants.h:11
T atan2(T number)
Definition: d0nt_math.hpp:72