All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CrtNoiseMonTool_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CrtNoiseMonTool
3 // Module Type: analyzer
4 // File: CrtNoiseMonTool_module.cc
5 // Description: Makes a tree with waveform information.
6 ////////////////////////////////////////////////////////////////////////
7 
8 #include "art/Framework/Core/EDAnalyzer.h"
9 #include "art/Framework/Core/ModuleMacros.h"
10 #include "art/Framework/Principal/Event.h"
11 #include "art/Framework/Principal/Handle.h"
12 
13 #include "canvas/Utilities/Exception.h"
14 
15 #include "sbndaq-artdaq-core/Overlays/Common/BernCRTFragment.hh"
16 #include "artdaq-core/Data/Fragment.hh"
17 #include "artdaq-core/Data/ContainerFragment.hh"
18 #include "sbndaq-artdaq-core/Overlays/FragmentType.hh"
19 
20 #include "sbndaq-artdaq-core/Overlays/Common/BernCRTTranslator.hh"
21 
22 //#include "art/Framework/Services/Optional/TFileService.h"
23 #include "art_root_io/TFileService.h"
24 #include "TH1F.h"
25 #include "TH2F.h"
26 #include "TNtuple.h"
27 #include "TGraph.h"
28 
29 #include <algorithm>
30 #include <cassert>
31 #include <cmath>
32 #include <fstream>
33 #include <iomanip>
34 #include <vector>
35 #include <iostream>
36 
37 namespace icarus {
38  namespace crt {
39  class CrtNoiseMonTool;
40  }
41 }
42 
43 class icarus::crt::CrtNoiseMonTool : public art::EDAnalyzer {
44 
45 public:
46  explicit CrtNoiseMonTool(fhicl::ParameterSet const & pset); // explicit doesn't allow for copy initialization
47  virtual ~CrtNoiseMonTool();
48 
49  virtual void endJob();
50  virtual void analyze(art::Event const & evt);
51 
52 
53 private:
54  // std::map<uint8_t,std::vector<float>> macToRate;
55  //std::map<uint8_t,std::vector<uint64_t>> macToRateTime;
56  bool IsSideCRT(icarus::crt::BernCRTTranslator & hit);
57  std::map<uint8_t,std::vector<float>> macToRateSide;
58  std::map<uint8_t,std::vector<uint64_t>> macToRateSideTime;
59  std::map<uint8_t,std::vector<float>> macToRateTop;
60  std::map<uint8_t,std::vector<uint64_t>> macToRateTopTime;
61 
62  //AA: this list will need to be updated (and perhaps moved to a separate file)
63  // const uint8_t macsWestIn[6] = {21,19,17,22,23,24};
64  //const uint8_t macsWestOut[6] = {13,11,10,14,15,16};
65  //const uint8_t macsNorthIn[4] = {7,6,3,1};
66  //const uint8_t macsNorthOut[4] = {9,8,5,4};
67 
68  //TTree * tree;
69  //std::map<std::pair<short,short>,TH1F*>> macChanToHist;
70  //std::map<short,short> macToIndex;
71 
72 //data payload
73  /* uint8_t mac5; //last 8 bits of FEB mac5 address
74  uint16_t flags;
75  uint16_t lostcpu;
76  uint16_t lostfpga;*/
77  //uint32_t ts0;
78  //uint32_t ts1;
79  //uint16_t adc[32];
80  //uint32_t coinc;
81 
82 //metadata
83  /*uint64_t run_start_time;
84  uint64_t this_poll_start;
85  uint64_t this_poll_end;
86  uint64_t last_poll_start;
87  uint64_t last_poll_end;
88  int32_t system_clock_deviation;
89  uint32_t feb_tree_per_poll;
90  uint32_t feb_event_number;*/
91 
92  //information from fragment header
93  //uint32_t sequence_id;
94  //uint64_t fragment_timestamp;
95 
96 };
97 
98 //Define the constructor
99 icarus::crt::CrtNoiseMonTool::CrtNoiseMonTool(fhicl::ParameterSet const & pset)
100  : EDAnalyzer(pset)
101 {
102 
103  //this is how you setup/use the TFileService
104  //I do this here in the constructor to setup the histogram just once
105  //but you can call the TFileService and make new objects from anywhere
106  //art::ServiceHandle<art::TFileService> tfs; //pointer to a file named tfs
107 
108  //tree = tfs->make<TTree>("tree", "FEB tree");
109 
110 //event data
111  /*tree->Branch("mac5", &mac5, "mac5/b");
112  tree->Branch("flags", &flags, "flags/s");
113  tree->Branch("lostcpu", &lostcpu, "lostcpu/s");
114  tree->Branch("lostfpga", &lostcpu, "lostfpga/s");
115  tree->Branch("ts0", &ts0, "ts0/i");
116  tree->Branch("ts1", &ts1, "ts1/i");
117  tree->Branch("adc", &adc, "adc[32]/s");
118  tree->Branch("coinc", &coinc, "coinc/i");*/
119 
120 //metadata
121  /*tree->Branch("run_start_time", &run_start_time, "run_start_time/l");
122  tree->Branch("this_poll_start", &this_poll_start, "this_poll_start/l");
123  tree->Branch("this_poll_end", &this_poll_end, "this_poll_end/l");
124  tree->Branch("last_poll_start", &last_poll_start, "last_poll_start/l");
125  tree->Branch("last_poll_end", &last_poll_end, "last_poll_end/l");
126  tree->Branch("system_clock_deviation", &system_clock_deviation, "system_clock_deviation/I");
127  tree->Branch("feb_tree_per_poll", &feb_tree_per_poll, "feb_tree_per_poll/i");
128  tree->Branch("feb_event_number", &feb_event_number, "feb_event_number/i");
129 
130  tree->Branch("sequence_id", &sequence_id, "sequence_id/i");
131  tree->Branch("fragment_timestamp", &fragment_timestamp, "fragment_timestamp/l");*/
132 }
133 
135 {
136 }
137 
138 bool icarus::crt::CrtNoiseMonTool::IsSideCRT(icarus::crt::BernCRTTranslator & hit) {
139  /**
140  * Fragment ID described in SBN doc 16111
141  **/
142  return (hit.fragment_ID & 0x3100) == 0x3100;
143 }
144 
145 void icarus::crt::CrtNoiseMonTool::analyze(art::Event const & evt) {
146 
147  //WK 09/02/21. Update to BernCRTTranslator in sbndaq_artdaq_core
148  std::vector<icarus::crt::BernCRTTranslator> hit_vector;
149 
150  auto fragmentHandles = evt.getMany<artdaq::Fragments>();
151  for (auto handle : fragmentHandles) {
152  if (!handle.isValid() || handle->size() == 0)
153  continue;
154 
155  auto this_hit_vector = icarus::crt::BernCRTTranslator::getCRTData(*handle);
156 
157  hit_vector.insert(hit_vector.end(),this_hit_vector.begin(),this_hit_vector.end());
158 
159  }
160 
161  for(auto & hit : hit_vector) {
162 
163  uint8_t mac5 = hit.mac5;
164 
165  if(hit.flags==7) continue; //AA: 1. is it a correct check? 2. does it even work? :(
166 
167  const uint64_t runtime = (hit.this_poll_start-hit.run_start_time) *1.0e-9; //time since run start [s]
168  const double dt = (hit.this_poll_start-hit.last_poll_start)*1.0e-9; //interval between consecutive polls [s]
169  const float rate = 1.0e-3*hit.hits_in_poll/dt; //poll-averaged event rate [kHz]
170 
171  //avoid processing the same poll multiple times (they all have the same
172  // if(macToRateTime[mac5].size()!=0 && macToRateTime[mac5].back() == runtime) continue;
173  if(IsSideCRT(hit)){
174  //avoid processing the same poll multiple times (they all have the same
175  if(macToRateSideTime[mac5].size()!=0 && macToRateSideTime[mac5].back() == runtime) continue;
176  macToRateSideTime[mac5].push_back(runtime);
177  macToRateSide[mac5].push_back(rate);
178  }else{
179  //avoid processing the same poll multiple times (they all have the same
180  if(macToRateTopTime[mac5].size()!=0 && macToRateTopTime[mac5].back() == runtime) continue;
181  macToRateTopTime[mac5].push_back(runtime);
182  macToRateTop[mac5].push_back(rate);
183  }
184  }
185 }//end analyze
186 
187 
189 
190  art::ServiceHandle<art::TFileService> tfs; //pointer to a file named tfs
191  /*
192  std::string hviewName = "hview_";
193  std::string hviewTitle = "noise rates: ";
194  char reg ='e';//, lay='e';
195  if(macToRate.find(macsWestIn[0])!=macToRate.end()) {
196  hviewName+="w_i";
197  hviewTitle+="West Inner";
198  reg = 'w';
199  //lay = 'i';
200  }
201  else if(macToRate.find(macsWestOut[0])!=macToRate.end()) {
202  hviewName+="w_o";
203  hviewTitle+="West Outer";
204  reg = 'w';
205  //lay = 'o';
206  }
207  else if(macToRate.find(macsNorthIn[0])!=macToRate.end()) {
208  hviewName+="n_i";
209  hviewTitle+="North Inner";
210  reg = 'n';
211  //lay = 'i';
212  }
213  else if(macToRate.find(macsNorthOut[0])!=macToRate.end()) {
214  hviewName+="n_o";
215  hviewTitle+="North Outer";
216  reg = 'n';
217  //lay = 'o';
218  }
219  else{
220  std::cout << "no matching mac5s found. reg and lay not set!" << std::endl;
221  }
222 
223 
224  size_t yhigh=0;
225  if(reg=='n') yhigh = 2;
226  else yhigh = 3;
227  TH2F *hview = tfs->make<TH2F>(hviewName.c_str(),hviewTitle.c_str(),2,0,2,yhigh,0,yhigh);
228  if(reg=='w'){
229  hview->GetXaxis()->SetBinLabel(1,"north");
230  hview->GetXaxis()->SetBinLabel(2,"south");
231  hview->GetYaxis()->SetBinLabel(1,"pit");
232  hview->GetYaxis()->SetBinLabel(2,"mezzpit");
233  hview->GetYaxis()->SetBinLabel(3,"mezz");
234  }
235  if(reg=='n'){
236  hview->GetXaxis()->SetBinLabel(1,"east");
237  hview->GetXaxis()->SetBinLabel(2,"west");
238  hview->GetYaxis()->SetBinLabel(1,"pit");
239  hview->GetYaxis()->SetBinLabel(2,"mezz");
240  }
241  */
242  //size_t imac=0;
243  for(auto it=macToRateSide.begin(); it!=macToRateSide.end(); it++){
244 
245  uint8_t mac = (*it).first;
246 
247  std::string hname = "hside_"+std::to_string(mac);
248  std::string htitle = "noise rate: mac5 "+std::to_string(mac)+" [Side CRT]";
249  TH1F *h = tfs->make<TH1F>(hname.c_str(),htitle.c_str(),4000,0.0,40.0);
250  h->GetXaxis()->SetTitle("rate [kHz]");
251  h->GetXaxis()->SetTitleSize(0.04);
252  h->GetXaxis()->SetLabelSize(0.04);
253  h->GetYaxis()->SetLabelSize(0.04);
254  h->SetLineWidth(2);
255 
256  const size_t n = macToRateSide[mac].size();
257  float* rate = new float[n];
258  float* time = new float[n];
259  float meanRate =0.0;
260  for(size_t i=0; i<n; i++){
261  rate[i] = macToRateSide[mac][i];
262  time[i] = macToRateSideTime[mac][i];
263  meanRate+=rate[i];
264 
265  h->Fill(rate[i]);
266  }
267  std::cout << "mean rate before division: " << meanRate << std::endl;
268  meanRate*=1.0/n;
269  std::cout << "mean rate averaged over " << n << " points: " << meanRate << std::endl;
270  meanRate*=1.0/30;
271 
272  std::cout << "relative (to 30kHz) average "+htitle << " = " << meanRate << std::endl;
273 
274  std::string gname = "gside_"+std::to_string(mac);
275  std::string gtitle = "noise rate: mac5 "+std::to_string(mac)+" [Side CRT]";
276  TGraph *g = tfs->make<TGraph>(n,time,rate);
277  g->SetName(gname.c_str());
278  g->SetTitle(gtitle.c_str());
279  g->GetXaxis()->SetTitle("run time [s]");
280  g->GetYaxis()->SetTitle("rate [kHz]");
281  g->GetXaxis()->SetTitleSize(0.04);
282  g->GetYaxis()->SetTitleSize(0.04);
283  g->GetXaxis()->SetLabelSize(0.04);
284  g->GetYaxis()->SetLabelSize(0.04);
285 
286  g->SetMarkerStyle(8);
287 
288  h->Write();
289  g->Write();
290 
291  delete[] rate;
292  delete[] time;
293  }
294 
295  for(auto it=macToRateTop.begin(); it!=macToRateTop.end(); it++){
296 
297  uint8_t mac = (*it).first;
298 
299  std::string hname = "htop_"+std::to_string(mac);
300  std::string htitle = "noise rate: mac5 "+std::to_string(mac)+" [Top CRT]";
301  TH1F *h = tfs->make<TH1F>(hname.c_str(),htitle.c_str(),4000,0.0,40.0);
302  h->GetXaxis()->SetTitle("rate [kHz]");
303  h->GetXaxis()->SetTitleSize(0.04);
304  h->GetXaxis()->SetLabelSize(0.04);
305  h->GetYaxis()->SetLabelSize(0.04);
306  h->SetLineWidth(2);
307 
308  const size_t n = macToRateTop[mac].size();
309  float* rate = new float[n];
310  float* time = new float[n];
311  float meanRate =0.0;
312  for(size_t i=0; i<n; i++){
313  rate[i] = macToRateTop[mac][i];
314  time[i] = macToRateTopTime[mac][i];
315  meanRate+=rate[i];
316 
317  h->Fill(rate[i]);
318  }
319  std::cout << "mean rate before division: " << meanRate << std::endl;
320  meanRate*=1.0/n;
321  std::cout << "mean rate averaged over " << n << " points: " << meanRate << std::endl;
322  meanRate*=1.0/30;
323 
324  std::cout << "relative (to 30kHz) average "+htitle << " = " << meanRate << std::endl;
325 
326  std::string gname = "gtop_"+std::to_string(mac);
327  std::string gtitle = "noise rate: mac5 "+std::to_string(mac)+" [Top CRT]";
328  TGraph *g = tfs->make<TGraph>(n,time,rate);
329  g->SetName(gname.c_str());
330  g->SetTitle(gtitle.c_str());
331  g->GetXaxis()->SetTitle("run time [s]");
332  g->GetYaxis()->SetTitle("rate [kHz]");
333  g->GetXaxis()->SetTitleSize(0.04);
334  g->GetYaxis()->SetTitleSize(0.04);
335  g->GetXaxis()->SetLabelSize(0.04);
336  g->GetYaxis()->SetLabelSize(0.04);
337 
338  g->SetMarkerStyle(8);
339 
340  h->Write();
341  g->Write();
342 
343  delete[] rate;
344  delete[] time;
345  }
346 
347  // hview->Write();
348 
349 } //endJob
350 
351 DEFINE_ART_MODULE(icarus::crt::CrtNoiseMonTool)
352 //this is where the name is specified
CrtNoiseMonTool(fhicl::ParameterSet const &pset)
std::map< uint8_t, std::vector< uint64_t > > macToRateSideTime
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
process_name hit
Definition: cheaterreco.fcl:51
BEGIN_PROLOG g
std::map< uint8_t, std::vector< float > > macToRateSide
while getopts h
virtual void analyze(art::Event const &evt)
std::map< uint8_t, std::vector< uint64_t > > macToRateTopTime
bool IsSideCRT(icarus::crt::BernCRTTranslator &hit)
std::string to_string(WindowPattern const &pattern)
do i e
art::ServiceHandle< art::TFileService > tfs
TCEvent evt
Definition: DataStructs.cxx:8
std::map< uint8_t, std::vector< float > > macToRateTop
process_name crt
BEGIN_PROLOG could also be cout