DetectorRateShutOff_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file DetectorRateShutOff_module.cc
3 /// \module analyzer
4 /// \brief Extract, format and save all relevant information for studies
5 /// related to rates of particles passing through the detector
6 /// and the associated effects on the detector, e.g. FEB shut-offs.
7 /// \author Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
8 /// \date August 2020
9 ////////////////////////////////////////////////////////////////////////////
10 
11 // framework
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // nova
26 #include "DAQDataFormats/NanoSliceVersionConvention.h"
27 #include "Geometry/Geometry.h"
29 #include "NovaDAQConventions/DAQConventions.h"
30 #include "RawData/DAQHeader.h"
31 #include "RawData/RawDigit.h"
32 #include "RawData/RawTrigger.h"
33 #include "RecoBase/HoughResult.h"
34 #include "DAQChannelMap/DAQChannelMap.h"
35 #include "DAQChannelMap/HardwareDisplay.h"
37 #include "NovaTimingUtilities/TimingUtilities.h"
39 //#include "SummaryData/TBSpillCounters.h"
40 
41 // c++
42 #include <vector>
43 #include <set>
44 
45 // ROOT includes
46 #include "TH1I.h"
47 #include "TH1D.h"
48 #include "TGraph.h"
49 #include "TMultiGraph.h"
50 #include "TProfile.h"
51 #include "TCanvas.h"
52 #include "TMarker.h"
53 #include "TFile.h"
54 
55 // -----------------------------------------------------------------------
56 namespace tbana {
57  class DetectorRateShutOff;
58 }
59 
60 // -----------------------------------------------------------------------
62 
63 public:
64 
65  explicit DetectorRateShutOff(const fhicl::ParameterSet& pset);
66  // virtual ~DetectorRateShutOff();
67 
68  void reconfigure(const fhicl::ParameterSet& pset);
69  void analyze(const art::Event& evt);
70  void beginJob();
71  void endJob();
72 
73  void StartSpill();
74  void EndSpill();
75 
76 private:
77 
78  // methods
79 
80  // labels
85 
86  // ana
89  float fTimeBin;
90 
91  double fMC6IC;
92  double fMC6Col;
93 
94  std::vector<std::vector<unsigned long long> > fHitTimes;
95  std::vector<std::vector<unsigned long long> > fShutOffTimes;
96  std::vector<double> fMuonCounters;
97  std::vector<std::pair<double, double> > fHoughX;
98  std::vector<std::pair<double, double> > fHoughY;
99 
100  // data
101 
102  // general
103  TH1F* hMC6IC;
104  TH1F* hMC6Col;
105 
106  // shut-off proxy
112 
113  // hough rates
114  TH1D* hHoughRate;
129 
130  // shut-offs mc6ic
131  TGraph* hShutOffsMC6IC;
141 
142  // counter shut-offs
148 
149  // hough shut-offs
170 
171  // profiles
172  std::map<unsigned int, TH1I*> hFEBSpillProfile;
173 
174  // benford
175  TMultiGraph* hBenfordLaw;
176  TMultiGraph* hBenfordLawNoShutOff;
177  TMultiGraph* hBenfordLawShutOff;
180 
181  // services
185 
186 };
187 
188 // -----------------------------------------------------------------------
190 
191 // -----------------------------------------------------------------------
193  this->reconfigure(pset);
194 }
195 
196 // // -----------------------------------------------------------------------
197 // tbana::DetectorRateShutOff::~DetectorRateShutOff() {
198 // }
199 
200 // -----------------------------------------------------------------------
202  fRawDataLabel = pset.get<std::string>("RawDataLabel");
203  fSpillInfoLabel = pset.get<std::string>("SpillInfoLabel");
204  fSpillCountersLabel = pset.get<std::string>("SpillCountersLabel");
205  fHoughModuleLabel = pset.get<std::string>("HoughModuleLabel");
206  fTimeBin = pset.get<float>("TimeBin"); // ms
207 }
208 
209 // -----------------------------------------------------------------------
212  fSpillNumber = -1;
213 
214  // general
215  hMC6IC = fFileService->make<TH1F>("MC6IC", ";MC6IC;", 100, 1e8, 1e11);
216  hMC6Col = fFileService->make<TH1F>("MC6Col", ";MC6 Collimator Aperture (mm);", 30, 0, 30);
217 
218  // shut-off proxy
220  ("ShutOffProxyDiff", ";Shut-Off - Proxy Shut-Off;", 20, -10, 10);
221  hShutOffProxyDiff->GetXaxis()->SetNdivisions(20);
222  hShutOffProxyDiff->GetXaxis()->CenterLabels();
224  ("WestShutOffProxyDiff", ";Shut-Off - Proxy Shut-Off;", 20, -10, 10);
225  hWestShutOffProxyDiff->GetXaxis()->SetNdivisions(20);
226  hWestShutOffProxyDiff->GetXaxis()->CenterLabels();
228  ("EastShutOffProxyDiff", ";Shut-Off - Proxy Shut-Off;", 20, -10, 10);
229  hEastShutOffProxyDiff->GetXaxis()->SetNdivisions(20);
230  hEastShutOffProxyDiff->GetXaxis()->CenterLabels();
232  ("TopShutOffProxyDiff", ";Shut-Off - Proxy Shut-Off;", 20, -10, 10);
233  hTopShutOffProxyDiff->GetXaxis()->SetNdivisions(20);
234  hTopShutOffProxyDiff->GetXaxis()->CenterLabels();
236  ("BottomShutOffProxyDiff", ";Shut-Off - Proxy Shut-Off;", 20, -10, 10);
237  hBottomShutOffProxyDiff->GetXaxis()->SetNdivisions(20);
238  hBottomShutOffProxyDiff->GetXaxis()->CenterLabels();
239 
240  // hough rates
241  hHoughRate = fFileService->make<TH1D>("HoughRate",";Hough Rate;", 100, 0, 150000);
242  hEastHoughRate = fFileService->make<TH1D>("EastHoughRate",";Hough Rate (East);", 100, 0, 70000);
243  hWestHoughRate = fFileService->make<TH1D>("WestHoughRate",";Hough Rate (West);", 100, 0, 150000);
244  hTopHoughRate = fFileService->make<TH1D>("TopHoughRate",";Hough Rate (Top);", 100, 0, 150000);
245  hBottomHoughRate = fFileService->make<TH1D>("BottomHoughRate",";Hough Rate (Bottom);", 100, 0, 100000);
246  hHoughRateMC6IC = fFileService->makeAndRegister<TGraph>("HoughRateMC6IC",";MC6IC;Hough Rate;");
247  hEastHoughRateMC6IC = fFileService->makeAndRegister<TGraph>("EastHoughRateMC6IC",";MC6IC;Hough Rate (East);");
248  hWestHoughRateMC6IC = fFileService->makeAndRegister<TGraph>("WestHoughRateMC6IC",";MC6IC;Hough Rate (West);");
249  hTopHoughRateMC6IC = fFileService->makeAndRegister<TGraph>("TopHoughRateMC6IC",";MC6IC;Hough Rate (Top);");
250  hBottomHoughRateMC6IC = fFileService->makeAndRegister<TGraph>("BottomHoughRateMC6IC",";MC6IC;Hough Rate (Bottom);");
251  hHoughRateMC6Col = fFileService->makeAndRegister<TGraph>("HoughRateMC6Col",";MC6 Collimator Aperature (mm);Hough Rate;");
252  hEastHoughRateMC6Col = fFileService->makeAndRegister<TGraph>("EastHoughRateMC6Col",";MC6 Collimator Aperature (mm);Hough Rate (East);");
253  hWestHoughRateMC6Col = fFileService->makeAndRegister<TGraph>("WestHoughRateMC6Col",";MC6 Collimator Aperature (mm);Hough Rate (West);");
254  hTopHoughRateMC6Col = fFileService->makeAndRegister<TGraph>("TopHoughRateMC6Col",";MC6 Collimator Aperature (mm);Hough Rate (Top);");
255  hBottomHoughRateMC6Col = fFileService->makeAndRegister<TGraph>("BottomHoughRateMC6Col",";MC6 Collimator Aperature (mm);Hough Rate (Bottom);");
256 
257  // shut-offs mc6ic
259  ("ShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
261  ("WestShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
263  ("EastShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
265  ("TopShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
267  ("BottomShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
269  ("ProxyShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
271  ("WestProxyShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
273  ("EastProxyShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
275  ("TopProxyShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
277  ("BottomProxyShutOffsMC6IC", ";MC6IC;Number of Shut-Offs;");
278 
279  // counter shut-offs
281  ("ShutOffsMuonCounters", ";Counts on Front Face (Multiple of East Side Counters);Number of Shut-Offs;");
283  ("WestCountersShutOffs", ";Counts on Front Face (West) (Multiple of East Side Counters);Number of Shut-Offs;");
285  ("EastCountersShutOffs", ";Counts on Front Face (East) (Multiple of East Side Counters);Number of Shut-Offs;");
287  ("TopCountersShutOffs", ";Counts on Front Face (Top) (Multiple of East Side Counters);Number of Shut-Offs;");
289  ("BottomCountersShutOffs", ";Counts on Front Face (Bottom) (Multiple of East Side Counters);Number of Shut-Offs;");
290 
291  // hough shut-offs
292  // shut-offs
294  ("ShutOffsHoughTracks", ";Hough Tracks;Number of Shut-Offs;");
296  ("WestHoughsShutOffs", ";Hough Tracks (West);Number of Shut-Offs;");
298  ("EastHoughsShutOffs", ";Hough Tracks (East);Number of Shut-Offs;");
300  ("TopHoughsShutOffs", ";Hough Tracks (Top);Number of Shut-Offs;");
302  ("BottomHoughsShutOffs", ";Hough Tracks (Bottom);Number of Shut-Offs;");
304  ("ShutOffsRelHoughTracks", ";Hough Tracks (Multiple of East Quadrant);Number of Shut-Offs;");
306  ("WestRelHoughsShutOffs", ";Hough Tracks (West) (Multiple of East Quadrant);Number of Shut-Offs;");
308  ("EastRelHoughsShutOffs", ";Hough Tracks (East) (Multiple of East Quadrant);Number of Shut-Offs;");
310  ("TopRelHoughsShutOffs", ";Hough Tracks (Top) (Multiple of East Quadrant);Number of Shut-Offs;");
312  ("BottomRelHoughsShutOffs", ";Hough Tracks (Bottom) (Multiple of East Quadrant);Number of Shut-Offs;");
313  // proxy
315  ("ProxyShutOffsHoughTracks", ";Hough Tracks;Number of Shut-Offs;");
317  ("WestHoughsProxyShutOffs", ";Hough Tracks (West);Number of Shut-Offs;");
319  ("EastHoughsProxyShutOffs", ";Hough Tracks (East);Number of Shut-Offs;");
321  ("TopHoughsProxyShutOffs", ";Hough Tracks (Top);Number of Shut-Offs;");
323  ("BottomHoughsProxyShutOffs", ";Hough Tracks (Bottom);Number of Shut-Offs;");
325  ("ProxyShutOffsRelHoughTracks", ";Hough Tracks (Multiple of East Quadrant);Number of Shut-Offs;");
327  ("WestRelHoughsProxyShutOffs", ";Hough Tracks (West) (Multiple of East Quadrant);Number of Shut-Offs;");
329  ("EastRelHoughsProxyShutOffs", ";Hough Tracks (East) (Multiple of East Quadrant);Number of Shut-Offs;");
331  ("TopRelHoughsProxyShutOffs", ";Hough Tracks (Top) (Multiple of East Quadrant);Number of Shut-Offs;");
333  ("BottomRelHoughsProxyShutOffs", ";Hough Tracks (Bottom) (Multiple of East Quadrant);Number of Shut-Offs;");
334 
335  // profiles
336  for (const unsigned int& feb : {0,16,32,48})
338  (Form("FEB%dSpillProfile", feb), Form(";Time in Spill (ns);Track Rate in FEB %2d (%f ms);", feb, fTimeBin), 5e9/(fTimeBin*1.e6), 0, 5e9);
339 
340  // benford
341  hBenfordLaw = fFileService->makeAndRegister<TMultiGraph>("BenfordLaw",";First Digit in Time Between Hits (ns);Frequency;");
342  hBenfordLawNoShutOff = fFileService->makeAndRegister<TMultiGraph>("BenfordLawNoShutOff",";First Digit in Time Between Hits (ns);Frequency;");
343  hBenfordLawShutOff = fFileService->makeAndRegister<TMultiGraph>("BenfordLawShutOff",";First Digit in Time Between Hits (ns);Frequency;");
344  hBenfordLawExpoConstNoShutOff = fFileService->make<TH1F>("BenfordLawExpoConstNoShutOff",";Exponential Constant;",100,7,10);
345  hBenfordLawExpoConstShutOff = fFileService->make<TH1F>("BenfordLawExpoConstShutOff",";Exponential Constant;",100,7,10);
346 }
347 
348 // -----------------------------------------------------------------------
350  //this->EndSpill();
351  // TCanvas* canv = new TCanvas("canv2", "", 800, 600);
352  // hBenfordLaw->Draw("AP");
353  // canv->SaveAs("BenfordLaw.png");
354  // delete canv;
355 
356  return;
357 }
358 
359 // -----------------------------------------------------------------------
361 
362  // Get ifdb info
363  art::Handle<sumdata::MCenterData> spillInfoHandle;
364  if (evt.getByLabel(fSpillInfoLabel, spillInfoHandle)) {
365  fMC6IC = spillInfoHandle->mc6int;
366  fMC6Col = spillInfoHandle->mc6col;
367  }
368  // art::Handle<sumdata::TBSpillCounters> spillCountersHandle;
369  // sumdata::TBSpillCounters spillCounters;
370  // if (evt.getByLabel(fSpillCountersLabel, spillCountersHandle))
371  // spillCounters = *spillCountersHandle;
372 
373  // Get hits
375  std::vector<art::Ptr<rawdata::RawDigit> > digits;
376  if (evt.getByLabel(fRawDataLabel, digitHandle))
377  art::fill_ptr_vector(digits, digitHandle);
378 
379  // Get trigger information
381  std::vector<art::Ptr<rawdata::RawTrigger> > triggers;
382  if (evt.getByLabel(fRawDataLabel, triggerHandle))
383  art::fill_ptr_vector(triggers, triggerHandle);
384  const art::Ptr<rawdata::RawTrigger> trigger = triggers.at(0);
385 
386  // Get Hough tracks
388  std::vector<art::Ptr<rb::HoughResult> > houghs;
389  if (evt.getByLabel(fHoughModuleLabel, houghHandle))
390  art::fill_ptr_vector(houghs, houghHandle);
391 
392  long long int trigger_master = trigger->fTriggerHeader_MasterTriggerNumber;
393  unsigned long long trigger_number = trigger->fTriggerHeader_TriggerNumber;
394  unsigned long long trigger_seq = trigger_number - trigger_master;
395  uint32_t trigger_len = trigger->fTriggerRange_TriggerLength;
396  // unsigned long long trigger_time = trigger->fTriggerTimingMarker_TimeStart; // subtrigger time in NOvA
397 
398  // Time of this subtrigger wrt the start of the spill
399  unsigned long long subtrigger_time = trigger_len * trigger_seq * 500; // ns
400 
401  // New spill
402  if (trigger_master != fPreviousTriggerMaster) {
403  if (fSpillNumber >= 0)
404  this->EndSpill();
405  ++fSpillNumber;
406  fPreviousTriggerMaster = trigger_master;
407  this->StartSpill();
408  }
409 
410  // // Muon counters
411  // if (fMuonCounters[0] < 0 and spillCounters.detface1 > 0)
412  // fMuonCounters = {spillCounters.detface1, spillCounters.detface2, spillCounters.detface3, spillCounters.detface4};
413 
414  // Look at hits
415  for (std::vector<art::Ptr<rawdata::RawDigit> >::const_iterator digitIt = digits.begin();
416  digitIt != digits.end(); ++digitIt) {
417 
418  // geometry
420  double pos[3], dpos[3];
421  uint plane = fChannelMap->GetPlane(digitIt->get());
422  uint cell = fChannelMap->GetCell(digitIt->get());
423  fGeo->CellInfo(plane, cell, &view, pos, dpos);
424 
425  // channel mapping
426  uint32_t dchan = (*digitIt)->DaqChannel();
427  uint8_t status = (*digitIt)->fFEBStatus;
428  int feb = fChannelMap->getFEB(dchan);
429  int dcm = fChannelMap->getDCM(dchan);
430  int feb_det = feb+(64*(dcm-1));
431 
432  // time
433  int32_t hit_time = (*digitIt)->TDC() * 15.625; // ns
434  unsigned long long hit_time_spill = hit_time + subtrigger_time; // ns
435 
436  fHitTimes.at(feb_det).push_back(hit_time_spill);
437  if (status == 32)
438  fShutOffTimes.at(feb_det).push_back(hit_time_spill);
439 
440  }
441 
442  // Look at Hough lines
443  for (std::vector<art::Ptr<rb::HoughResult> >::const_iterator houghIt = houghs.begin();
444  houghIt != houghs.end(); ++houghIt) {
445 
446  geo::View_t view = (*houghIt)->fView;
447  double time = (*houghIt)->fTNSlo + subtrigger_time; // ns
448 
449  double Z0 = (*houghIt)->fZ0;
450  double XY0 = (*houghIt)->fXY0;
451 
452  // get the peaks
453  const std::vector<rb::HoughPeak> peaks = (*houghIt)->fPeak;
454  for (std::vector<rb::HoughPeak>::const_iterator peakIt = peaks.begin(); peakIt != peaks.end(); ++peakIt) {
455 
456  double rho = peakIt->fRho;
457  double theta = peakIt->fTheta;
458 
459  // calculate slope and intercept
460  double c = cos(theta);
461  double s = sin(theta);
462  if(s == 0) s = 1.0E-9;
463  double m = -c/s;
464  double b = rho/s - m*Z0 + XY0;
465 
466  if (view == geo::kX)
467  fHoughX.push_back(std::make_pair(time, b));
468  if(view == geo::kY)
469  fHoughY.push_back(std::make_pair(time, b));
470 
471  }
472  }
473 
474  return;
475 
476 }
477 
478 // -----------------------------------------------------------------------
480 
481  fHitTimes.clear();
482  fHitTimes.resize(136);
483  fShutOffTimes.clear();
484  fShutOffTimes.resize(136);
485  fMuonCounters.clear();
486  fMuonCounters.resize(4, -1);
487  fHoughX.clear();
488  fHoughY.clear();
489 
490  return;
491 }
492 
493 // -----------------------------------------------------------------------
495 
496  std::cout << "End of spill" << std::endl;
497 
498  hMC6IC->Fill(fMC6IC);
499  hMC6Col->Fill(fMC6Col);
500 
501  // determine a rough number of total FEB shut-offs (combining hits from multiple cells)
502  std::vector<std::vector<unsigned long long> > shutOffs(136);
503  for (unsigned int feb = 0; feb < 136; ++feb) {
504  std::sort(fShutOffTimes.at(feb).begin(), fShutOffTimes.at(feb).end());
505  for (std::vector<unsigned long long>::const_iterator hittimeIt = fShutOffTimes.at(feb).begin();
506  hittimeIt != fShutOffTimes.at(feb).end(); ++hittimeIt) {
507  bool newshutoff = true;
508  for (std::vector<unsigned long long>::const_iterator prevShutOffTimeIt = shutOffs.at(feb).begin();
509  prevShutOffTimeIt != shutOffs.at(feb).end(); ++prevShutOffTimeIt)
510  if (*hittimeIt-*prevShutOffTimeIt < 1e6) // 1ms
511  newshutoff = false;
512  if (newshutoff)
513  shutOffs.at(feb).push_back(*hittimeIt);
514  }
515  std::cout << "FEB " << feb << " had " << shutOffs.at(feb).size() << " shutoffs" << std::endl;
516  }
517 
518  // sum up all shut-offs in each FEB region
519  std::vector<unsigned int> west_febs = {16};
520  std::vector<unsigned int> east_febs = {0};
521  std::vector<unsigned int> top_febs = {32};
522  std::vector<unsigned int> bottom_febs = {48};
523  unsigned int west_shutoffs = 0, east_shutoffs = 0, top_shutoffs = 0, bottom_shutoffs = 0;
524  for (std::vector<unsigned int>::const_iterator feb = west_febs.begin(); feb != west_febs.end(); ++feb)
525  west_shutoffs += shutOffs.at(*feb).size();
526  for (std::vector<unsigned int>::const_iterator feb = east_febs.begin(); feb != east_febs.end(); ++feb)
527  east_shutoffs += shutOffs.at(*feb).size();
528  for (std::vector<unsigned int>::const_iterator feb = top_febs.begin(); feb != top_febs.end(); ++feb)
529  top_shutoffs += shutOffs.at(*feb).size();
530  for (std::vector<unsigned int>::const_iterator feb = bottom_febs.begin(); feb != bottom_febs.end(); ++feb)
531  bottom_shutoffs += shutOffs.at(*feb).size();
532 
533  // shut-offs vs MC6IC
534  hShutOffsMC6IC->SetPoint(hShutOffsMC6IC->GetN(), fMC6IC, west_shutoffs);
535  hShutOffsMC6IC->SetPoint(hShutOffsMC6IC->GetN(), fMC6IC, east_shutoffs);
536  hShutOffsMC6IC->SetPoint(hShutOffsMC6IC->GetN(), fMC6IC, top_shutoffs);
537  hShutOffsMC6IC->SetPoint(hShutOffsMC6IC->GetN(), fMC6IC, bottom_shutoffs);
538  hWestShutOffsMC6IC->SetPoint(hWestShutOffsMC6IC->GetN(), fMC6IC, west_shutoffs);
539  hEastShutOffsMC6IC->SetPoint(hEastShutOffsMC6IC->GetN(), fMC6IC, east_shutoffs);
540  hTopShutOffsMC6IC->SetPoint(hTopShutOffsMC6IC->GetN(), fMC6IC, top_shutoffs);
541  hBottomShutOffsMC6IC->SetPoint(hBottomShutOffsMC6IC->GetN(), fMC6IC, bottom_shutoffs);
542 
543  // get the Hough intercepts at the front face
544  std::vector<double> west_hough_times, east_hough_times, top_hough_times, bottom_hough_times;
545  for (std::vector<std::pair<double, double> >::const_iterator houghIt = fHoughX.begin();
546  houghIt != fHoughX.end(); ++houghIt) {
547  if (houghIt->second > 0)
548  west_hough_times.push_back(houghIt->first);
549  else
550  east_hough_times.push_back(houghIt->first);
551  }
552  for (std::vector<std::pair<double, double> >::const_iterator houghIt = fHoughY.begin();
553  houghIt != fHoughY.end(); ++houghIt) {
554  if (houghIt->second > 0)
555  top_hough_times.push_back(houghIt->first);
556  else
557  bottom_hough_times.push_back(houghIt->first);
558  }
559 
560  // fill spill profiles
561  for (std::vector<double>::const_iterator timeIt = west_hough_times.begin(); timeIt != west_hough_times.end(); ++timeIt)
562  hFEBSpillProfile[16]->Fill(*timeIt);
563  for (std::vector<double>::const_iterator timeIt = east_hough_times.begin(); timeIt != east_hough_times.end(); ++timeIt)
564  hFEBSpillProfile[0]->Fill(*timeIt);
565  for (std::vector<double>::const_iterator timeIt = top_hough_times.begin(); timeIt != top_hough_times.end(); ++timeIt)
566  hFEBSpillProfile[32]->Fill(*timeIt);
567  for (std::vector<double>::const_iterator timeIt = bottom_hough_times.begin(); timeIt != bottom_hough_times.end(); ++timeIt)
568  hFEBSpillProfile[48]->Fill(*timeIt);
569 
570  // draw profiles with shut-offs
571  TCanvas* canv = new TCanvas("canv", "", 800, 600);
572  TLine line;
573  line.SetLineColor(kRed);
574  line.SetLineWidth(2);
575  line.SetLineStyle(kDashed);
576  for (const unsigned int& feb : {0,16,32,48}) {
577  canv->cd();
578  canv->Clear();
579  hFEBSpillProfile[feb]->Draw("hist");
580  canv->Modified(); canv->Update();
581  for (std::vector<unsigned long long>::const_iterator shutOffIt = shutOffs.at(feb).begin();
582  shutOffIt != shutOffs.at(feb).end(); ++shutOffIt)
583  line.DrawLine(*shutOffIt, 0, *shutOffIt, canv->GetUymax());
584  canv->SaveAs(Form("figs/Spill%dFEB%02dSpillProfile.png", fSpillNumber, feb));
585  canv->cd();
586  canv->Clear();
587  TH1* scaled = (TH1*)hFEBSpillProfile[feb]->Clone(Form("%sScaled", hFEBSpillProfile[feb]->GetName()));
588  scaled->Divide(hFEBSpillProfile[0]);
589  scaled->GetYaxis()->SetTitle(Form("Hits on FEB %02d / Hits on FEB 00", feb));
590  scaled->Draw("hist");
591  canv->Modified(); canv->Update();
592  for (std::vector<unsigned long long>::const_iterator shutOffIt = shutOffs.at(feb).begin();
593  shutOffIt != shutOffs.at(feb).end(); ++shutOffIt)
594  line.DrawLine(*shutOffIt, 0, *shutOffIt, canv->GetUymax());
595  canv->SaveAs(Form("figs/Spill%dFEB%02dSpillProfileNorm.png", fSpillNumber, feb));
596  delete scaled;
597  }
598  delete canv;
599 
600  float west_houghs = west_hough_times.size();
601  float east_houghs = east_hough_times.size();
602  float top_houghs = top_hough_times.size();
603  float bottom_houghs = bottom_hough_times.size();
604 
605  // hough rates
606  hHoughRate->Fill(east_houghs);
607  hHoughRate->Fill(west_houghs);
608  hHoughRate->Fill(top_houghs);
609  hHoughRate->Fill(bottom_houghs);
610  hEastHoughRate->Fill(east_houghs);
611  hWestHoughRate->Fill(west_houghs);
612  hTopHoughRate->Fill(top_houghs);
613  hBottomHoughRate->Fill(bottom_houghs);
614 
615  // hough rates mc6ic
616  hHoughRateMC6IC->SetPoint(hHoughRateMC6IC->GetN(), fMC6IC, east_houghs);
617  hHoughRateMC6IC->SetPoint(hHoughRateMC6IC->GetN(), fMC6IC, west_houghs);
618  hHoughRateMC6IC->SetPoint(hHoughRateMC6IC->GetN(), fMC6IC, top_houghs);
619  hHoughRateMC6IC->SetPoint(hHoughRateMC6IC->GetN(), fMC6IC, bottom_houghs);
620  hEastHoughRateMC6IC->SetPoint(hEastHoughRateMC6IC->GetN(), fMC6IC, east_houghs);
621  hWestHoughRateMC6IC->SetPoint(hEastHoughRateMC6IC->GetN(), fMC6IC, west_houghs);
622  hTopHoughRateMC6IC->SetPoint(hTopHoughRateMC6IC->GetN(), fMC6IC, top_houghs);
623  hBottomHoughRateMC6IC->SetPoint(hBottomHoughRateMC6IC->GetN(), fMC6IC, bottom_houghs);
624 
625  // hough rates mc6col
626  hHoughRateMC6Col->SetPoint(hHoughRateMC6Col->GetN(), fMC6Col, east_houghs);
627  hHoughRateMC6Col->SetPoint(hHoughRateMC6Col->GetN(), fMC6Col, west_houghs);
628  hHoughRateMC6Col->SetPoint(hHoughRateMC6Col->GetN(), fMC6Col, top_houghs);
629  hHoughRateMC6Col->SetPoint(hHoughRateMC6Col->GetN(), fMC6Col, bottom_houghs);
630  hEastHoughRateMC6Col->SetPoint(hEastHoughRateMC6Col->GetN(), fMC6Col, east_houghs);
631  hWestHoughRateMC6Col->SetPoint(hEastHoughRateMC6Col->GetN(), fMC6Col, west_houghs);
632  hTopHoughRateMC6Col->SetPoint(hTopHoughRateMC6Col->GetN(), fMC6Col, top_houghs);
633  hBottomHoughRateMC6Col->SetPoint(hBottomHoughRateMC6Col->GetN(), fMC6Col, bottom_houghs);
634 
635  // hough rates shut-offs
636  hShutOffsHoughTracks->SetPoint(hShutOffsHoughTracks->GetN(), west_houghs, west_shutoffs);
637  hShutOffsHoughTracks->SetPoint(hShutOffsHoughTracks->GetN(), east_houghs, east_shutoffs);
638  hShutOffsHoughTracks->SetPoint(hShutOffsHoughTracks->GetN(), top_houghs, top_shutoffs);
639  hShutOffsHoughTracks->SetPoint(hShutOffsHoughTracks->GetN(), bottom_houghs, bottom_shutoffs);
640  hWestHoughsShutOffs->SetPoint(hWestHoughsShutOffs->GetN(), west_houghs, west_shutoffs);
641  hEastHoughsShutOffs->SetPoint(hEastHoughsShutOffs->GetN(), east_houghs, east_shutoffs);
642  hTopHoughsShutOffs->SetPoint(hTopHoughsShutOffs->GetN(), top_houghs, top_shutoffs);
643  hBottomHoughsShutOffs->SetPoint(hBottomHoughsShutOffs->GetN(), bottom_houghs, bottom_shutoffs);
644 
645  float west_houghs_sc = west_houghs/east_houghs;
646  float east_houghs_sc = east_houghs/east_houghs;
647  float top_houghs_sc = top_houghs/east_houghs;
648  float bottom_houghs_sc = bottom_houghs/east_houghs;
649 
650  hShutOffsRelHoughTracks->SetPoint(hShutOffsRelHoughTracks->GetN(), west_houghs_sc, west_shutoffs);
651  hShutOffsRelHoughTracks->SetPoint(hShutOffsRelHoughTracks->GetN(), east_houghs_sc, east_shutoffs);
652  hShutOffsRelHoughTracks->SetPoint(hShutOffsRelHoughTracks->GetN(), top_houghs_sc, top_shutoffs);
653  hShutOffsRelHoughTracks->SetPoint(hShutOffsRelHoughTracks->GetN(), bottom_houghs_sc, bottom_shutoffs);
654  hWestRelHoughsShutOffs->SetPoint(hWestRelHoughsShutOffs->GetN(), west_houghs_sc, west_shutoffs);
655  hEastRelHoughsShutOffs->SetPoint(hEastRelHoughsShutOffs->GetN(), east_houghs_sc, east_shutoffs);
656  hTopRelHoughsShutOffs->SetPoint(hTopRelHoughsShutOffs->GetN(), top_houghs_sc, top_shutoffs);
657  hBottomRelHoughsShutOffs->SetPoint(hBottomRelHoughsShutOffs->GetN(), bottom_houghs_sc, bottom_shutoffs);
658 
659  // hit-level analysis
660  std::vector<std::vector<unsigned long long> > shutOffsProxy(136);
661  for (unsigned int feb = 0; feb < 136; ++feb) {
662  std::sort(fHitTimes.at(feb).begin(), fHitTimes.at(feb).end());
663  std::vector<unsigned int> firstDigits(10, 0);
664  for (unsigned int hit = 1; hit < fHitTimes.at(feb).size(); ++hit) {
665  unsigned long long diff = fHitTimes.at(feb).at(hit)-fHitTimes.at(feb).at(hit-1);
666  if (diff > 20e6)
667  shutOffsProxy.at(feb).push_back(fHitTimes.at(feb).at(hit-1));
668  while (diff >= 10)
669  diff = diff / 10;
670  ++firstDigits.at(diff);
671  }
672  // Benford's law
673  TGraph* gr = new TGraph(9);
674  for (unsigned int firstDigit = 1; firstDigit <= 9; ++firstDigit)
675  gr->SetPoint(firstDigit-0, firstDigit, firstDigits.at(firstDigit));
676  gr->Fit("expo", "q");
677  if (shutOffs.at(feb).size())
678  hBenfordLawExpoConstShutOff->Fill(gr->GetFunction("expo")->GetParameter(0));
679  else
680  hBenfordLawExpoConstNoShutOff->Fill(gr->GetFunction("expo")->GetParameter(0));
681  gr->SetMarkerStyle(8);
682  if (shutOffs.at(feb).size()) {
683  gr->SetMarkerColor(kRed);
684  gr->SetLineColor(kRed);
685  hBenfordLawShutOff->Add((TGraph*)gr->Clone());
686  }
687  else {
688  gr->SetMarkerColor(kBlack);
689  gr->SetLineColor(kBlack);
690  hBenfordLawNoShutOff->Add((TGraph*)gr->Clone());
691  }
692  hBenfordLaw->Add(gr);
693  gr = nullptr;
694  }
695 
696  // sum up proxy shut-offs in each FEB region
697  unsigned int west_proxy_shutoffs = 0, east_proxy_shutoffs = 0, top_proxy_shutoffs = 0, bottom_proxy_shutoffs = 0;
698  for (std::vector<unsigned int>::const_iterator feb = west_febs.begin(); feb != west_febs.end(); ++feb)
699  west_proxy_shutoffs += shutOffsProxy.at(*feb).size();
700  for (std::vector<unsigned int>::const_iterator feb = east_febs.begin(); feb != east_febs.end(); ++feb)
701  east_proxy_shutoffs += shutOffsProxy.at(*feb).size();
702  for (std::vector<unsigned int>::const_iterator feb = top_febs.begin(); feb != top_febs.end(); ++feb)
703  top_proxy_shutoffs += shutOffsProxy.at(*feb).size();
704  for (std::vector<unsigned int>::const_iterator feb = bottom_febs.begin(); feb != bottom_febs.end(); ++feb)
705  bottom_proxy_shutoffs += shutOffsProxy.at(*feb).size();
706 
707  // difference between real shut-offs
708  hShutOffProxyDiff->Fill(west_shutoffs-west_proxy_shutoffs);
709  hShutOffProxyDiff->Fill(east_shutoffs-east_proxy_shutoffs);
710  hShutOffProxyDiff->Fill(top_shutoffs-top_proxy_shutoffs);
711  hShutOffProxyDiff->Fill(bottom_shutoffs-bottom_proxy_shutoffs);
712  hWestShutOffProxyDiff->Fill(west_shutoffs-west_proxy_shutoffs);
713  hEastShutOffProxyDiff->Fill(east_shutoffs-east_proxy_shutoffs);
714  hTopShutOffProxyDiff->Fill(top_shutoffs-top_proxy_shutoffs);
715  hBottomShutOffProxyDiff->Fill(bottom_shutoffs-bottom_proxy_shutoffs);
716 
717  // shut-offs vs MC6IC
718  hProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, west_proxy_shutoffs);
719  hProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, east_proxy_shutoffs);
720  hProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, top_proxy_shutoffs);
721  hProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, bottom_proxy_shutoffs);
722  hWestProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, west_proxy_shutoffs);
723  hEastProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, east_proxy_shutoffs);
724  hTopProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, top_proxy_shutoffs);
725  hBottomProxyShutOffsMC6IC->SetPoint(hProxyShutOffsMC6IC->GetN(), fMC6IC, bottom_proxy_shutoffs);
726 
727  // shut-offs vs houghs
728  hProxyShutOffsHoughTracks->SetPoint(hProxyShutOffsHoughTracks->GetN(), west_houghs, west_proxy_shutoffs);
729  hProxyShutOffsHoughTracks->SetPoint(hProxyShutOffsHoughTracks->GetN(), east_houghs, east_proxy_shutoffs);
730  hProxyShutOffsHoughTracks->SetPoint(hProxyShutOffsHoughTracks->GetN(), top_houghs, top_proxy_shutoffs);
731  hProxyShutOffsHoughTracks->SetPoint(hProxyShutOffsHoughTracks->GetN(), bottom_houghs, bottom_proxy_shutoffs);
732  hWestHoughsProxyShutOffs->SetPoint(hWestHoughsProxyShutOffs->GetN(), west_houghs, west_proxy_shutoffs);
733  hEastHoughsProxyShutOffs->SetPoint(hEastHoughsProxyShutOffs->GetN(), east_houghs, east_proxy_shutoffs);
734  hTopHoughsProxyShutOffs->SetPoint(hTopHoughsProxyShutOffs->GetN(), top_houghs, top_proxy_shutoffs);
735  hBottomHoughsProxyShutOffs->SetPoint(hBottomHoughsProxyShutOffs->GetN(), bottom_houghs, bottom_proxy_shutoffs);
736 
737  hProxyShutOffsRelHoughTracks->SetPoint(hProxyShutOffsRelHoughTracks->GetN(), west_houghs_sc, west_proxy_shutoffs);
738  hProxyShutOffsRelHoughTracks->SetPoint(hProxyShutOffsRelHoughTracks->GetN(), east_houghs_sc, east_proxy_shutoffs);
739  hProxyShutOffsRelHoughTracks->SetPoint(hProxyShutOffsRelHoughTracks->GetN(), top_houghs_sc, top_proxy_shutoffs);
740  hProxyShutOffsRelHoughTracks->SetPoint(hProxyShutOffsRelHoughTracks->GetN(), bottom_houghs_sc, bottom_proxy_shutoffs);
741  hWestRelHoughsProxyShutOffs->SetPoint(hWestRelHoughsProxyShutOffs->GetN(), west_houghs_sc, west_proxy_shutoffs);
742  hEastRelHoughsProxyShutOffs->SetPoint(hEastRelHoughsProxyShutOffs->GetN(), east_houghs_sc, east_proxy_shutoffs);
743  hTopRelHoughsProxyShutOffs->SetPoint(hTopRelHoughsProxyShutOffs->GetN(), top_houghs_sc, top_proxy_shutoffs);
744  hBottomRelHoughsProxyShutOffs->SetPoint(hBottomRelHoughsProxyShutOffs->GetN(), bottom_houghs_sc, bottom_proxy_shutoffs);
745 
746  // // get the scintillator counts from the front face of the detector
747  // if (fMuonCounters.at(0) < 0)
748  // return;
749  // float west_counters = fMuonCounters.at(0)+fMuonCounters.at(1);
750  // float east_counters = fMuonCounters.at(2)+fMuonCounters.at(3);
751  // float top_counters = fMuonCounters.at(0)+fMuonCounters.at(2);
752  // float bottom_counters = fMuonCounters.at(1)+fMuonCounters.at(3);
753 
754  // float west_counters_sc = west_counters/east_counters;
755  // float east_counters_sc = east_counters/east_counters;
756  // float top_counters_sc = top_counters/east_counters;
757  // float bottom_counters_sc = bottom_counters/east_counters;
758 
759  // hShutOffsMuonCounters->SetPoint(hShutOffsMuonCounters->GetN(), west_counters_sc, west_shutoffs);
760  // hShutOffsMuonCounters->SetPoint(hShutOffsMuonCounters->GetN(), east_counters_sc, east_shutoffs);
761  // hShutOffsMuonCounters->SetPoint(hShutOffsMuonCounters->GetN(), top_counters_sc, top_shutoffs);
762  // hShutOffsMuonCounters->SetPoint(hShutOffsMuonCounters->GetN(), bottom_counters_sc, bottom_shutoffs);
763  // hWestCountersShutOffs->SetPoint(hWestCountersShutOffs->GetN(), west_counters_sc, west_shutoffs);
764  // hEastCountersShutOffs->SetPoint(hEastCountersShutOffs->GetN(), east_counters_sc, east_shutoffs);
765  // hTopCountersShutOffs->SetPoint(hTopCountersShutOffs->GetN(), top_counters_sc, top_shutoffs);
766  // hBottomCountersShutOffs->SetPoint(hBottomCountersShutOffs->GetN(), bottom_counters_sc, bottom_shutoffs);
767 
768  return;
769 }
void analyze(const art::Event &evt)
enum BeamMode kRed
int status
Definition: fabricate.py:1613
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
void reconfigure(const fhicl::ParameterSet &pset)
TCanvas * canv
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
double mc6int
MC6 intensity [ppp].
Definition: MCenterData.h:25
const XML_Char * s
Definition: expat.h:262
unsigned long long fTriggerHeader_MasterTriggerNumber
Definition: RawTrigger.h:35
art::ServiceHandle< geo::Geometry > fGeo
std::string GetName(int i)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
T get(std::string const &key) const
Definition: ParameterSet.h:231
double mc6col
MC6 vertical collimator opening [mm].
Definition: MCenterData.h:26
int evt
art::ServiceHandle< art::TFileService > fFileService
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
T * makeAndRegister(char const *name, char const *title, ARGS...args) const
OStream cout
Definition: OStream.cxx:6
std::vector< std::vector< unsigned long long > > fShutOffTimes
int getDCM(daqchannelmap::dchan)
Definition: CMap.cxx:191
int getFEB(daqchannelmap::dchan)
Definition: CMap.cxx:206
std::map< unsigned int, TH1I * > hFEBSpillProfile
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
Data resulting from a Hough transform on the cell hit positions.
unsigned long long fTriggerHeader_TriggerNumber
Definition: RawTrigger.h:34
DetectorRateShutOff(const fhicl::ParameterSet &pset)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::vector< std::pair< double, double > > fHoughY
Definition: structs.h:12
std::vector< std::pair< double, double > > fHoughX
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
uint32_t fTriggerRange_TriggerLength
Definition: RawTrigger.h:40
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
uint32_t dchan
< DAQ Channel Map Package
std::vector< std::vector< unsigned long long > > fHitTimes
Definition: fwd.h:28
art::ServiceHandle< cmap::CMap > fChannelMap
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string
unsigned int uint