All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CrtCalAnalysis_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CrtCalAnalyzer
3 // Module Type: analyzer
4 // File: CrtCalAnalyzer_module.cc
5 // Description: Makes a tree with waveform information.
6 ////////////////////////////////////////////////////////////////////////
7 
8 //framwork includes
9 #include "art/Framework/Core/EDAnalyzer.h"
10 #include "art/Framework/Core/ModuleMacros.h"
11 #include "art/Framework/Principal/Event.h"
12 #include "art/Framework/Principal/Handle.h"
13 
14 #include "fhiclcpp/ParameterSet.h"
15 #include "fhiclcpp/types/Table.h"
16 #include "fhiclcpp/types/Atom.h"
17 
18 #include "canvas/Utilities/Exception.h"
19 
20 //local includes
21 #include "sbndaq-artdaq-core/Overlays/Common/BernCRTFragment.hh"
22 #include "artdaq-core/Data/Fragment.hh"
23 #include "artdaq-core/Data/ContainerFragment.hh"
24 #include "sbndaq-artdaq-core/Overlays/FragmentType.hh"
25 #include "icaruscode/CRT/CRTDecoder/CrtCal.h"
26 
27 #include "sbndaq-artdaq-core/Overlays/Common/BernCRTTranslator.hh"
28 
29 //ROOT includes
30 #include "art_root_io/TFileService.h"
31 #include "TH1F.h"
32 #include "TNtuple.h"
33 
34 //c++ includes
35 #include <algorithm>
36 #include <cassert>
37 #include <cmath>
38 #include <fstream>
39 #include <iomanip>
40 #include <vector>
41 #include <iostream>
42 #include <map>
43 #include <string>
44 
45 namespace icarus {
46  namespace crt {
47  class CrtCalAnalyzer;
48  }
49 }
50 
51 using std::map;
52 using std::vector;
53 using std::string;
54 using std::to_string;
55 using icarus::crt::CrtCal;
56 
57 class icarus::crt::CrtCalAnalyzer : public art::EDAnalyzer {
58 
59 public:
60  struct Config {
61 
62  // Save some typing:
63  using Name = fhicl::Name;
64  using Comment = fhicl::Comment;
65 
66  // One Atom for each parameter
67 
68  fhicl::Sequence<uint8_t> Mac {
69  Name("mac5"),
70  Comment("mac addresses for each FEB in the diasychain")
71  };
72  };
73 
74  using Parameters = art::EDAnalyzer::Table<Config>;
75 
76  explicit CrtCalAnalyzer(Parameters const& config); // explicit doesn't allow for copy initialization
77  virtual ~CrtCalAnalyzer();
78 
79  virtual void beginJob();
80  virtual void analyze(art::Event const & evt);
81  virtual void endJob();
82 
83 
84 private:
85  art::ServiceHandle<art::TFileService> tfs;
86 
87  vector<uint8_t> macs;
88  map<uint8_t,vector<TH1F*>*> macToHistos;
89  TTree* calTree;
90 
91 //vars for calTree
92 
93  uint8_t fMac5;
94  bool fActive[32];
95  float fGain[32];
96  float fGainErr[32];
97  float fGainXsqr[32];
98  short fGainNdf[32];
99  float fGainPed[32];
100  float fGainPedErr[32];
101  short fNpeak[32];
102  float fPed[32];
103  float fPedErr[32];
104  float fPedXsqr[32];
105  short fPedNdf[32];
106  float fPedSigma[32];
107  float fPedSigmaErr[32];
108  float fPedNorm[32];
109  float fPedNormErr[32];
110  int fThreshADC[32];
111  float fThreshPE[32];
112  int fNabove[32];
113  float fPeakNorm[32][5];
114  float fPeakNormErr[32][5];
115  float fPeakSigma[32][5];
116  float fPeakSigmaErr[32][5];
117  float fPeakMean[32][5];
118  float fPeakMeanErr[32][5];
119  float fPeakXsqr[32][5];
120  short fPeakNdf[32][5];
121 };
122 
123 //Define the constructor
124 icarus::crt::CrtCalAnalyzer::CrtCalAnalyzer(Parameters const& config)
125  : EDAnalyzer(config) , macs(config().Mac())
126 {
127  fMac5 = 0;
128  for(int i=0; i<(int)macs.size(); i++){
129  macToHistos[macs[i]] = new vector<TH1F*>();
130  for(int ch=0; ch<32; ch++){
131  string hname = "hadc_"+to_string(macs[i])+"_"+to_string(ch);
132  string htitle = "raw charge: mac5 "+to_string(macs[i])+", ch. "+to_string(ch);
133  macToHistos[macs[i]]->push_back(tfs->make<TH1F>(hname.c_str(),htitle.c_str(),4100,0,4100));
134  }
135  }
136 
137  for(int i=0; i<32; i++) {
138  fActive[i] = false;
139  fGain[i] = FLT_MAX;
140  fGainErr[i] = FLT_MAX;
141  fGainXsqr[i] = FLT_MAX;
142  fGainNdf[i] = SHRT_MAX;
143  fGainPed[i] = FLT_MAX;
144  fGainPedErr[i] = FLT_MAX;
145  fNpeak[i] = SHRT_MAX;
146  fPed[i] = FLT_MAX;
147  fPedErr[i] = FLT_MAX;
148  fPedXsqr[i] = FLT_MAX;
149  fPedNdf[i] = SHRT_MAX;
150  fPedSigma[i] = FLT_MAX;
151  fPedSigmaErr[i] = FLT_MAX;
152  fPedNorm[i] = FLT_MAX;
153  fPedNormErr[i] = FLT_MAX;
154  fThreshADC[i] = INT_MAX;
155  fThreshPE[i] = FLT_MAX;
156  fNabove[i] = INT_MAX;
157 
158  for(size_t j=0; j<5; j++){
159  fPeakNorm[i][j] = FLT_MAX;
160  fPeakNormErr[i][j] = FLT_MAX;
161  fPeakSigma[i][j] = FLT_MAX;
162  fPeakSigmaErr[i][j] = FLT_MAX;
163  fPeakMean[i][j] = FLT_MAX;
164  fPeakMeanErr[i][j] = FLT_MAX;
165  fPeakXsqr[i][j] = FLT_MAX;
166  fPeakNdf[i][j] = SHRT_MAX;
167  }
168  }
169 
170 }
171 
172 icarus::crt::CrtCalAnalyzer::~CrtCalAnalyzer()
173 {
174 }
175 
176 void icarus::crt::CrtCalAnalyzer::beginJob(){
177 
178  calTree = tfs->make<TTree>("calAnalyzerTree", "SiPM/FEB channel calibration data");
179 
180  calTree->Branch("mac5", &fMac5, "mac5/b");
181  calTree->Branch("active", fActive, "active[32]/O");
182  calTree->Branch("gain", fGain, "gain[32]/F");
183  calTree->Branch("gainErr", fGainErr, "gainErr[32]/F");
184  calTree->Branch("gainXsqr", fGainXsqr, "gainXsqr[32]/F");
185  calTree->Branch("gainNdf", fGainNdf, "gainNdf[32]/s");
186  calTree->Branch("gainPed", fGainPed, "gainPed[32]/F");
187  calTree->Branch("gainPedErr", fGainPedErr, "gainPedErr[32]/F");
188  calTree->Branch("nPeak", fNpeak, "nPeak[32]/s");
189  calTree->Branch("peakXsqr", fPeakXsqr, "peakXsqr[32][5]/F");
190  calTree->Branch("peakNdf", fPeakNdf, "peakNdf[32][5]/s");
191  calTree->Branch("peakMean", fPeakMean, "peakMean[32][5]/F");
192  calTree->Branch("peakMeanErr", fPeakMeanErr, "peakMeanErr[32][5]/F");
193  calTree->Branch("peakNorm", fPeakNorm, "peakNorm[32][5]/F");
194  calTree->Branch("peakNormErr", fPeakNormErr, "peakNormErr[32][5]/F");
195  calTree->Branch("peakSigma", fPeakSigma, "peakSigma[32][5]/F");
196  calTree->Branch("peakSigmaErr", fPeakSigmaErr, "peakSigmaErr[32][5]/F");
197  calTree->Branch("ped", fPed, "ped[32]/F");
198  calTree->Branch("pedErr", fPedErr, "pedErr[32]/F");
199  calTree->Branch("pedXsqr", fPedXsqr, "pedXsqr[32]/F");
200  calTree->Branch("pedNdf", fPedNdf, "pedNdf[32]/s");
201  calTree->Branch("pedSigma", fPedSigma, "pedSigma[32]/F");
202  calTree->Branch("pedSigmaErr", fPedSigmaErr, "pedSigmaErr[32]/F");
203  calTree->Branch("pedNorm", fPedNorm, "pedNorm[32]/F");
204  calTree->Branch("pedNormErr", fPedNormErr, "pedNormErr[32]/F");
205  calTree->Branch("threshAdc", fThreshADC, "threshAdc[32]/I");
206  calTree->Branch("threshPe", fThreshPE, "threshPe[32]/F");
207  calTree->Branch("nAbove", fNabove, "nAbove[32]/I");
208 
209 }
210 
211 void icarus::crt::CrtCalAnalyzer::analyze(art::Event const & evt)
212 {
213  //WK 09/02/21. Update to BernCRTTranslator in sbndaq_artdaq_core
214  std::vector<icarus::crt::BernCRTTranslator> hit_vector;
215 
216  auto fragmentHandles = evt.getMany<artdaq::Fragments>();
217  for (auto handle : fragmentHandles) {
218  if (!handle.isValid() || handle->size() == 0)
219  continue;
220 
221  auto this_hit_vector = icarus::crt::BernCRTTranslator::getCRTData(*handle);
222 
223  hit_vector.insert(hit_vector.end(),this_hit_vector.begin(),this_hit_vector.end());
224 
225  }
226 
227  for(auto & hit : hit_vector) {
228 
229  fMac5 = hit.mac5;
230 
231  if(macToHistos.find(fMac5) == macToHistos.end()) return;
232 
233  for(int ch=0; ch<32; ch++) {
234  macToHistos[fMac5]->at(ch)->Fill( hit.adc[ch] );
235  }
236  }
237 }//end analyze
238 
239 
240 void icarus::crt::CrtCalAnalyzer::endJob(){
241 
242  std::cout << "done filling histograms..." << std::endl;
243  std::cout << "found " << macToHistos.size() << " FEBs" << std::endl;
244  if(!macToHistos.begin()->second->empty()){
245  std::cout << "first histo size: " << macToHistos.begin()->second->at(0)->Integral()
246  << " entries" << std::endl;
247  }
248  else {
249  std::cout << "hist vect is empty!" << std::endl;
250  }
251 
252  for(auto const& macHist : macToHistos){
253 
254  uint8_t mac = macHist.first;
255  vector<TH1F*>* hvec = macHist.second;
256 
257  std::cout << "construct instance of CrtCal for mac5 " << (short)mac << std::endl;
258  CrtCal* cal = new CrtCal(hvec);
259  cal->Cal();
260 
261  std::cout << "retreiving cal data..." << std::endl;
262  fMac5 = macHist.first;
263  bool* ptrActive = cal->GetActive();
264  float* ptrGain = cal->GetGain();
265  float* ptrGainErr = cal->GetGainErr();
266  float* ptrGainXsqr = cal->GetGainXsqr();
267  short* ptrGainNdf = cal->GetGainNdf();
268  float* ptrGainPed = cal->GetGainPed();
269  float* ptrGainPedErr = cal->GetGainPedErr();
270  short* ptrNpeak = cal->GetNpeak();
271  float** ptrPeakNorm = cal->GetPeakNorm();
272  float** ptrPeakNormErr = cal->GetPeakNormErr();
273  float** ptrPeakSigma = cal->GetPeakSigma();
274  float** ptrPeakSigmaErr = cal->GetPeakSigmaErr();
275  float** ptrPeakMean = cal->GetPeakMean();
276  float** ptrPeakMeanErr = cal->GetPeakMeanErr();
277  float** ptrPeakXsqr = cal->GetPeakXsqr();
278  short** ptrPeakNdf = cal->GetPeakNdf();
279  float* ptrPed = cal->GetPed();
280  float* ptrPedErr = cal->GetPedErr();
281  float* ptrPedXsqr = cal->GetPedXsqr();
282  short* ptrPedNdf = cal->GetPedNdf();
283  float* ptrPedSigma = cal->GetPedSigma();
284  float* ptrPedSigmaErr = cal->GetPedSigmaErr();
285  float* ptrPedNorm = cal->GetPedNorm();
286  float* ptrPedNormErr = cal->GetPedNormErr();
287  int* ptrThreshADC = cal->GetThreshADC();
288  float* ptrThreshPE = cal->GetThreshPE();
289  int* ptrNabove = cal->GetNabove();
290 
291  //copy values to local arrays
292  for(size_t i=0; i<32; i++){
293  if(ptrActive!=nullptr) fActive[i] = ptrActive[i];
294  if(ptrGain!=nullptr) fGain[i] = ptrGain[i];
295  if(ptrGainErr!=nullptr) fGainErr[i] = ptrGainErr[i];
296  if(ptrGainXsqr!=nullptr) fGainXsqr[i] = ptrGainXsqr[i];
297  if(ptrGainNdf!=nullptr) fGainNdf[i] = ptrGainNdf[i];
298  if(ptrGainPed!=nullptr) fGainPed[i] = ptrGainPed[i];
299  if(ptrGainPedErr!=nullptr) fGainPedErr[i] = ptrGainPedErr[i];
300  if(ptrNpeak!=nullptr) fNpeak[i] = ptrNpeak[i];
301  if(ptrPed!=nullptr) fPed[i] = ptrPed[i];
302  if(ptrPedErr!=nullptr) fPedErr[i] = ptrPedErr[i];
303  if(ptrPedXsqr!=nullptr) fPedXsqr[i] = ptrPedXsqr[i];
304  if(ptrPedNdf!=nullptr) fPedNdf[i] = ptrPedNdf[i];
305  if(ptrPedSigma!=nullptr) fPedSigma[i] = ptrPedSigma[i];
306  if(ptrPedSigmaErr!=nullptr) fPedSigmaErr[i] = ptrPedSigmaErr[i];
307  if(ptrPedNorm!=nullptr) fPedNorm[i] = ptrPedNorm[i];
308  if(ptrPedNormErr!=nullptr) fPedNormErr[i] = ptrPedNormErr[i];
309  if(ptrThreshADC!=nullptr) fThreshADC[i] = ptrThreshADC[i];
310  if(ptrThreshPE!=nullptr) fThreshPE[i] = ptrThreshPE[i];
311  if(ptrNabove!=nullptr) fNabove[i] = ptrNabove[i];
312 
313  for(size_t j=0; j<5; j++){
314  if(ptrPeakNorm!=nullptr) fPeakNorm[i][j] = ptrPeakNorm[i][j];
315  if(ptrPeakNormErr!=nullptr) fPeakNormErr[i][j] = ptrPeakNormErr[i][j];
316  if(ptrPeakSigma!=nullptr) fPeakSigma[i][j] = ptrPeakSigma[i][j];
317  if(ptrPeakSigmaErr!=nullptr) fPeakSigmaErr[i][j] = ptrPeakSigmaErr[i][j];
318  if(ptrPeakMean!=nullptr) fPeakMean[i][j] = ptrPeakMean[i][j];
319  if(ptrPeakMeanErr!=nullptr) fPeakMeanErr[i][j] = ptrPeakMeanErr[i][j];
320  if(ptrPeakXsqr!=nullptr) fPeakXsqr[i][j] = ptrPeakXsqr[i][j];
321  if(ptrPeakNdf!=nullptr) fPeakNdf[i][j] = ptrPeakNdf[i][j];
322  }
323  }
324 
325  std::cout << "fill tree event" << std::endl;
326  calTree->Fill();
327 
328  delete cal;
329  }
330 
331  //calTree->Write();
332 
333 }
334 
335 DEFINE_ART_MODULE(icarus::crt::CrtCalAnalyzer)
336 //this is where the name is specified
hit
process_name hit
Definition: cheaterreco.fcl:51
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
Name
BEGIN_PROLOG vertical distance to the surface Name
Definition: geometry_icarus.fcl:82
icarus::crt::CrtCalAnalyzer::macToHistos
map< uint8_t, vector< TH1F * > * > macToHistos
Definition: CrtCalAnalysis_module.cc:88
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
icarus::crt::CrtCalAnalyzer::analyze
virtual void analyze(art::Event const &evt)
Definition: CrtCalAnalysis_module.cc:211
icarus::crt::CrtCalAnalyzer::Parameters
art::EDAnalyzer::Table< Config > Parameters
Definition: CrtCalAnalysis_module.cc:74
icarus::crt::CrtCalAnalyzer::CrtCalAnalyzer
CrtCalAnalyzer(Parameters const &config)
Definition: CrtCalAnalysis_module.cc:124
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
icarus::crt::CrtCalAnalyzer::tfs
art::ServiceHandle< art::TFileService > tfs
Definition: CrtCalAnalysis_module.cc:85
crt
process_name crt
Definition: run_genericCRT.fcl:82
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
Generated by   doxygen 1.8.5