MoonShadowAna_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 
18 
22 #include "fhiclcpp/ParameterSet.h"
23 
24 #include <iostream>
25 #include <set>
26 #include <cmath>
27 #include <algorithm>
28 #include <boost/math/special_functions.hpp>
29 #include <boost/math/special_functions/gamma.hpp>
30 #include "TNtuple.h"
31 
32 // To convert TDC to NOvA time
33 #include <NovaTimingUtilities/TimingUtilities.h>
34 
35 #include "Geometry/Geometry.h"
39 #include "MEFinder/MEClusters.h"
40 #include "MCCheater/BackTracker.h"
41 #include "RecoBase/CellHit.h"
42 #include "RecoBase/Track.h"
43 #include "RecoBase/Prong.h"
44 #include "RecoBase/Cluster.h"
45 #include "RecoBase/RecoHit.h"
46 #include "RecoBase/Vertex.h"
47 #include "Calibrator/Calibrator.h"
48 #include "Utilities/func/MathUtil.h"
49 #include "Simulation/FLSHit.h"
50 #include "Simulation/FLSHitList.h"
51 #include "Simulation/Particle.h"
55 
57 
58 #include <TH2D.h>
59 #include <iostream>
60 #include <stdlib.h>
61 #include <fstream>
62 
63 namespace moonshadowana {
64  class MoonShadowAna;
65 
66  const double PLANE_WIDTH = 6.6522; // cm
67  const double CELL_WIDTH = 3.9887; // cm
68  const double TDC_TICK = 15.625; // ns
69  const double SPEED_OF_LIGHT = 30.0; // cm / ns
70 
71  const unsigned N_CELLS = 384; // cells
72  const unsigned N_PLANES = 896; // planes
73 }
74 
76 public:
77  explicit MoonShadowAna(fhicl::ParameterSet const & p);
78  virtual ~MoonShadowAna();
79 
80  void analyze(art::Event const & e) override;
81 
82  void beginJob() override;
83 
84  TVector3 vec(double zen, double azi);
85  double GetDot(double zen1, double azi1, double zen2, double azi2) const;
86  TVector3 AnglesToVector(double zen, double azi) const;
87  void InitLookupTable();
88  double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2);
89  bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2);
90  double GetRate(double zen, double azi);
91  double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2);
92  double GetSmartPrescale(double zen, double azi);
93  bool GetSmartPrescaleDecision(double prescale);
94  void GetHisto(std::string filename, int xmax, int ymax, TH2D* hist);
96 
97  void WriteEventInfo(const art::Event& event);
98  void WriteTrackInfo(const rb::Track& track, long timeStart);
99  void WriteHitListInfo(const art::PtrVector<rb::CellHit> hlist, long timeStart, TVector3 start, TVector3 end);
100  int GetX(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end);
101  int GetY(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end);
102  TVector3 GetDetVec(TVector3 vec) const;
103 
104 private:
105  TNtuple* fNtuple;
108  // Lookup table for introducing "smart prescale"
110  int bin_size;
111 
112  // cut on angle
113  const double fDotCut = cos(5*M_PI/180);
114  // Limit on what rate is safe for detector hardware
115  double max_rate;
116  // We will cut more tracks on the border of 10 degree window than in the middle.
118  // file name of lookup table
120 
121  std::ofstream fOutputEvd;
122  bool init;
123 
124 };
125 
127  : EDAnalyzer (p)
128  , fTrackModuleLabel (p.get<std::string>("TrackModuleLabel")),
129  max_rate(p.get<double>("AvgRate")),
130  angular_stdev(p.get<double>("AngularSuppression")),
132 {
133  InitLookupTable();
134  bin_size = 1;
135 
136  init = false;
137 
138  //fOutputEvd.open("OFFEvd_Data.txt");
139 }
140 
142 {
143 fOutputEvd.close();
144 }
145 
147 {
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");
150 }
151 
152  // write event info
154 {
155  if(!init)
156  {
157  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"));
158  init = true;
159  }
160 
161  // Event info
162  fOutputEvd << "E " << event.event() << " ";
163  fOutputEvd << event.time().timeHigh() << std::endl;
164 
165  // Tracks
167  event.getByLabel(fTrackModuleLabel, "", tracks);
168 
169  for(unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx)
170  {
171  const rb::Track& trk = (*tracks)[trkIdx];
172  if(trk.Is2D())continue;
173 
174 
175 
176  WriteTrackInfo(trk, event.time().timeHigh());
177  }
178 }
179 
181 {
182  TVector3 start = GetDetVec(track.Start());
183  TVector3 end = GetDetVec(track.Stop());
184  if(start.Y() < end.Y()) std::swap(start, end); // all tracks are down-going
185 
186  //Check whether 3D track vector makes sense
187  if(
188  (start.X() == 0 && end.X() == 0) ||
189  (start.Y() == 0 && end.Y() == 0) ||
190  (start.Z() == 0 && end.Z() == 0))
191  {
192  return;
193  }
194 
195  // Track info
196  fOutputEvd << "T ";
197  fOutputEvd << "(" << start.X() << "," << start.Y() << "," << start.Z() << ") ";
198  fOutputEvd << "(" << end.X() << "," << end.Y() << "," << end.Z() << ")";
200 
201  art::PtrVector<rb::CellHit> hlist = track.AllCells();
202 
203  WriteHitListInfo(hlist, timeStart, start, end);
204 
205 }
206 void moonshadowana::MoonShadowAna::WriteHitListInfo(const art::PtrVector<rb::CellHit> hlist, long timeStart, TVector3 start, TVector3 end)
207 {
208  for (size_t j=0; j < hlist.size(); ++j) {
209  art::Ptr<rb::CellHit> h = hlist.at(j);
210 
211  if(j != 0)
212  fOutputEvd << ",";
213 
214  int view = (int)h->View();
215  fOutputEvd << ((long)h->TDC())/64 << "/";
216  fOutputEvd << (view == 0 ? (int)h->Cell() : GetX(h, start, end)) << "/";
217  fOutputEvd << (view == 1 ? (int)h->Cell() : GetY(h, start, end)) << "/";
218  fOutputEvd << (int)h->Plane() << "/";
219  fOutputEvd << view + 1;
220  }
221 
223 }
225 {
226  int z = h->Plane();
227 
228  double xSpread = end.X() - start.X();
229  double zSpread = end.Z() - start.Z();
230 
231  double zSpreadLoc = z - start.Z();
232 
233  double x = start.X() + xSpread * zSpreadLoc/zSpread;
234 
235  return round(x);
236 }
238 {
239  int z = h->Plane();
240 
241  double ySpread = end.Y() - start.Y();
242  double zSpread = end.Z() - start.Z();
243 
244  double zSpreadLoc = z - start.Z();
245 
246  double y = start.Y() + ySpread * zSpreadLoc/zSpread;
247 
248  return round(y);
249 }
250 
252 {
253  int x = int((vec.X() + CELL_WIDTH * N_CELLS /2) / CELL_WIDTH);
254  int y = int((vec.Y() + CELL_WIDTH * N_CELLS /2) / CELL_WIDTH);
255  int z = int(vec.Z() / PLANE_WIDTH);
256 
257  return TVector3(x,y,z);
258 }
259 
261 
262  if(e.event() == 486130)
263  WriteEventInfo(e);
264 
265  //fOutputEvd << "E " << e.event() << std::endl;
266  art::Timestamp ts = e.time();
267  time_t timeSec = ts.timeHigh();
268  // std::cout<<"timeSec "<<timeSec<<std::endl;
269 
270  double zen_moon, azi_moon, zen_sun, azi_sun;
271  loc->GetMoonPosition_FD(timeSec,zen_moon, azi_moon);
272  loc->GetSunPosition_FD(timeSec,zen_sun, azi_sun);
273 
275  e.getByLabel(fTrackModuleLabel, "", tracks);
276 
277  // float vars[10];
278  float vars[9];
279  vars[0] = zen_moon; // zen_moon
280  vars[1] = azi_moon; // azi_moon
281 
282  double smartPrescale = GetSmartPrescale(zen_moon, azi_moon, zen_sun, azi_sun);
283 
284  double zen_trk, azi_trk, cos_trk_moon, cos_trk_sun;
285  for(unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx){
286  const rb::Track& trk = (*tracks)[trkIdx];
287  if(trk.Is2D())continue;
288 
289  //const TVector3 dir = (trk.Start().Y() > trk.Stop().Y()) ? -trk.Dir() : trk.StopDir();
290 
291  TVector3 start = trk.Start();
292  TVector3 end = trk.Stop();
293  if(start.Y() < end.Y()) std::swap(start, end);
294  const TVector3 dir = (start-end).Unit();
295 
296  // test print out params of the track
297  //std::cout << "Track #" << trkIdx << std::endl;
298  //std::cout << start.X() << "," << start.Y() << "," << start.Z() << " || ";
299  //std::cout << end.X() << "," << end.Y() << "," << end.Z() << std::endl;
300 
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;
305 
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));
308  (void)cos_trk_sun; // Suppress unused variable warning
309 
310  vars[2] = zen_trk; // zen_trk
311  vars[3] = azi_trk; // azi_trk
312  vars[4] = trk.NCell(); // nhit
313  vars[5] = trk.TotalLength(); // len
314  vars[6] = cos_trk_moon; // cosine angle trk and moon
315 
316  double angularPrescale = GetAngularPrescale(zen_moon, azi_moon, zen_trk, azi_trk);
317 
318  vars[7] = smartPrescale;
319  vars[8] = angularPrescale;
320 
321  /*if(true)//cos_trk_moon > cos(5*3.14159/180))
322  {
323  fOutputEvd << "T " << "(" << start.X() << "," << start.Y() << "," << start.Z() << ") ";
324  fOutputEvd << "(" << end.X() << "," << end.Y() << "," << end.Z() << ") ";
325  fOutputEvd << std::endl;
326  }*/
327 
328  fNtuple->Fill(vars);
329  }
330 
331 }
332 
333 TVector3 moonshadowana::MoonShadowAna::vec(double zen, double azi){
334  zen *= M_PI / 180;
335  azi *= M_PI / 180;
336  return TVector3(sin(zen)*cos(azi), sin(zen)*sin(azi), cos(zen));
337 }
338 
339  double moonshadowana::MoonShadowAna::GetRate(double zen, double azi)
340  {
341  if(zen > 90)
342  return 0;
343  // Smart prescale is ratio of maximal rate to expected rate divided by angular prescale factor for given position
344  int bin_azi = round(azi/bin_size);
345  int bin_zen = round(zen/bin_size);
346 
347  //return 0;
348  return rate_collected_angular->GetBinContent(bin_azi, bin_zen);
349  }
350 
351  double moonshadowana::MoonShadowAna::GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)
352  {
353  double dTheta = acos(GetDot(zen1, azi1, zen2, azi2)) * 180 / M_PI;
354  return exp(-dTheta*dTheta/(angular_stdev*angular_stdev));
355  }
356  bool moonshadowana::MoonShadowAna::GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
357  {
358  double prescale = GetAngularPrescale(zen1, azi1, zen2, azi2);
359 
360  double r = (double)rand() / (double)RAND_MAX;
361 
362  return r < prescale;
363  }
364 
366  {
367  // Init lookup table and make 2d hist with expected rate
368  // Get histogram with full rates in 2degx2deg bins
369  int azi_max = 360;
370  int zen_max = 90;
371 
372  rate_collected_angular = new TH2D("LookupPrescale", "", azi_max, 0, 360, zen_max, 0, 90);
373 
374  GetHisto(fLookupFile, azi_max, zen_max, rate_collected_angular);
375  }
376 
378 {
379  std::ifstream data;
380  data.open(filename.c_str());
381  for(int j=0; j<ymax;j++)
382  {
383  for(int i=0; i<xmax; i++)
384  {
385  double val;
386  data >> val;
387  hist->SetBinContent(i,j, val);
388  }
389  }
391  TH2D* htmp = tfs->make<TH2D>("LookupPrescale_stored", "", 360, 0, 360, 90, 0, 90);
392  htmp = (TH2D*)hist->Clone();
393  htmp->Write();
394  data.close();
395 }
396 
397 double moonshadowana::MoonShadowAna::GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
398  {
399  double rate1 = GetRate(zen1, azi1);
400  double rate2 = GetRate(zen2, azi2);
401 
402  if(rate1 == 0 && rate2 == 0)
403  return 0;
404 
405  return max_rate / (rate1 + rate2);
406  }
407  double moonshadowana::MoonShadowAna::GetSmartPrescale(double zen, double azi)
408  {
409  double rate = GetRate(zen, azi);
410 
411  if(rate == 0)
412  return 0;
413 
414  return max_rate / rate;
415  }
416 
418  {
419  // Is it good enough for small prescale factors ~1/100?
420  double r = (double)rand() / (double)RAND_MAX;
421  return r < prescale;
422  }
423 
424 double moonshadowana::MoonShadowAna::GetDot(double zen1, double azi1, double zen2, double azi2) const
425  {
426  if(zen1 == zen2 && azi1 == azi2)
427  return 1;
428 
429  zen1 *= M_PI / 180;
430  zen2 *= M_PI / 180;
431  azi1 *= M_PI / 180;
432  azi2 *= M_PI / 180;
433  double dot = sin(zen1)*cos(azi1)*sin(zen2)*cos(azi2) + sin(zen1)*sin(azi1)*sin(zen2)*sin(azi2) + cos(zen1)*cos(zen2);
434 
435  return dot;
436  }
437 
438 // find lookup table file
440 {
441  cet::search_path sp("FW_SEARCH_PATH");
442  std::string lookupFname;
443  sp.find_file(("Eval/PrescaleLookup.txt"), lookupFname);
444  assert(!lookupFname.empty());
445  return lookupFname;
446 }
447 
SubRunNumber_t subRun() const
Definition: Event.h:72
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.
Definition: Cluster.cxx:134
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.
Definition: RawDigit.h:94
const double CELL_WIDTH
const double TDC_TICK
bool Is2D() const
Definition: Cluster.h:96
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
void GetMoonPosition_FD(time_t time, double &zen, double &azi)
const char * p
Definition: xmltok.h:285
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
Definition: event.h:19
bool GetSmartPrescaleDecision(double prescale)
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
void WriteEventInfo(const art::Event &event)
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_t ymax
Definition: plot.C:25
Locate positions of celestial bodies.
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
const unsigned N_CELLS
double GetDot(double zen1, double azi1, double zen2, double azi2) const
const double SPEED_OF_LIGHT
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
MoonShadowAna(fhicl::ParameterSet const &p)
art::ServiceHandle< locator::CelestialLocator > loc
const std::map< std::pair< std::string, std::string >, Variable > vars
const double PLANE_WIDTH
const double j
Definition: BetheBloch.cxx:29
EventNumber_t event() const
Definition: Event.h:67
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)
TVector3 AnglesToVector(double zen, double azi) const
TVector3 Unit() const
Vertex location in position and time.
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
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
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
T sin(T number)
Definition: d0nt_math.hpp:132
bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
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 vec(double zen, double azi)
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
Timestamp time() const
Definition: Event.h:61
TVector3 GetDetVec(TVector3 vec) const
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
void WriteTrackInfo(const rb::Track &track, long timeStart)
const unsigned N_PLANES
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)
T atan2(T number)
Definition: d0nt_math.hpp:72
double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)