MoonShadowAnaHough_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //// Class: MoonShadowAna
3 //// Module Type: analyzer
4 //// File: MoonShadowAna_module.cc
5 ////
6 //////////////////////////////////////////////////////////////////////////
7 
8 
13 
21 #include "fhiclcpp/ParameterSet.h"
22 
27 
28 #include <iostream>
29 #include <set>
30 #include <cmath>
31 #include <algorithm>
32 #include <boost/math/special_functions.hpp>
33 #include <boost/math/special_functions/gamma.hpp>
34 #include "TNtuple.h"
35 
36 // To convert TDC to NOvA time
37 #include <NovaTimingUtilities/TimingUtilities.h>
38 
39 #include "Geometry/Geometry.h"
43 #include "MEFinder/MEClusters.h"
44 #include "MCCheater/BackTracker.h"
45 #include "RecoBase/CellHit.h"
46 #include "RecoBase/Track.h"
47 #include "RecoBase/Prong.h"
48 #include "RecoBase/Cluster.h"
49 #include "RecoBase/RecoHit.h"
50 #include "RecoBase/Vertex.h"
51 #include "Calibrator/Calibrator.h"
52 #include "Utilities/func/MathUtil.h"
53 #include "Simulation/FLSHit.h"
54 #include "Simulation/FLSHitList.h"
55 #include "Simulation/Particle.h"
59 
61 
62 #include <TH2D.h>
63 #include <iostream>
64 #include <stdlib.h>
65 #include <fstream>
66 
67 namespace moonshadowana {
68  class MoonShadowAnaHough ;
69 
70  const double PLANE_WIDTH = 6.6522; // cm
71  const double CELL_WIDTH = 3.9887; // cm
72  const double TDC_TICK = 15.625; // ns
73  const double SPEED_OF_LIGHT = 30.0; // cm / ns
74 
75  const unsigned N_CELLS = 384; // cells
76  const unsigned N_PLANES = 896; // planes
77 }
78 
80 public:
81  explicit MoonShadowAnaHough (fhicl::ParameterSet const & p);
82  virtual ~MoonShadowAnaHough ();
83 
84  void analyze(art::Event const & e) override;
85 
86  void beginJob() override;
87 
88  TVector3 vec(double zen, double azi);
89  double GetDot(double zen1, double azi1, double zen2, double azi2) const;
90  TVector3 AnglesToVector(double zen, double azi) const;
91  void InitLookupTable();
92  double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2);
93  bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2);
94  double GetRate(double zen, double azi);
95  double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2);
96  double GetSmartPrescale(double zen, double azi);
97  bool GetSmartPrescaleDecision(double prescale);
98  void GetHisto(std::string filename, int xmax, int ymax, TH2D* hist);
100 
101  void WriteEventInfo(const art::Event& event);
102  void WriteTrackInfo(const rb::Track& track, long timeStart);
103  void WriteHitListInfo(const art::PtrVector<rb::CellHit> hlist, long timeStart, TVector3 start, TVector3 end);
104  int GetX(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end);
105  int GetY(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end);
106  TVector3 GetDetVec(TVector3 vec) const;
107  double GetPhysLength(TVector3 vec) const;
108  TVector3 GetPhysVec(TVector3 vec) const;
109 
110 
111 private:
112  TNtuple* fNtuple;
115  // Lookup table for introducing "smart prescale"
117  int bin_size;
118 
119  // cut on angle
120  const double fDotCut = cos(5*M_PI/180);
121  // Limit on what rate is safe for detector hardware
122  double max_rate;
123  // We will cut more tracks on the border of 10 degree window than in the middle.
125  // file name of lookup table
127  // track 3d
129  // min track length
130  double fLengthCut;
131  // min # of hits in a track
132  double fBinCut;
133 
134  std::ofstream fOutputEvd;
135  bool init;
136 
137 };
138 
140  : EDAnalyzer (p),
141  //, fTrackModuleLabel (p.get<std::string>("TrackModuleLabel")),
142  max_rate(p.get<double>("AvgRate")),
143  angular_stdev(p.get<double>("AngularSuppression")),
145  fTrack3DLabel(p.get<std::string>("Track3DLabel")),
146  fLengthCut(p.get<double>("LengthCut")),
147  fBinCut(p.get<double>("BinCut"))
148 {
149  InitLookupTable();
150  bin_size = 1;
151 
152  init = false;
153 
154  //fOutputEvd.open("OFFEvd_Data.txt");
155 }
156 
158 {
159 fOutputEvd.close();
160 }
161 
163 {
165  fNtuple = tfs->make<TNtuple>("moon_ntuple", "Moon Shadow Ntuple", "zen_moon:azi_moon:zen_trk:azi_trk:nhit:len:dot:smart_pre:ang_pre");
166 }
167 
168  // write event info
170 {
171  if(!init)
172  {
173  fOutputEvd.open(std::string("OFFEvd_Data_") + std::string(std::to_string(event.run())) + std::string("_") + std::string(std::to_string(event.subRun())) + std::string(".txt"));
174  init = true;
175  }
176 
177  // Event info
178  fOutputEvd << "E " << event.event() << " ";
179  fOutputEvd << event.time().timeHigh() << std::endl;
180 
181  // Tracks
183  event.getByLabel(fTrack3DLabel, "", tracks);
184 
185  for(unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx)
186  {
187  const rb::Track& trk = (*tracks)[trkIdx];
188  if(trk.Is2D())continue;
189 
190 
191 
192  WriteTrackInfo(trk, event.time().timeHigh());
193  }
194 }
195 
197 {
198  TVector3 start = GetDetVec(track.Start());
199  TVector3 end = GetDetVec(track.Stop());
200  if(start.Y() < end.Y()) std::swap(start, end); // all tracks are down-going
201 
202  //Check whether 3D track vector makes sense
203  if(
204  (start.X() == 0 && end.X() == 0) ||
205  (start.Y() == 0 && end.Y() == 0) ||
206  (start.Z() == 0 && end.Z() == 0))
207  {
208  return;
209  }
210 
211  // Track info
212  fOutputEvd << "T ";
213  fOutputEvd << "(" << start.X() << "," << start.Y() << "," << start.Z() << ") ";
214  fOutputEvd << "(" << end.X() << "," << end.Y() << "," << end.Z() << ")";
216 
217  art::PtrVector<rb::CellHit> hlist = track.AllCells();
218 
219  WriteHitListInfo(hlist, timeStart, start, end);
220 
221 }
223 {
224  for (size_t j=0; j < hlist.size(); ++j) {
225  art::Ptr<rb::CellHit> h = hlist.at(j);
226 
227  if(j != 0)
228  fOutputEvd << ",";
229 
230  int view = (int)h->View();
231  fOutputEvd << ((long)h->TDC())/64 << "/";
232  fOutputEvd << (view == 0 ? (int)h->Cell() : GetX(h, start, end)) << "/";
233  fOutputEvd << (view == 1 ? (int)h->Cell() : GetY(h, start, end)) << "/";
234  fOutputEvd << (int)h->Plane() << "/";
235  fOutputEvd << view + 1;
236  }
237 
239 }
241 {
242  int z = h->Plane();
243 
244  double xSpread = end.X() - start.X();
245  double zSpread = end.Z() - start.Z();
246 
247  double zSpreadLoc = z - start.Z();
248 
249  double x = start.X() + xSpread * zSpreadLoc/zSpread;
250 
251  return round(x);
252 }
254 {
255  int z = h->Plane();
256 
257  double ySpread = end.Y() - start.Y();
258  double zSpread = end.Z() - start.Z();
259 
260  double zSpreadLoc = z - start.Z();
261 
262  double y = start.Y() + ySpread * zSpreadLoc/zSpread;
263 
264  return round(y);
265 }
266 
268 {
269  int x = int((vec.X() + CELL_WIDTH * N_CELLS /2) / CELL_WIDTH);
270  int y = int((vec.Y() + CELL_WIDTH * N_CELLS /2) / CELL_WIDTH);
271  int z = int(vec.Z() / PLANE_WIDTH);
272 
273  return TVector3(x,y,z);
274 }
275 
277  {
278  double x = vec.X() * CELL_WIDTH;
279  double y = vec.Y() * CELL_WIDTH;
280  double z = vec.Z() * PLANE_WIDTH;
281 
282  return sqrt(x*x + y*y + z*z);
283  }
285 {
286  double x = vec.X() * CELL_WIDTH;
287  double y = vec.Y() * CELL_WIDTH;
288  double z = vec.Z() * PLANE_WIDTH;
289 
290  return TVector3(x,y,z);
291 }
292 
294 
295  // Grab the tracks from the event
297  e.getByLabel(fTrack3DLabel, tracks);
298 
299  // We need the hitlist just so we can get the times for the trigger
301 
302  art::Timestamp ts = e.time();
303  time_t timeSec = ts.timeHigh();
304 
305  double zen_moon, azi_moon, zen_sun, azi_sun;
306  loc->GetMoonPosition_FD(timeSec,zen_moon, azi_moon);
307  loc->GetSunPosition_FD(timeSec,zen_sun, azi_sun);
308 
309  float vars[9];
310  vars[0] = zen_moon; // zen_moon
311  vars[1] = azi_moon; // azi_moon
312 
313  double smartPrescale = GetSmartPrescale(zen_moon, azi_moon, zen_sun, azi_sun);
314 
315  double cos_trk_moon, cos_trk_sun;
316  int trackIndex = 0;
317  for(const novaddt::Track3D& track: *tracks)
318  {
319  const novaddt::HitList& hlist = *hls.at(trackIndex);
320 
321 trackIndex++;
322 
323  TVector3 start = track.Start();
324  TVector3 end = track.End();
325  if(start.Y() < end.Y()) std::swap(start, end);
326 
327  // WORK
328  //Check whether 3D track vector makes sense
329  if(
330  (start.X() == 0 && end.X() == 0) ||
331  (start.Y() == 0 && end.Y() == 0) ||
332  (start.Z() == 0 && end.Z() == 0))
333  {
334  continue;
335  }
336 
337  // Apply length cuts
338  if(GetPhysLength(start-end) < fLengthCut || (int)hlist.size() < fBinCut)
339  {
340  continue;
341  }
342 
343  // Point back along the track towards the moon
344  const TVector3 dir = GetPhysVec((start-end)).Unit();
345 
346  const double zen_trk = acos(dir.Y()) * 180 / M_PI;
347  // docdb 5485 / numix 17 says FD is oriented -27o51'26'' from true North
348  // azi_trk = atan2(-dir.X(), dir.Z()) * 180 / M_PI - 27.857222;
349  // Or: 332o03'58.071769" (clockwise from North). Email from Virgil Bocean
350  // to Alec Habig.
351  double azi_trk = atan2(-dir.X(), dir.Z()) * 180 / M_PI + 332.142778;
352  if(azi_trk < 0) azi_trk += 360;
353  if(azi_trk > 360) azi_trk -= 360;
354 
355  cos_trk_moon = vec(zen_moon, azi_moon).Dot(vec(zen_trk,azi_trk));
356  cos_trk_sun = vec(zen_sun, azi_sun).Dot(vec(zen_trk,azi_trk));
357  (void)cos_trk_sun; // Suppresss "unused variable" warnings
358 
359  vars[2] = zen_trk; // zen_trk
360  vars[3] = azi_trk; // azi_trk
361  vars[4] = hlist.size(); // nhit
362  vars[5] = GetPhysLength(start-end); // len
363  vars[6] = cos_trk_moon; // cosine angle trk and moon
364 
365  double angularPrescale = GetAngularPrescale(zen_moon, azi_moon, zen_trk, azi_trk);
366 
367  vars[7] = smartPrescale;
368  vars[8] = angularPrescale;
369 
370  fNtuple->Fill(vars);
371 
372  } // end of tracks
373 
374 }
375 
376 TVector3 moonshadowana::MoonShadowAnaHough::vec(double zen, double azi){
377  zen *= M_PI / 180;
378  azi *= M_PI / 180;
379  return TVector3(sin(zen)*cos(azi), sin(zen)*sin(azi), cos(zen));
380 }
381 
382  double moonshadowana::MoonShadowAnaHough::GetRate(double zen, double azi)
383  {
384  if(zen > 90)
385  return 0;
386  // Smart prescale is ratio of maximal rate to expected rate divided by angular prescale factor for given position
387  int bin_azi = round(azi/bin_size);
388  int bin_zen = round(zen/bin_size);
389 
390  //return 0;
391  return rate_collected_angular->GetBinContent(bin_azi, bin_zen);
392  }
393 
394  double moonshadowana::MoonShadowAnaHough::GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)
395  {
396  double dTheta = acos(GetDot(zen1, azi1, zen2, azi2)) * 180 / M_PI;
397  return exp(-dTheta*dTheta/(angular_stdev*angular_stdev));
398  }
399  bool moonshadowana::MoonShadowAnaHough::GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
400  {
401  double prescale = GetAngularPrescale(zen1, azi1, zen2, azi2);
402 
403  double r = (double)rand() / (double)RAND_MAX;
404 
405  return r < prescale;
406  }
407 
409  {
410  // Init lookup table and make 2d hist with expected rate
411  // Get histogram with full rates in 2degx2deg bins
412  int azi_max = 360;
413  int zen_max = 90;
414 
415  rate_collected_angular = new TH2D("LookupPrescale", "", azi_max, 0, 360, zen_max, 0, 90);
416 
417  GetHisto(fLookupFile, azi_max, zen_max, rate_collected_angular);
418  }
419 
421 {
422  std::ifstream data;
423  data.open(filename.c_str());
424  for(int j=0; j<ymax;j++)
425  {
426  for(int i=0; i<xmax; i++)
427  {
428  double val;
429  data >> val;
430  hist->SetBinContent(i,j, val);
431  }
432  }
434  TH2D* htmp = tfs->make<TH2D>("LookupPrescale_stored", "", 360, 0, 360, 90, 0, 90);
435  htmp = (TH2D*)hist->Clone();
436  htmp->Write();
437  data.close();
438 }
439 
440 double moonshadowana::MoonShadowAnaHough::GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
441  {
442  double rate1 = GetRate(zen1, azi1);
443  double rate2 = GetRate(zen2, azi2);
444 
445  if(rate1 == 0 && rate2 == 0)
446  return 0;
447 
448  return max_rate / (rate1 + rate2);
449  }
451  {
452  double rate = GetRate(zen, azi);
453 
454  if(rate == 0)
455  return 0;
456 
457  return max_rate / rate;
458  }
459 
461  {
462  // Is it good enough for small prescale factors ~1/100?
463  double r = (double)rand() / (double)RAND_MAX;
464  return r < prescale;
465  }
466 
467 double moonshadowana::MoonShadowAnaHough::GetDot(double zen1, double azi1, double zen2, double azi2) const
468  {
469  if(zen1 == zen2 && azi1 == azi2)
470  return 1;
471 
472  zen1 *= M_PI / 180;
473  zen2 *= M_PI / 180;
474  azi1 *= M_PI / 180;
475  azi2 *= M_PI / 180;
476  double dot = sin(zen1)*cos(azi1)*sin(zen2)*cos(azi2) + sin(zen1)*sin(azi1)*sin(zen2)*sin(azi2) + cos(zen1)*cos(zen2);
477 
478  return dot;
479  }
480 
481 // find lookup table file
483 {
484  cet::search_path sp("FW_SEARCH_PATH");
485  std::string lookupFname;
486  sp.find_file(("Eval/PrescaleLookup.txt"), lookupFname);
487  assert(!lookupFname.empty());
488  return lookupFname;
489 }
490 
int GetY(art::Ptr< rb::CellHit > h, TVector3 start, TVector3 end)
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
std::map< std::string, double > xmax
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
const double CELL_WIDTH
const double TDC_TICK
bool Is2D() const
Definition: Cluster.h:96
void WriteHitListInfo(const art::PtrVector< rb::CellHit > hlist, long timeStart, TVector3 start, TVector3 end)
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
std::vector< DAQHit > HitList
Definition: HitList.h:15
void GetMoonPosition_FD(time_t time, double &zen, double &azi)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
Definition: event.h:19
TVector3 GetPhysVec(TVector3 vec) const
string filename
Definition: shutoffs.py:106
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
T acos(T number)
Definition: d0nt_math.hpp:54
std::string find_file(std::string const &filename) const
#define M_PI
Definition: SbMath.h:34
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double GetRate(double zen, double azi)
Double_t ymax
Definition: plot.C:25
Locate positions of celestial bodies.
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
double GetDot(double zen1, double azi1, double zen2, double azi2) const
const unsigned N_CELLS
void WriteEventInfo(const art::Event &event)
const double SPEED_OF_LIGHT
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
void GetHisto(std::string filename, int xmax, int ymax, TH2D *hist)
int GetX(art::Ptr< rb::CellHit > h, TVector3 start, TVector3 end)
const std::map< std::pair< std::string, std::string >, Variable > vars
const double PLANE_WIDTH
const double j
Definition: BetheBloch.cxx:29
Eigen::VectorXd vec
reference at(size_type n)
Definition: PtrVector.h:365
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void GetSunPosition_FD(time_t time, double &zen, double &azi)
Vertex location in position and time.
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TVector3 vec(double zen, double azi)
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
void analyze(art::Event const &e) override
void WriteTrackInfo(const rb::Track &track, long timeStart)
TVector3 GetDetVec(TVector3 vec) const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TVector3 AnglesToVector(double zen, double azi) const
art::ServiceHandle< locator::CelestialLocator > loc
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
Timestamp time() const
Definition: Event.h:61
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
const unsigned N_PLANES
typedef void(XMLCALL *XML_ElementDeclHandler)(void *userData
MoonShadowAnaHough(fhicl::ParameterSet const &p)
Encapsulate the geometry of one entire detector (near, far, ndos)
double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
T atan2(T number)
Definition: d0nt_math.hpp:72
bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)