ChannelPlots_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief TODO
3 /// \author shanahan@fnal.gov
4 ////////////////////////////////////////////////////////////////////////
5 #include <iostream>
6 #include <map>
7 #include <vector>
8 
9 // Root
10 #include "TH1F.h"
11 #include "TH2F.h"
12 
13 // Framework includes
22 
23 // NOvA includes
24 #include "Geometry/Geometry.h"
27 #include "RecoBase/Track.h"
28 #include "RecoBase/CellHit.h"
31 
32 
33 namespace dprf {
34 
35  class RawFebPlots;
36 
38  {
39  public:
40  explicit ChannelPlots(fhicl::ParameterSet const &p);
41  ~ChannelPlots();
42 
43  void beginRun(const art::Run& run);
44  virtual void analyze(art::Event const&);
45  void endJob();
46 
47  private:
48 
49  // configurables
51  const std::vector< unsigned int > fHumanReadableFEBs;
52  const unsigned int fMinPlaneTrack,fMinCellTrack,fMinCellTrackView; //track cuts
54  typedef enum _cutlevel { kAll, kPass, kNumLevel} Cuts_t;
55 
56  struct ViewH {
57  TH2F* hists[2];
58  };
59  std::map<unsigned int, ViewH> fHBookedChanCP;// map of booked channels in cell-plane, vs group
60  TH2F* fHRawDigits[geo::kXorY][kNumLevel]; // plots of raw digits
61  TH1F* fHExtentPlaneTrk[kNumLevel]; // plots of raw digits
62  TH1F* fHNCellTrk[kNumLevel]; // plots of raw digits
63  TH1F* fHNCellViewTrk[geo::kXorY][kNumLevel]; // plots of raw digits
64 
65  typedef std::map< unsigned int, RawFebPlots* > FebPlotGroup;
66  typedef std::map< unsigned int, FebPlotGroup > FebPlotSet;
67 
68  FebPlotSet fFebRawAll;
69  FebPlotSet fFebRawOnTrack;
70 
71 
72  void FillFebRawPlots(FebPlotSet&, art::Ptr<rawdata::RawDigit> );
74  std::string HRFebToStringLabel(unsigned int,
75  int,
76  unsigned int&, unsigned int&, unsigned int&, unsigned int&);
77 
78  void ResetEventPlots();
79  void FillEventPlots();
80 
81  };
82 
83 
84 }//namespace
85 
86 //////////////////////////////////////////////////////////////////////
87 namespace dprf {
88 
90  EDAnalyzer(pset),
91  fPSet(pset),
92  fHumanReadableFEBs(pset.get< std::vector< unsigned int > >("FEBs")),
93  fMinPlaneTrack(pset.get< unsigned int >("MinPlaneTrack")),
94  fMinCellTrack(pset.get< unsigned int >("MinCellTrack")),
95  fMinCellTrackView(pset.get< unsigned int >("MinCellTrackView")),
96  fRawDigitLabel(pset.get< std::string >("RawDigitLabel")),
97  fCellHitLabel(pset.get< std::string >("CellHitLabel")),
98  fTrackLabel(pset.get< std::string >("TrackLabel"))
99  {
100  }
101 
103  {
104  }
105 
107  {
108 
109  // std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ geo::kXorY = " << "Here1" << std::endl;
110 
111  // if(fHBookedChanCP.empty()) return; // Don't create them all twice
112 
116 
117 
118 
119  for (unsigned int ichan=0; ichan<fHumanReadableFEBs.size(); ichan++) {
120  const int det=geom->DetId();
121 
122  // std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ geo::kXorY = " << "Here1" << std::endl;
123 
124 
125  unsigned int group,diblock,ldcm,lfeb;
127  group,diblock,ldcm,lfeb);
128  const daqchannelmap::dchan chan = cmap->Map()->encodeDChan(det,diblock,ldcm,lfeb,0);
129 
130  bool detailed=true; //need to put in config...
131 
132  (fFebRawAll[group])[chan]=
133  new RawFebPlots(label,"All",detailed,fPSet);
134  (fFebRawOnTrack[group])[chan]=
135  new RawFebPlots(label,"OnTracks",detailed,fPSet);
136 
137  //book and fill histogram showing which channels were booked
138  if (fHBookedChanCP.find(group) == fHBookedChanCP.end())
139  {
140  ViewH hist={{0,0}};
141  for (unsigned int iv=geo::kX; iv<geo::kXorY; iv++)
142  {
143 
145  std::string viewName(view == geo::kX ? "X" : "Y");
146  const std::set<unsigned int>& planes=geom->GetPlanesByView(view);
147  const geo::PlaneGeo* plane=(geom->Plane(*(planes.begin())));
148  //(fHBookedChanCP[group])[iv]=
149  hist.hists[iv]=
150  tfs->make<TH2F>(Form("HBookedChanCP_%s_%d",viewName.c_str(),group),
151  Form("Booked Channel, %s, group %d",viewName.c_str(),group),
152  geom->NPlanes()+1,-.5,geom->NPlanes()+.5,
153  plane->Ncells()+1,-.5,plane->Ncells()+.5);
154  }
155  fHBookedChanCP.insert(std::pair<unsigned int, ViewH>
156  (group,hist));
157  }
158 
159 
160  for (unsigned int ipix=0;
161  ipix<daqchannelmap::NUM_PIXELS; ipix++) {
162  const unsigned int chan = cmap->Map()->encodeDChan(det,diblock,ldcm,lfeb,ipix);
163  const daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
164  const unsigned int plannum = cmap->Map()->getPlane(logchan);
165  const geo::PlaneGeo* plane = geom->Plane(plannum);
166  const unsigned int view=(unsigned int)plane->View();
167 
168  (fHBookedChanCP[group]).hists[view]
169  ->Fill(cmap->Map()->getPlane(logchan),
170  cmap->Map()->getCell(logchan));
171  }//pixels
172 
173  }//booked channels
174 
175  // book general hit histograms
176  std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ geo::kXorY = " << "Here" << std::endl;
177 
178  for (unsigned int iv=geo::kX; iv<geo::kXorY; iv++)
179  {
180  std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ LOOPING $$$$$$$$ "<< std::endl;
181 
183  std::string viewName(view == geo::kX ? "X" : "Y");
184  const std::set<unsigned int>& planes=geom->GetPlanesByView(view);
185  const geo::PlaneGeo* plane=(geom->Plane(*(planes.begin())));
186 
187  for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
188  //for (Cuts_t cut=kAll; cut<kNumLevel; cut++) {
189  std::string cutname(cut==kAll ? "All" : "On Good Tracks");
190  fHRawDigits[iv][cut]=
191  tfs->make<TH2F>(Form("HRawDigits_%s_%d",viewName.c_str(),cut),
192  Form("Raw Digits, %s, %s",viewName.c_str(),cutname.c_str()),
193  geom->NPlanes()+1,-.5,geom->NPlanes()+.5,
194  plane->Ncells()+1,-.5,plane->Ncells()+.5);
195 
196  }//Cut levels
197 
198  }//view
199 
200  // book track cut histograms
201  for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
202  std::string cutname(cut==kAll ? "All Tracks" : "Good Tracks");
204  tfs->make<TH1F>(Form("HExtentPlaneTrk_%d",cut),
205  Form("Track Extent, Planes, %s",cutname.c_str()),
206  geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
207  fHNCellTrk[cut]=
208  tfs->make<TH1F>(Form("HNCellTrk_%d",cut),
209  Form("Track Number of Cells, %s",cutname.c_str()),
210  geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
212  tfs->make<TH1F>(Form("HNCellViewTrk_X_%d",cut),
213  Form("Track Number of Cells, X, %s",cutname.c_str()),
214  geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
216  tfs->make<TH1F>(Form("HNCellViewTrk_Y_%d",cut),
217  Form("Track Number of Cells, Y, %s",cutname.c_str()),
218  geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
219  }// track cuts
220 
221  }
223  {
224  }
225 
226 
228  unsigned int hrfeb, int det,
229  unsigned int &group, unsigned int &diblock, unsigned int &ldcm,
230  unsigned int &lfeb)
231  {
232  // return a useful label from a human-readable FEB id. And return
233  // other useful info
234  group= (hrfeb/1000000)%100;
235  diblock=(hrfeb/ 10000)%100;
236  ldcm= (hrfeb/ 100)%100;
237  lfeb= (hrfeb/ 1)%100;
238 
239  std::string label(Form("g%d_%d_%2.2d_%2.2d_%2.2d",group,det,diblock,ldcm,lfeb));
240  return label;
241 
242  }
243 
245  {
246 
247  ResetEventPlots();
248 
249  // look for RawDigits. Don't continue if none
251  try {
252  evt.getByLabel(fRawDigitLabel,digs);
253  }
254  catch (...) {
255  mf::LogWarning("CPLTWarn") << "Failed to find a RawDigit, bail";
256  return ;
257  }
258 
261 
262  // make "all hits" plots
263  for (unsigned int id=0; id<digs->size(); id++)
264  {
265 
266  art::Ptr<rawdata::RawDigit> dig(digs,id);
267  unsigned int chan=dig->DaqChannel();
268 
269  // fill general plots
270  // need to use "offline" view, not online....
271  daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
272  const unsigned int plannum = cmap->Map()->getPlane(logchan);
273  const geo::PlaneGeo* plane = geom->Plane(plannum);
274  const unsigned int view = (unsigned int)plane->View();
275 
276 
277  fHRawDigits[view][kAll]->Fill(plannum, cmap->Map()->getCell(logchan));
278 
280 
281 
282  }// digit id
283 
284  // look for Tracks. OK to return w/out error if none
286  try {
287  evt.getByLabel(fTrackLabel,tracks);
288  }
289  catch (...) {
290  FillEventPlots();
291  return ;
292  }
293 
294  // can it exist but have zero size?
295  if (!tracks->size()) {
296  FillEventPlots();
297  return;
298  }
299 
300  // get the CellHits, so we can find RawHits on tracks
302  try {
303  evt.getByLabel(fCellHitLabel,cells);
304  }
305  catch (...) {
306  mf::LogError("CPLTWarn") << "\nFound Tracks but no CellHits!\n";
307  assert(0);
308  }
309 
310  for(unsigned int tr = 0; tr < tracks->size(); tr++)
311  {
312  art::Ptr<rb::Track> track(tracks,tr);
313 
314  // does the track have a decent number of hits and planes in each view?
315 
316  if (TrackOK(track)) {
317  // loop over CellHits on tracks
318  for (unsigned int ic=0; ic<track->NCell(); ic++)
319  {
320  art::Ptr<rb::CellHit> cell = track->Cell(ic);
321  const unsigned int chan = cell->DaqChannel();
322  daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
323  const unsigned int plannum = cmap->Map()->getPlane(logchan);
324  const geo::PlaneGeo* plane = geom->Plane(plannum);
325  unsigned int view=(unsigned int)plane->View();
326 
327  fHRawDigits[view][kPass]->Fill(cmap->Map()->getPlane(logchan),
328  cmap->Map()->getCell(logchan));
329 
331  }//cell
332  } // track is OK
333 
334  } // tracks
335 
336  FillEventPlots();
337 
338  } //ana
339 
342  {
343 
344  //Fill all FebRawPlots for this channel, in all groups
347 
348  unsigned int chan=dig->DaqChannel();
349 
350  //key for the FEB is pixel 0 dchan for that FEB.
351  unsigned int lookupchan=chan &
353 
354  //loop over groups, and fill histograms if this channel is booked
355  for (FebPlotSet::iterator iter=plots.begin(); iter!=plots.end(); iter++)
356  {
357 
358  //unsigned int group=iter->first;
359  if (iter->second.find(lookupchan) != iter->second.end()) {
360  iter->second[lookupchan]->AddHit(dig->ADC(),
361  cmap->Map()->getPixel(chan));
362  }
363  } // groups
364  } // FillFebRawPlots
365 
367  {
368  //some histos to monitor the cuts
369  fHExtentPlaneTrk[kAll]->Fill(track->ExtentPlane());
370  fHNCellTrk[kAll]->Fill(track->NCell());
371  fHNCellViewTrk[geo::kX][kAll]->Fill(track->NXCell());
372  fHNCellViewTrk[geo::kY][kAll]->Fill(track->NYCell());
373 
374  // enough planes (extent)
375  if (track->ExtentPlane() < fMinPlaneTrack) return false;
376 
377  // enough cells
378  if (track->NCell() < fMinCellTrack) return false;
379 
380  // enough cells in each view
381  if (track->NXCell() < fMinCellTrackView ||
382  track->NYCell() < fMinCellTrackView ) return false;
383 
384  fHExtentPlaneTrk[kPass]->Fill(track->ExtentPlane());
385  fHNCellTrk[kPass]->Fill(track->NCell());
386  fHNCellViewTrk[geo::kX][kPass]->Fill(track->NXCell());
387  fHNCellViewTrk[geo::kY][kPass]->Fill(track->NYCell());
388  return true;
389  }
391  {
392  // Fill the "per-event" plots
393  FebPlotSet::iterator siter=fFebRawAll.begin();
394  for (;siter!=fFebRawAll.end();siter++) {
395  FebPlotGroup::iterator giter=(siter->second).begin();
396  for (;giter!=(siter->second).end();giter++)
397  (giter->second)->FinishEvent();
398  }
399 
400  siter=fFebRawOnTrack.begin();
401  for (;siter!=fFebRawOnTrack.end();siter++) {
402  FebPlotGroup::iterator giter=(siter->second).begin();
403  for (;giter!=(siter->second).end();giter++)
404  (giter->second)->FinishEvent();
405  }
406  }
408  {
409  // Reset the "per-event" counters
410  FebPlotSet::iterator siter=fFebRawAll.begin();
411  for (;siter!=fFebRawAll.end();siter++) {
412  FebPlotGroup::iterator giter=(siter->second).begin();
413  for (;giter!=(siter->second).end();giter++)
414  (giter->second)->Reset();
415  }
416 
417  siter=fFebRawOnTrack.begin();
418  for (;siter!=fFebRawOnTrack.end();siter++) {
419  FebPlotGroup::iterator giter=(siter->second).begin();
420  for (;giter!=(siter->second).end();giter++)
421  (giter->second)->Reset();
422  }
423  }
424 } //namespace dprf
425 
426 ///////////////////////////////////////////////////////////////
427 namespace dprf
428 {
430 }
void beginRun(const art::Run &run)
const std::string fRawDigitLabel
diblock
print "ROW IS " print row
Definition: geo2elec.py:31
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::map< unsigned int, FebPlotGroup > FebPlotSet
std::map< unsigned int, RawFebPlots * > FebPlotGroup
APD Pixel Number (8bit), valid range 0-31.
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
Definition: event.h:19
TH1F * fHNCellTrk[kNumLevel]
const daqchannelmap::DAQChannelMap * Map() const
Definition: CMap.h:57
const unsigned int fMinCellTrack
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
TString hists[nhists]
Definition: bdt_com.C:3
DEFINE_ART_MODULE(TestTMapFile)
cell_t getCell(lchan logicalchan) const
Decode the cell number from an lchan.
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
const char * label
const std::vector< unsigned int > fHumanReadableFEBs
std::string HRFebToStringLabel(unsigned int, int, unsigned int &, unsigned int &, unsigned int &, unsigned int &)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
uint32_t DaqChannel() const
Definition: RawDigit.h:85
const std::string fTrackLabel
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
const std::vector< Plot > plots
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
enum dprf::ChannelPlots::_cutlevel Cuts_t
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
virtual void analyze(art::Event const &)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
pixel_t getPixel(dchan daqchan) const
Decode the pixel id from a dchan.
TH2F * fHRawDigits[geo::kXorY][kNumLevel]
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
std::map< unsigned int, ViewH > fHBookedChanCP
T * make(ARGS...args) const
const unsigned int fMinCellTrackView
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
void FillFebRawPlots(FebPlotSet &, art::Ptr< rawdata::RawDigit >)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
const unsigned int fMinPlaneTrack
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
const fhicl::ParameterSet & fPSet
TH1F * fHNCellViewTrk[geo::kXorY][kNumLevel]
plane_t getPlane(lchan logicalchan) const
Decode the plane number from an lchan.
assert(nhit_max >=nhit_nbins)
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
ChannelPlots(fhicl::ParameterSet const &p)
unsigned int NPlanes() const
const std::string fCellHitLabel
uint32_t dchan
< DAQ Channel Map Package
Float_t track
Definition: plot.C:35
Definition: fwd.h:28
bool TrackOK(art::Ptr< rb::Track >)
Encapsulate the geometry of one entire detector (near, far, ndos)
TH1F * fHExtentPlaneTrk[kNumLevel]