ValidationPlotter.cpp
Go to the documentation of this file.
1 //File: ValidationPlotter.cpp
2 //Brief: An "algorithm" class for creating online monitoring plots in a way that works
3 // either with or without LArSoft. Should start out by plotting diagnostics from
4 // perl plotting code, then branch out as we identify more pathologies.
5 //Author: Andrew Olivier aolivier@ur.rochester.edu
6 
7 #ifndef CRT_ONLINEPLOTTER_CPP
8 #define CRT_ONLINEPLOTTER_CPP
9 
10 //crt-core includes
12 
13 //crt-alg includes
15 
16 //ROOT includes
17 #include "TDirectory.h"
18 #include "TH1D.h"
19 
20 //c++ includes
21 #include <map>
22 #include <string>
23 #include <algorithm> //For std::max_element()
24 
25 //TODO: Remove me
26 #include <iostream>
27 
28 namespace CRT
29 {
30  template <class TFS>
31  //TFS is an object that defines a function like:
32  //
33  //template <class HIST, class ARGS...> HIST* make(ARGS... args)
34  //
35  //Looks like ART's TFileService.
37  {
38  public:
39 
40  ValidationPlotter(TFS& tfs, std::unique_ptr<CRT::Geometry>&& geom): fFileService(tfs),
41  fGeom(std::move(geom))
42  {
43  }
44 
45  virtual ~ValidationPlotter() = default; //All pointers kept by this class are owned
46  //by another object.
47 
48 
49  void AnalyzeEvent(const std::vector<CRT::Trigger>& triggers) //Make plots from CRT::Triggers
50  {
51  std::cout << "New event\n";
52  //First, fill a nested map to ADC values to find overlaps.
53  CRT::geoMap<int16_t> stripToHits;
54 
55  for(const auto& trigger: triggers)
56  {
57  std::cout << "On module " << trigger.Channel() << ":\n";
58  const auto& hits = trigger.Hits();
59  for(const auto& hit: hits)
60  {
61  std::cout << "\tGot a hit on strip " << hit.Channel() << "\n";
62  const auto id = fGeom->StripID(trigger.Channel(), hit.Channel());
63 
64  //Create a place to put plots from this (strip, module) pair in case it will be needed later.
65  auto foundModule = fModuleToDir.find(id);
66  if(foundModule == fModuleToDir.end()) fModuleToDir.emplace(id, fFileService->mkdir("Module"+std::to_string(trigger.Channel())));
67 
68  auto foundStrip = fStripToChannel.find(id);
69  if(foundStrip == fStripToChannel.end()) fStripToChannel.emplace(id, hit.Channel());
70 
71  if(hit.ADC() > stripToHits[id]) stripToHits[id] = hit.ADC(); //TODO: Assuming for now that there is no more
72  // than 1 hit for each strip in an event.
73  // If there were more than one hit, I'd
74  // keep the largest here because I am
75  // looking for the largest ADC on each
76  // plane for matching later.
77  }
78  }
79 
80  //Next, find the highest ADC hit in each layer and look for an overlap with the other layer in its module.
81  //Fill a map with all such overlaps sorted by plane so that I can find 2D pixels made of pairs of highest-energy
82  //hits.
84  for(const auto& framePair: stripToHits)
85  {
86  //If there are 2 planes in this frame with hits, then
87  //there could be 2D overlaps in this frame.
88  if(framePair.second.size() > 1)
89  {
90  for(const auto& planePair: framePair.second)
91  {
92  for(const auto& modulePair: planePair.second)
93  {
94  //There have to be 2 layers with hits in this plane for overlaps to be possible.
95  if(modulePair.second.size() > 1)
96  {
97  auto layerPair = modulePair.second.begin();
98  auto& stripPairs = layerPair->second;
99  const auto firstMax = std::max_element(stripPairs.begin(), stripPairs.end(), [](const auto& first, const auto& second)
100  {
101  return first.second < second.second;
102  });
103 
104  ++layerPair; //Go to the other layer in this plane
105  auto& secondStripPairs = layerPair->second;
106  const auto secondMax = std::max_element(secondStripPairs.begin(), secondStripPairs.end(), [](const auto& first, const auto& second)
107  {
108  return first.second < second.second;
109  });
110 
111  if(firstMax->first.Overlaps(secondMax->first))
112  {
113  const auto& firstID = firstMax->first;
114  stripToOverlappingHits[firstID].push_back(*firstMax);
115 
116  const auto& secondID = secondMax->first;
117  stripToOverlappingHits[secondID].push_back(*secondMax);
118  } //If firstMax->first overlaps secondMax->first
119  } //If more than one layer has hits
120  } //For each module
121  } //For each plane in this frame
122  } //If more than one plane has hits
123  } //For each frame with hits
124 
125  //Finally, plot the ADCs of strips that are part of 2D pixels
126  for(const auto& framePair: stripToOverlappingHits)
127  {
128  if(framePair.second.size() > 1) //If there is more than one plane with overlapping hit pairs within a frame,
129  //that plane has 2D pixels.
130  {
131  for(const auto& planePair: framePair.second)
132  {
133  const std::vector<std::pair<CRT::StripID, int64_t>>& strips = planePair.second;
134  for(const auto& stripPair: strips) //TODO: Why doesn't planePair.second work here?!
135  {
136  //Look for plots for this StripID and create them if necessary. This way,
137  //I'll only get plots for modules that were part of overlapping pairs.
138  const auto& id = stripPair.first;
139  auto foundStrip = fStrips.find(stripPair.first);
140  if(foundStrip == fStrips.end())
141  {
142  foundStrip = fStrips.emplace(stripPair.first, StripPlots(fModuleToDir.at(id), fStripToChannel.at(id))).first;
143  }
144 
145  auto& strip = foundStrip->second;
146  strip.fOverlappingADCs->Fill(stripPair.second);
147  } //For each strip with an overlap in this plane
148  } //For each plane in this frame
149  } //If more than one plane
150  } //For each frame with pairs of overlapping strips
151  } //AnalyzeEvent()
152 
153  private:
154  TFS fFileService; //Handle to tool for writing out files
155  using DIRECTORY = decltype(fFileService->mkdir("null"));
156 
157  //Group plots for a module into one object
158  struct StripPlots
159  {
160  StripPlots(DIRECTORY& dir, const size_t channel)
161  {
162  fOverlappingADCs = dir.template make<TH1D>(("channel"+std::to_string(channel)+"OverlappingADCs").c_str(),
163  "ADC Values for 2D-Matched Hits;ADC;Hits",
164  1024, 0, 4095);
165  }
166 
167  TH1D* fOverlappingADCs; //ADC values for this strip when it overlaps with something
168  };
169 
170  //Plots produced
171  std::map<CRT::StripID, StripPlots> fStrips; //Plots for each strip
172 
173  //Configuation state
174  std::unique_ptr<CRT::Geometry> fGeom; //Handle to geometrical description of CRT
175 
176  //Hack to get around lack of backwards mapping in Geometry interface
177  std::map<CRT::ModuleID, DIRECTORY> fModuleToDir; //DIRECTORY-channel number pair for each StripID processed so far.
178  std::map<CRT::StripID, size_t> fStripToChannel; //Mapping from strip number to electronics channel number within its board
179  };
180 }
181 
182 #endif //CRT_ONLINEPLOTTER_CPP
std::map< CRT::ModuleID, DIRECTORY > fModuleToDir
void AnalyzeEvent(const std::vector< CRT::Trigger > &triggers)
detail::Node< void, uint8_t > FrameID
Definition: CRTID.h:120
STL namespace.
uint8_t channel
Definition: CRTFragment.hh:201
string dir
decltype(fFileService->mkdir("null")) DIRECTORY
std::map< CRT::StripID, StripPlots > fStrips
def move(depos, offset)
Definition: depos.py:107
Detector simulation of raw signals on wires.
virtual ~ValidationPlotter()=default
std::unique_ptr< CRT::Geometry > fGeom
StripPlots(DIRECTORY &dir, const size_t channel)
std::map< CRT::StripID, size_t > fStripToChannel
ValidationPlotter(TFS &tfs, std::unique_ptr< CRT::Geometry > &&geom)
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34