AdcDetectorPlotter_tool.cc
Go to the documentation of this file.
1 // AdcDetectorPlotter_tool.cc
2 
3 #include "AdcDetectorPlotter.h"
4 #include <iostream>
5 #include <sstream>
14 #include "TH2F.h"
15 #include "TCanvas.h"
16 #include "TColor.h"
17 #include "TStyle.h"
18 #include "TDirectory.h"
19 #include "TFile.h"
20 #include "TGraph.h"
21 
22 using std::string;
23 using std::cout;
24 using std::cin;
25 using std::endl;
26 using std::ostringstream;
27 
28 using Tick = AdcSignalVector::size_type;
29 using geo::GeometryCore;
30 using std::vector;
31 
32 //**********************************************************************
33 // Local methods.
34 //**********************************************************************
35 
36 namespace {
37 
38 void initializeState(AdcDetectorPlotter::State& state, const AdcChannelData& acd) {
39  state.reportCount = 0;
40  state.channelCount = 0;
41  state.run = acd.run();
42  state.subrun = acd.subRun();
43  state.event = acd.event();
44  state.ppad.reset(nullptr);
45 }
46 
47 } // end unnamed namespace
48 
49 //**********************************************************************
50 // Class methods.
51 //**********************************************************************
52 
54 : m_LogLevel(ps.get<int>("LogLevel")),
55  m_WireAngle(ps.get<float>("WireAngle")),
56  m_DataType(ps.get<int>("DataType")),
57  m_Tick0(ps.get<float>("Tick0")),
58  m_DriftSpeed(ps.get<float>("DriftSpeed")),
59  m_XMin(ps.get<float>("XMin")),
60  m_XMax(ps.get<float>("XMax")),
61  m_ZMin(ps.get<float>("ZMin")),
62  m_ZMax(ps.get<float>("ZMax")),
63  m_SignalThreshold(ps.get<float>("SignalThreshold")),
64  m_SkipBadChannels(ps.get<bool>("SkipBadChannels")),
65  m_ShowAllTicks(ps.get<bool>("ShowAllTicks")),
66  m_FirstTick(ps.get<unsigned long>("FirstTick")),
67  m_LastTick(ps.get<unsigned long>("LastTick")),
68  m_ShowWires(ps.get<bool>("ShowWires")),
69  m_ShowCathode(ps.get<bool>("ShowCathode")),
70  m_ShowTpcSets(ps.get<IndexVector>("ShowTpcSets")),
71  m_ShowGrid(ps.get<bool>("ShowGrid")),
72  m_Title(ps.get<string>("Title")),
73  m_PlotTitle(ps.get<string>("PlotTitle")),
74  m_FileName(ps.get<string>("FileName")),
75  m_pChannelStatusProvider(nullptr),
76  m_state(new State) {
77  const string myname = "AdcDetectorPlotter::ctor: ";
79  string stringBuilder = "adcStringBuilder";
80  m_adcStringBuilder = ptm->getShared<AdcChannelStringTool>(stringBuilder);
81  if ( m_adcStringBuilder == nullptr ) {
82  cout << myname << "WARNING: AdcChannelStringTool not found: " << stringBuilder << endl;
83  }
84  if ( m_SkipBadChannels ) {
85  if ( m_LogLevel >= 1 ) cout << myname << "Fetching channel status service." << endl;
87  if ( m_pChannelStatusProvider == nullptr ) {
88  cout << myname << "WARNING: Channel status provider not found." << endl;
89  }
90  }
91  if ( m_LogLevel ) {
92  cout << myname << "Configuration: " << endl;
93  cout << myname << " LogLevel: " << m_LogLevel << endl;
94  cout << myname << " WireAngle: " << m_WireAngle << endl;
95  cout << myname << " DataType: " << m_DataType << endl;
96  cout << myname << " Tick0: " << m_Tick0 << endl;
97  cout << myname << " DriftSpeed: " << m_DriftSpeed << " cm/tick" << endl;
98  cout << myname << " XMin: " << m_XMin << " cm" << endl;
99  cout << myname << " XMax: " << m_XMax << " cm" << endl;
100  cout << myname << " ZMin: " << m_ZMin << " cm" << endl;
101  cout << myname << " ZMax: " << m_ZMax << " cm" << endl;
102  cout << myname << " SignalThreshold: " << m_SignalThreshold << endl;
103  cout << myname << " SkipBadChannels: " << (m_SkipBadChannels ? "true" : "false") << endl;
104  cout << myname << " ShowAllTicks: " << m_ShowAllTicks << endl;
105  cout << myname << " FirstTick: " << m_FirstTick << endl;
106  cout << myname << " LastTick: " << m_LastTick << endl;
107  cout << myname << " ShowWires: " << m_ShowWires << endl;
108  cout << myname << " ShowCathode: " << m_ShowCathode << endl;
109  cout << myname << " ShowTpcSets: [";
110  bool first = true;
111  for ( Index itps : m_ShowTpcSets ) {
112  if ( first ) first = false;
113  else cout << ", ";
114  cout << itps;
115  }
116  cout << "]" << endl;
117  cout << myname << " ShowGrid: " << m_ShowGrid << endl;
118  cout << myname << " Title: " << m_Title << endl;
119  cout << myname << " PlotTitle: " << m_PlotTitle << endl;
120  cout << myname << " FileName: " << m_FileName << endl;
121  }
122  WireSelector& sel = getState()->sel;
125  //for ( Index itps : m_ShowTpcSets ) sel.selectTpcSet(itps);
126  const WireSelector::WireInfoVector& wdat = sel.fillData();
127  const WireSelector::WireInfoMap& wmap = sel.fillDataMap();
128  const WireSelector::WireSummary& wsum = sel.fillWireSummary();
129  const GeometryCore* pgeo = sel.geometry();
130  string gname = pgeo == nullptr ? "NONE" : pgeo->DetectorName();
131  if ( m_LogLevel ) {
132  cout << myname << "Derived: " << endl;
133  cout << myname << " Geometry name: " << gname << endl;
134  cout << myname << " # selected planes: " << sel.planeIDs().size() << endl;
135  cout << myname << " # selected wires: " << wdat.size() << endl;
136  cout << myname << " # selected wire points: " << wsum.size() << endl;
137  cout << myname << " # selected channels: " << wmap.size() << endl;
138  }
139 }
140 
141 //**********************************************************************
142 
144 }
145 
146 //**********************************************************************
147 
149  const string myname = "AdcDetectorPlotter::view: ";
150  DataMap ret;
151  State& state = *getState();
152  if ( m_LogLevel >= 2 ) cout << myname << "Begin call " << state.reportCount
153  << "/" << state.jobCount << "." << endl;
154  if ( acds.size() == 0 ) {
155  cout << myname << "WARNING: Channel map is empty. No data extracted." << endl;
156  return ret.setStatus(1);
157  }
158  const AdcChannelData& acdFirst = acds.begin()->second;
159  const AdcChannelData& acdLast = acds.rbegin()->second;
160  string hname = "hdet";
161  int npadx = 1200;
162  int npady = 1000;
163  double xmin = m_XMin;
164  double xmax = m_XMax;
165  string sttlx = "Drift coordinate [cm]";
166  double xsign = 1.0;
167  if ( xmax < xmin ) {
168  xmin = m_XMax;
169  xmax = m_XMin;
170  sttlx = "-" + sttlx;
171  xsign = -1.0;
172  }
173  string sttly = "Wire coordinate [cm]";
174  if ( state.reportCount ) {
175  if ( acdFirst.run() != state.run ||
176  acdFirst.subRun() != state.subrun ||
177  acdFirst.event() != state.event ) {
178  cout << myname << "ERROR: Received unexpected event ID. Clearing data." << endl;
179  cout << myname << "State: " << state.event << "-" << state.subrun << "-" << state.event << endl;
180  cout << myname << " Data: " << acdFirst.event() << "-" << acdFirst.subRun() << "-" << acdFirst.event() << endl;
181  state.reportCount = 0;
182  }
183  }
184  if ( state.reportCount == 0 ) {
185  if ( m_LogLevel >= 2 ) cout << myname << " Starting new event." << endl;
186  initializeState(state, acdFirst);
187  string sttl = AdcChannelStringTool::build(m_adcStringBuilder, acdFirst, m_Title);
190  // Create graph.
191  state.ppad.reset(new TPadManipulator(npadx, npady));
192  TGraph* pg = new TGraph;
193  pg->SetTitle(sttl.c_str());
194  pg->GetXaxis()->SetTitle(sttlx.c_str());
195  pg->GetYaxis()->SetTitle(sttly.c_str());
196  state.ppad->add(pg, "P");
197  state.ppad->graph()->Expand(512); // Allocate 1024 points.
198  state.ppad->setRangeX(xmin, xmax);
199  state.ppad->setRangeY(m_ZMin, m_ZMax);
200  if ( m_ShowGrid ) state.ppad->setGrid();
201  LineColors cols;
202  const WireSelector::WireInfoVector& wins = state.sel.fillData();
203  if ( m_ShowWires ) {
204  TGraph* pgw = new TGraph;
205  pgw->SetMarkerColor(cols.red());
206  pgw->Expand(wins.size());
207  for ( Index iwin=0; iwin<wins.size(); ++iwin ) {
208  const WireSelector::WireInfo& win = wins[iwin];
209  pgw->SetPoint(iwin, xsign*win.x, win.z);
210  }
211  state.ppad->add(pgw, "P");
212  }
213  if ( m_ShowCathode ) {
214  TGraph* pgc = new TGraph;
215  pgc->SetMarkerColor(cols.green());
216  pgc->Expand(wins.size());
217  for ( Index iwin=0; iwin<wins.size(); ++iwin ) {
218  const WireSelector::WireInfo& win = wins[iwin];
219  pgc->SetPoint(iwin, xsign*(win.x + win.driftMax), win.z);
220  }
221  state.ppad->add(pgc, "P");
222  }
223  // Add lower left label.
224  if ( spttl.size() ) {
225  state.pttl.reset(new TLatex(0.01, 0.015, spttl.c_str()));
226  state.pttl->SetNDC();
227  state.pttl->SetTextFont(42);
228  state.pttl->SetTextSize(0.030);
229  state.ppad->add(state.pttl.get());
230  }
231  ++state.jobCount;
232  } else {
233  TGraph* pgr = state.ppad->graph();
234  if ( pgr == nullptr ) {
235  cout << "ERROR: Graph not found." << endl;
236  return ret;
237  }
238  if ( m_LogLevel >= 2 ) {
239  cout << myname << " Adding to existing event. Graph point count is " << pgr->GetN() << endl;
240  }
241  }
242  ++state.reportCount;
243  Tick maxtick = 0;
244  for ( const AdcChannelDataMap::value_type& iacd : acds ) {
245  if ( iacd.first == AdcChannelData::badChannel() ) {
246  cout << myname << "WARNING: Channel map has invalid channels. No plot is created." << endl;
247  }
248  Tick ntick = iacd.second.samples.size();
249  if ( ntick > maxtick ) maxtick = ntick;
250  }
251  AdcIndex chanFirst = acdFirst.channel();
252  AdcIndex chanLast = acdLast.channel();
253  AdcIndex nchan = chanLast + 1 - chanFirst;
254  if ( m_LogLevel >= 2 ) cout << myname << " Input channel count is " << nchan << endl;
255  // Fill graph.
256  for ( const AdcChannelDataMap::value_type& iacd : acds ) {
257  if ( m_LogLevel >= 3 ) cout << myname << " Filling with channel " << iacd.first << endl;
258  const AdcChannelData& acd = iacd.second;
259  if ( m_SkipBadChannels && m_pChannelStatusProvider != nullptr &&
261  if ( m_LogLevel >= 3 ) cout << myname << " Skipping bad channel " << acd.channel() << endl;
262  } else {
263  if ( m_LogLevel >= 3 ) cout << myname << " Adding channel " << acd.channel() << endl;
264  addChannel(acd, xsign);
265  }
266  }
267  if ( state.ppad->graph()->GetN() == 0 ) {
268  cout << myname << "Graph has no points. Adding one to avoid root exception." << endl;
269  state.ppad->graph()->SetPoint(0, m_XMin, m_ZMin);
270  }
271  if ( m_LogLevel >= 2 ) cout << myname << " Graph point count: " << state.ppad->graph()->GetN() << endl;
272  return ret;
273 }
274 
275 //**********************************************************************
276 
278  const string myname = "AdcDetectorPlotter::endEvent: ";
279  DataMap ret;
280  State& state = *getState();
281  if ( state.reportCount == 0 ) {
282  cout << myname << "WARNING: No data recorded." << endl;
283  return ret;
284  }
285  if ( evi.run != state.run ) {
286  cout << myname << "ERROR: Received request for unexpected run: " << evi.run << endl;
287  return ret.setStatus(1);
288  }
289  if ( evi.event != state.event ) {
290  cout << myname << "ERROR: Received request for unexpected event: " << evi.event << endl;
291  return ret.setStatus(2);
292  }
293  state.ppad->print(state.ofname);
294  state.reportCount = 0;
295  return ret;
296 }
297 
298 //**********************************************************************
299 
300 int AdcDetectorPlotter::addChannel(const AdcChannelData& acd, double xsign) const {
301  const string myname = "AdcDetectorPlotter::addChannel: ";
302  bool isRaw = m_DataType == 1;
303  bool isPrep = m_DataType == 0;
304  if ( !isRaw && !isPrep ) {
305  cout << myname << "ERROR: Invalid data type: " << m_DataType << endl;
306  return 2;
307  }
308  State& state = *getState();
309  TGraph* pg = state.ppad->graph();
310  if ( pg == nullptr ) {
311  cout << myname << "ERROR: Graph is missing." << endl;
312  return 1;
313  }
314  AdcChannel icha = acd.channel();
315  auto rng = state.sel.dataMap().equal_range(icha);
316  Index nsam = isRaw ? acd.raw.size() : acd.samples.size();
317  if ( m_LogLevel >= 4 ) {
318  cout << myname << " Adding " << nsam << (isRaw ? " raw" : " processed")
319  << " sample" << (nsam==1 ? "" : "s") << endl;
320  }
321  for ( auto ient=rng.first; ient!=rng.second; ++ient) {
322  const WireSelector::WireInfo& win = *(ient->second);
323  float z = win.z;
324  float x1 = win.x1();
325  float x2 = win.x2();
326  float driftVelocity = win.driftSign()*m_DriftSpeed;
327  Index isam1 = 0;
328  Index isam2 = nsam;
329  if ( m_LastTick > m_FirstTick ) {
330  if ( m_FirstTick > isam1 ) isam1 = m_FirstTick;
331  if ( m_LastTick < isam2 ) isam2 = m_LastTick;
332  }
333  for ( Index isam=isam1; isam<isam2; ++isam ) {
334  float sig = isRaw ? acd.raw[isam] - acd.pedestal : acd.samples[isam];
335  if ( sig > m_SignalThreshold ) {
336  float x = win.x + driftVelocity*(isam - m_Tick0);
337  if ( ! m_ShowAllTicks ) {
338  if ( x < x1 ) continue;
339  if ( x > x2 ) continue;
340  }
341  Index ipt = pg->GetN();
342  pg->SetPoint(ipt, xsign*x, z);
343  if ( m_LogLevel >= 4 ) {
344  ostringstream sout;
345  sout.precision(2);
346  sout << myname << "Added point " << ipt << " (" << x << ", " << z << ")";
347  cout << sout.str() << endl;
348  }
349  } else if ( m_LogLevel >= 5 ) {
350  cout << myname << "Skipped sample " << isam << " with signal " << sig << endl;
351  }
352  }
353  }
354 /*
355  AdcChannel chan = acd.channel();
356  const AdcSignalVector& sams = acd.samples;
357  const AdcCountVector& raw = acd.raw;
358  AdcSignal ped = 0.0;
359  bool isRawPed = false;
360  if ( isRaw ) {
361  ped = acd.pedestal;
362  isRawPed = ped != AdcChannelData::badSignal;
363  }
364  unsigned int ibin = ph->GetBin(1, chan-chanFirst+1);
365  for ( Tick isam=0; isam<sams.size(); ++isam, ++ibin ) {
366  unsigned int isig = isam + m_FirstTick;
367  if ( isPrep ) {
368  if ( isig < sams.size() ) ph->SetBinContent(ibin, sams[isig]);
369  } else if ( isRawPed ) {
370  if ( isig < raw.size() ) ph->SetBinContent(ibin, raw[isig] - ped);
371  } else {
372  cout << myname << "Fill failed for bin " << ibin << endl;
373  }
374  }
375 */
376  return 0;
377 }
378 
379 //**********************************************************************
380 
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
const WireInfoMap & dataMap() const
Definition: WireSelector.h:128
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
float x1() const
Definition: WireSelector.h:59
AdcIndex subRun() const
unsigned int Index
const lariov::ChannelStatusProvider * m_pChannelStatusProvider
DataMap & setStatus(int stat)
Definition: DataMap.h:130
std::string string
Definition: nybbler.cc:12
const PlaneIDVector & planeIDs() const
Definition: WireSelector.h:123
static ColorType red()
Definition: LineColors.h:27
struct vector vector
DataMap viewMap(const AdcChannelDataMap &acds) const override
const GeometryCore * geometry() const
Definition: WireSelector.h:106
art framework interface to geometry description
static ColorType green()
Definition: LineColors.h:28
StatePtr getState() const
DataMap endEvent(const DuneEventInfo &) const override
int addChannel(const AdcChannelData &acd, double xfac) const
float x2() const
Definition: WireSelector.h:60
void selectWireAngle(double wireAngle, double tol=0.001)
void selectTpcSets(const IndexVector &itpss)
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
const AdcChannelStringTool * m_adcStringBuilder
static std::string build(const AdcChannelStringTool *ptool, const AdcChannelData &acd, const DataMap &dm, std::string spat)
std::vector< WireInfo > WireInfoVector
Definition: WireSelector.h:97
AdcIndex run() const
AdcIndex event() const
AdcSignalVector::size_type Tick
AdcDetectorPlotter(fhicl::ParameterSet const &ps)
static constexpr double ps
Definition: Units.h:99
AdcCountVector raw
Description of geometry of one entire detector.
unsigned int AdcIndex
Definition: AdcTypes.h:15
const WireSummary & fillWireSummary()
std::vector< Index > IndexVector
Channel channel() const
AdcSignal pedestal
const WireInfoVector & fillData()
float driftSign() const
Definition: WireSelector.h:65
unsigned int AdcChannel
Definition: AdcTypes.h:50
std::multimap< Index, const WireInfo * > WireInfoMap
Definition: WireSelector.h:98
Interface for experiment-specific channel quality info provider.
list x
Definition: train.py:276
std::map< AdcChannel, AdcChannelData > AdcChannelDataMap
static Index badChannel()
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
Interface for experiment-specific service for channel quality info.
const WireInfoMap & fillDataMap()
int bool
Definition: qglobal.h:345
static DuneToolManager * instance(std::string fclname="", int dbg=1)
AdcSignalVector samples
T * getShared(std::string name)
QTextStream & endl(QTextStream &s)