BeamStructureAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file BeamStructureAna_module.cc
3 /// \module analyzer
4 /// \brief Basic analyzer module to look at the raw data within a file.
5 /// \author Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
6 /// \date February 2020
7 ////////////////////////////////////////////////////////////////////////////
8 
9 // framework
13 #include "fhiclcpp/ParameterSet.h"
21 
22 // nova
23 #include "RawData/RawBeamline.h"
25 #include "BeamlineRecoBase/ToF.h"
27 
28 // stl
29 #include <iostream>
30 #include <bitset>
31 
32 // root
33 #include "TStyle.h"
34 #include "TFile.h"
35 #include "TH2I.h"
36 #include "TProfile.h"
37 #include "TLine.h"
38 #include "TText.h"
39 #include "TCanvas.h"
40 
41 // -----------------------------------------------------------------------
42 namespace tbana {
43 
45 
46  public:
47 
49 
50  void reconfigure(const fhicl::ParameterSet& pset);
51  void analyze(const art::Event& evt);
52  void beginJob();
53  void endJob();
54  void beginSubRun(const art::SubRun& sr);
55 
56  private:
57 
58  // config
65 
66  // analysis
67  TH1I* hNumTCHits;
70  TH1I* hBLReco;
71  std::vector<std::vector<TH1I*> > hWCProfiles;
72  std::vector<std::vector<TH1I*> > hCumulativeWCProfiles;
73  unsigned int fTotSpills;
74  unsigned int fTotTriggers;
75 
76  // output
77  TFile* fOutFile;
78 
79  };
80 
81 }
82 
83 // -----------------------------------------------------------------------
85  this->reconfigure(pset);
86 }
87 
88 // -----------------------------------------------------------------------
90  fRawWCDataLabel = pset.get<std::string>("RawWCDataLabel");
91  fRecoTCDigitLabel = pset.get<std::string>("RecoTCDigitLabel");
92  fRecoWCDigitLabel = pset.get<std::string>("RecoWCDigitLabel");
93  fRecoToFDigitLabel = pset.get<std::string>("RecoToFDigitLabel");
94  fRecoToFLabel = pset.get<std::string>("RecoToFLabel");
95  fRecoWCTrackLabel = pset.get<std::string>("RecoWCTrackLabel");
96 }
97 
98 // -----------------------------------------------------------------------
100  gStyle->SetOptStat(0);
101  fTotSpills = 0;
102  fTotTriggers = 0;
103  fOutFile = new TFile("BeamStructureAna.root", "RECREATE");
104  hNumTCHits = new TH1I("NumTCHits", ";Number of Target Counter Hits;", 10, 0, 10);
105  hNumTCHits->GetXaxis()->SetNdivisions(10);
106  hNumTCHits->GetXaxis()->CenterLabels();
107  hNumTCHits->SetMinimum(0);
108  hNumBLDetHits = new TH1I("NumBLDetHits", ";Beamline Detector;Number of Hits;", 14, 0, 14);
109  hNumBLDetHits->GetXaxis()->SetBinLabel(1, "TC");
110  hNumBLDetHits->GetXaxis()->SetBinLabel(2, "ToF-US-1");
111  hNumBLDetHits->GetXaxis()->SetBinLabel(3, "ToF-US-2");
112  hNumBLDetHits->GetXaxis()->SetBinLabel(4, "ToF-US-3");
113  hNumBLDetHits->GetXaxis()->SetBinLabel(5, "ToF-US-4");
114  hNumBLDetHits->GetXaxis()->SetBinLabel(6, "PAD1");
115  hNumBLDetHits->GetXaxis()->SetBinLabel(7, "WC1");
116  hNumBLDetHits->GetXaxis()->SetBinLabel(8, "PAD2");
117  hNumBLDetHits->GetXaxis()->SetBinLabel(9, "PAD3");
118  hNumBLDetHits->GetXaxis()->SetBinLabel(10, "PAD4");
119  hNumBLDetHits->GetXaxis()->SetBinLabel(11, "ToF-DS-1");
120  hNumBLDetHits->GetXaxis()->SetBinLabel(12, "ToF-DS-2");
121  hNumBLDetHits->GetXaxis()->SetBinLabel(13, "ToF-DS-3");
122  hNumBLDetHits->GetXaxis()->SetBinLabel(14, "ToF-DS-4");
123  //hNumBLDetHits->SetMinimum(0);
124  hCumulativeBLDetHits = new TH1I("CumulativeBLDetHits", ";Cumulative Detectors;Number of Triggers;", 8, 0, 8);
125  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(1, "TC");
126  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(2, "+ToF-US");
127  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(3, "+PAD1");
128  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(4, "+WC1");
129  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(5, "+PAD2");
130  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(6, "+PAD3");
131  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(7, "+PAD4");
132  hCumulativeBLDetHits->GetXaxis()->SetBinLabel(8, "+ToF-DS");
133  //hCumulativeBLDetHits->SetMinimum(0);
134  hBLReco = new TH1I("BLReco", ";Beamline Reco;", 4, 0, 4);
135  hBLReco->GetXaxis()->SetBinLabel(1, "None");
136  hBLReco->GetXaxis()->SetBinLabel(2, "ToF");
137  hBLReco->GetXaxis()->SetBinLabel(3, "WCTrack");
138  hBLReco->GetXaxis()->SetBinLabel(4, "Both");
139  for (unsigned int wc = 0; wc < 4; ++wc) {
140  std::vector<TH1I*> wcProfiles, cumulativeWCProfiles;
141  for (unsigned int plane = 0; plane < 2; ++plane) {
142  wcProfiles.push_back(new TH1I(Form("WC%dPlane%dProfile", wc, plane), Form(";WC %d, Plane %d", wc, plane), 128, 0, 128));
143  cumulativeWCProfiles.push_back(new TH1I(Form("WC%dPlane%dCumulativeProfile", wc, plane), Form(";WC %d, Plane %d", wc, plane), 128, 0, 128));
144  }
145  hWCProfiles.push_back(wcProfiles);
146  hCumulativeWCProfiles.push_back(cumulativeWCProfiles);
147  }
148 }
149 
150 // -----------------------------------------------------------------------
152  ++fTotSpills;
153  return;
154 }
155 
156 // -----------------------------------------------------------------------
158  fOutFile->cd();
159  hNumTCHits->Write();
160  hNumBLDetHits->Write();
161  hCumulativeBLDetHits->Write();
162  hBLReco->Write();
163  for (unsigned int wc = 0; wc < 4; ++wc) {
164  for (unsigned int plane = 0; plane < 2; ++plane) {
165  hWCProfiles[wc][plane]->Write();
166  hCumulativeWCProfiles[wc][plane]->Write();
167  TCanvas* canv = new TCanvas(Form("WC%dPlane%d", wc, plane), Form("WC %d, Plane %d", wc, plane), 1600, 600);
168  canv->Divide(2,1);
169  canv->cd(1);
170  hWCProfiles[wc][plane]->Draw("hist");
171  canv->cd(2);
172  hCumulativeWCProfiles[wc][plane]->Draw("hist");
173  canv->Write();
174  delete canv;
175  }
176  }
177  fOutFile->Close();
178  delete fOutFile;
179 
180  std::cout << "Total spills: " << fTotSpills << "; total triggers: " << fTotTriggers << std::endl;
181  for (int bin = 1; bin <= hNumTCHits->GetNbinsX(); ++bin)
182  std::cout << " Triggers with " << hNumTCHits->GetBinLowEdge(bin) << " TC hits: " << hNumTCHits->GetBinContent(bin) << std::endl;
183 }
184 
185 // -----------------------------------------------------------------------
187 
188  ++fTotTriggers;
189 
190  // Get the raw WC data
192  std::vector<art::Ptr<rawdata::RawBeamlineWC> > wcs;
193  if (evt.getByLabel(fRawWCDataLabel, wcHandle))
194  art::fill_ptr_vector(wcs, wcHandle);
195 
196  // Get the reconstructed digits for the target counter
198  std::vector<art::Ptr<brb::BeamlineDigit> > scDigits;
199  if (evt.getByLabel(fRecoTCDigitLabel, scDigitHandle))
200  art::fill_ptr_vector(scDigits, scDigitHandle);
201 
202  // Get the reconstructed wc digits
204  std::vector<art::Ptr<brb::BeamlineDigit> > wcDigits;
205  if (evt.getByLabel(fRecoWCDigitLabel, wcDigitHandle))
206  art::fill_ptr_vector(wcDigits, wcDigitHandle);
207 
208  // Get reco ToF
210  std::vector<art::Ptr<brb::ToF> > tofs;
211  if (evt.getByLabel(fRecoToFLabel, tofHandle))
212  art::fill_ptr_vector(tofs, tofHandle);
213 
214  // Get reco WCTrack
216  std::vector<art::Ptr<brb::WCTrack> > wcTracks;
217  if (evt.getByLabel(fRecoWCTrackLabel, wcTrackHandle))
218  art::fill_ptr_vector(wcTracks, wcTrackHandle);
219 
220  // Get the reconstructed ToF digits
222  std::vector<art::Ptr<brb::BeamlineDigit> > tofDigits;
223  if (evt.getByLabel(fRecoToFDigitLabel, tofDigitHandle))
224  art::fill_ptr_vector(tofDigits, tofDigitHandle);
225 
226  // Fill wire chamber profiles
227  std::vector<std::vector<rawdata::RawBeamlineWC::WCPulse> > x_pulses;
228  std::vector<std::vector<rawdata::RawBeamlineWC::WCPulse> > y_pulses;
229  for (unsigned int wc = 0; wc < wcs.size(); ++wc) {
230  x_pulses.push_back(wcs[wc]->XPulses());
231  y_pulses.push_back(wcs[wc]->YPulses());
232  }
233  for (unsigned int wc = 0; wc < wcs.size(); ++wc) {
234  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator xPulseIt = x_pulses[wc].begin();
235  xPulseIt != x_pulses[wc].end(); ++xPulseIt)
236  hWCProfiles[wc][0]->Fill(xPulseIt->Channel);
237  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator yPulseIt = y_pulses[wc].begin();
238  yPulseIt != y_pulses[wc].end(); ++yPulseIt)
239  hWCProfiles[wc][1]->Fill(yPulseIt->Channel);
240  }
241 
242  // Get the number of hits on each SC
243  // 0 is the TC, 1-4 are the paddles
244  std::map<unsigned int, unsigned int> scDigitHits;
245  for (std::vector<art::Ptr<brb::BeamlineDigit> >::const_iterator scDigitIt = scDigits.begin();
246  scDigitIt != scDigits.end(); ++scDigitIt)
247  ++scDigitHits[(*scDigitIt)->ChannelID().Channel];
248 
249  // Get the number of hits on each WC digit
250  std::map<unsigned int, unsigned int> wcDigitHits;
251  for (std::vector<art::Ptr<brb::BeamlineDigit> >::const_iterator wcDigitIt = wcDigits.begin();
252  wcDigitIt != wcDigits.end(); ++wcDigitIt)
253  ++wcDigitHits[(*wcDigitIt)->ChannelID().Channel];
254 
255  // Get the number of hits on each ToF channel
256  std::map<unsigned int, unsigned int> usToFHits, dsToFHits;
257  for (std::vector<art::Ptr<brb::BeamlineDigit> >::const_iterator tofDigitIt = tofDigits.begin();
258  tofDigitIt != tofDigits.end(); ++tofDigitIt) {
259  if ((*tofDigitIt)->ChannelID().Detector == beamlinegeo::ToFCounter::US)
260  ++usToFHits[(*tofDigitIt)->ChannelID().Channel];
261  if ((*tofDigitIt)->ChannelID().Detector == beamlinegeo::ToFCounter::DS)
262  ++dsToFHits[(*tofDigitIt)->ChannelID().Channel];
263  }
264 
265  hNumTCHits->Fill(scDigitHits[0]);
266 
267  hNumBLDetHits->Fill("TC", scDigitHits[0]);
268  hNumBLDetHits->Fill("ToF-US-1", usToFHits[0]);
269  hNumBLDetHits->Fill("ToF-US-2", usToFHits[1]);
270  hNumBLDetHits->Fill("ToF-US-3", usToFHits[2]);
271  hNumBLDetHits->Fill("ToF-US-4", usToFHits[3]);
272  hNumBLDetHits->Fill("PAD1", scDigitHits[1]);
273  hNumBLDetHits->Fill("WC1", wcDigitHits[0]);
274  hNumBLDetHits->Fill("PAD2", scDigitHits[2]);
275  hNumBLDetHits->Fill("PAD3", scDigitHits[3]);
276  hNumBLDetHits->Fill("PAD4", scDigitHits[4]);
277  hNumBLDetHits->Fill("ToF-DS-1", dsToFHits[0]);
278  hNumBLDetHits->Fill("ToF-DS-2", dsToFHits[1]);
279  hNumBLDetHits->Fill("ToF-DS-3", dsToFHits[2]);
280  hNumBLDetHits->Fill("ToF-DS-4", dsToFHits[3]);
281 
282  if (scDigitHits[0]) {
283  hCumulativeBLDetHits->Fill("TC", 1);
284  if (usToFHits[0] and usToFHits[1] and usToFHits[2] and usToFHits[3]) {
285  hCumulativeBLDetHits->Fill("+ToF-US", 1);
286  if (scDigitHits[1]) {
287  hCumulativeBLDetHits->Fill("+PAD1", 1);
288  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator xPulseIt = x_pulses[0].begin();
289  xPulseIt != x_pulses[0].end(); ++xPulseIt)
290  hCumulativeWCProfiles[0][0]->Fill(xPulseIt->Channel);
291  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator yPulseIt = y_pulses[0].begin();
292  yPulseIt != y_pulses[0].end(); ++yPulseIt)
293  hCumulativeWCProfiles[0][1]->Fill(yPulseIt->Channel);
294  if (wcDigitHits[0]) {
295  hCumulativeBLDetHits->Fill("+WC1", 1);
296  if (scDigitHits[2]) {
297  hCumulativeBLDetHits->Fill("+PAD2", 1);
298  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator xPulseIt = x_pulses[1].begin();
299  xPulseIt != x_pulses[1].end(); ++xPulseIt)
300  hCumulativeWCProfiles[1][0]->Fill(xPulseIt->Channel);
301  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator yPulseIt = y_pulses[1].begin();
302  yPulseIt != y_pulses[1].end(); ++yPulseIt)
303  hCumulativeWCProfiles[1][1]->Fill(yPulseIt->Channel);
304  if (scDigitHits[3]) {
305  hCumulativeBLDetHits->Fill("+PAD3", 1);
306  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator xPulseIt = x_pulses[2].begin();
307  xPulseIt != x_pulses[2].end(); ++xPulseIt)
308  hCumulativeWCProfiles[2][1]->Fill(xPulseIt->Channel);
309  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator yPulseIt = y_pulses[2].begin();
310  yPulseIt != y_pulses[2].end(); ++yPulseIt)
311  hCumulativeWCProfiles[2][1]->Fill(yPulseIt->Channel);
312  if (scDigitHits[4]) {
313  hCumulativeBLDetHits->Fill("+PAD4", 1);
314  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator xPulseIt = x_pulses[3].begin();
315  xPulseIt != x_pulses[3].end(); ++xPulseIt)
316  hCumulativeWCProfiles[3][1]->Fill(xPulseIt->Channel);
317  for (std::vector<rawdata::RawBeamlineWC::WCPulse>::const_iterator yPulseIt = y_pulses[3].begin();
318  yPulseIt != y_pulses[3].end(); ++yPulseIt)
319  hCumulativeWCProfiles[3][1]->Fill(yPulseIt->Channel);
320  if (dsToFHits[0] and dsToFHits[1] and dsToFHits[2] and dsToFHits[3])
321  hCumulativeBLDetHits->Fill("+ToF-DS", 1);
322  }
323  }
324  }
325  }
326  }
327  }
328  }
329 
330  if (!tofs.size() and !wcTracks.size())
331  hBLReco->Fill("None", 1);
332  else if (tofs.size() and !wcTracks.size())
333  hBLReco->Fill("ToF", 1);
334  else if (!tofs.size() and wcTracks.size())
335  hBLReco->Fill("WCTrack", 1);
336  else if (tofs.size() and wcTracks.size())
337  hBLReco->Fill("Both", 1);
338 
339  return;
340 
341 }
342 
void beginSubRun(const art::SubRun &sr)
std::vector< std::vector< TH1I * > > hCumulativeWCProfiles
TCanvas * canv
DEFINE_ART_MODULE(TestTMapFile)
void reconfigure(const fhicl::ParameterSet &pset)
Encapsulation of reconstructed digitizer &#39;hits&#39;. Used for ToF PMTs and SiPMs, and Cherenkov and Muon ...
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
T get(std::string const &key) const
Definition: ParameterSet.h:231
Encapsulation of reconstructed Time-of-Flight (ToF) information. Part of beamline reconstruction for ...
int evt
BeamStructureAna(const fhicl::ParameterSet &pset)
std::vector< std::vector< TH1I * > > hWCProfiles
Encapsulation of reconstructed Wire Chamber track. Part of beamline reconstruction for NOvA test beam...
caf::StandardRecord * sr
void analyze(const art::Event &evt)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Raw data definitions for beamline data used in NOvA test beam experiment.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Definition: fwd.h:28
enum BeamMode string