CosmicEff_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////////////////
2 /// \brief Calculate efficiency for cosmic rays to identify bad channels
3 /// \author Christopher Backhouse - bckhouse@caltech.edu
4 /////////////////////////////////////////////////////////////////////////////////////////
5 
6 #include <vector>
7 
8 // ROOT includes
9 #include "TH1D.h"
10 #include "TH2D.h"
11 
12 // Framework includes
18 #include "fhiclcpp/ParameterSet.h"
20 
21 // NOvA includes
23 #include "Geometry/Geometry.h"
25 #include "RecoBase/CellHit.h"
26 #include "RecoBase/Track.h"
27 #include "RecoBase/HitMap.h"
28 
29 namespace chaninfo
30 {
31  /// \brief Calculate efficiency for cosmic rays to identify bad channels
32  class CosmicEff: public art::EDAnalyzer
33  {
34  public:
35  explicit CosmicEff(const fhicl::ParameterSet& pset);
36  virtual ~CosmicEff();
37 
38  virtual void analyze(const art::Event& evt);
39  virtual void reconfigure(const fhicl::ParameterSet& pset);
40 
41  virtual void beginRun(const art::Run& run);
42  virtual void endJob();
43 
44  protected:
46 
47  TH2F *fHits, *fMiss;
52  };
53 
54  //......................................................................
55  // HitEfficiency::HitEfficiency(fhicl::ParameterSet const & p)
56  // : EDAnalyzer(p)
58  : EDAnalyzer(pset), fHits(0) // This one indicates all the others are uninitialized too
59  {
60  reconfigure(pset);
61  }
62 
63  //......................................................................
65  {
66  }
67 
68  //......................................................................
70  {
71  fTrackLabel = pset.get<std::string>("TrackLabel");
72  }
73 
74  //......................................................................
75  void CosmicEff::beginRun(const art::Run& /*run*/)
76  {
77  if(fHits) return; // Only make the plots once
78 
81 
82  // Make histograms big enough to hold all detector cells
83  const unsigned int maxCell = std::max(geom->Plane(0)->Ncells(), geom->Plane(1)->Ncells());
84  fHits = tfs->make<TH2F>("hits", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
85  fMiss = tfs->make<TH2F>("miss", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
86 
87  fHitsXYLoose = tfs->make<TH2F>("hitsXYLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
88  fMissXYLoose = tfs->make<TH2F>("missXYLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
89 
90  fHitsXYTight = tfs->make<TH2F>("hitsXYTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
91  fMissXYTight = tfs->make<TH2F>("missXYTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
92 
93  fHitsZLoose = tfs->make<TH2F>("hitsZLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
94  fMissZLoose = tfs->make<TH2F>("missZLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
95 
96  fHitsZTight = tfs->make<TH2F>("hitsZTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
97  fMissZTight = tfs->make<TH2F>("missZTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
98  }
99 
100  //......................................................................
102  {
104 
106  evt.getByLabel(fTrackLabel, trks);
107  if(trks->empty()) return;
108 
109  const unsigned int trkMax = trks->size();
110 
111  for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
112  art::Ptr<rb::Track> track(trks, trkIdx);
113 
114  rb::HitMap hmap(track.get());
115 
116  // Which cells does this track's trajectory pass through?
117  std::vector<geo::OfflineChan> xhits, yhits;
118  geom->CountCellsOnLineFast(track->Start(), track->Stop(), xhits, yhits);
119 
120  // If there's enough hits in this view, then count how many tracks were
121  // hit when they should be and how many were missed.
122  if(track->NXCell() > 2){
123  for(unsigned int n = 0; n < xhits.size(); ++n){
124  const unsigned int plane = xhits[n].Plane();
125  const unsigned int cell = xhits[n].Cell();
126  if(hmap.CellExists(plane, cell)){
127  fHits->Fill(plane, cell);
128  }
129  else{
130  fMiss->Fill(plane, cell);
131  }
132  }
133  }
134  if(track->NYCell() > 2){
135  for(unsigned int n = 0; n < yhits.size(); ++n){
136  const unsigned int plane = yhits[n].Plane();
137  const unsigned int cell = yhits[n].Cell();
138  if(hmap.CellExists(plane, cell)){
139  fHits->Fill(plane, cell);
140  }
141  else{
142  fMiss->Fill(plane, cell);
143  }
144  }
145  }
146 
147  // Min/Max cell number of hits on this plane
148  std::map<int, int> minsXY, maxsXY;
149  // Channels that have a hit below/above them
150  std::set<geo::OfflineChan> belows, aboves;
151  // Min/max plane number by view and cell index
152  std::map<int, int> minsZ[2], maxsZ[2];
153  // Channels that have a hit upstream/downstream of them
154  std::set<geo::OfflineChan> lefts, rights;
155 
156  for(unsigned int chitIdx = 0; chitIdx < track->NCell(); ++chitIdx){
157  const art::Ptr<rb::CellHit>& chit = track->Cell(chitIdx);
158  const int plane = chit->Plane();
159  const int cell = chit->Cell();
160  const geo::View_t view = chit->View();
161 
162  // Have to be more careful on min because value defaults to zero
163  if(minsXY.find(plane) == minsXY.end() || cell < minsXY[plane])
164  minsXY[plane] = cell;
165  if(cell > maxsXY[plane]) maxsXY[plane] = cell;
166 
167  if(minsZ[view].find(cell) == minsZ[view].end() || plane < minsZ[view][cell])
168  minsZ[view][cell] = plane;
169  if(plane > maxsZ[view][cell]) maxsZ[view][cell] = plane;
170 
171  belows.insert(geo::OfflineChan(plane, cell+1));
172  aboves.insert(geo::OfflineChan(plane, cell-1));
173 
174  const int nextplane = geom->NextPlaneInView(plane, +1);
175  lefts.insert(geo::OfflineChan(nextplane, cell));
176  const int prevplane = geom->NextPlaneInView(plane, -1);
177  rights.insert(geo::OfflineChan(prevplane, cell));
178  } // end for chitIdx
179 
180  for(std::map<int, int>::iterator it = minsXY.begin(); it != minsXY.end(); ++it){
181  const int plane = it->first;
182  const int from = it->second;
183  const int to = maxsXY[plane];
184 
185  // Expect a hit if there is a hit above and below somewhere on this plane
186  for(int cell = from+1; cell <= to-1; ++cell){
187  if(hmap.CellExists(plane, cell))
188  fHitsXYLoose->Fill(plane, cell);
189  else
190  fMissXYLoose->Fill(plane, cell);
191  }
192  } // end for it
193 
194  for(int view = 0; view < 2; ++view){
195  for(std::map<int, int>::iterator it = minsZ[view].begin(); it != minsZ[view].end(); ++it){
196  const int cell = it->first;
197  const int from = it->second;
198  const int to = maxsZ[view][cell];
199 
200  if(to == from) continue;
201 
202  // Except a hit if there's one on an earlier and later plane
203  int plane = geom->NextPlaneInView(from, +1);
204  while(plane != to){
205  if(hmap.CellExists(plane, cell))
206  fHitsZLoose->Fill(plane, cell);
207  else
208  fMissZLoose->Fill(plane, cell);
209 
210  plane = geom->NextPlaneInView(plane, +1);
211  }
212  } // end for it
213  } // end for view
214 
215  for(std::set<geo::OfflineChan>::iterator it = belows.begin(); it != belows.end(); ++it){
216  const int plane = it->Plane();
217  const int cell = it->Cell();
218  if(!aboves.count(geo::OfflineChan(plane, cell))) continue;
219 
220  // Expect a hit if both our neighbours are hit
221  if(hmap.CellExists(plane, cell))
222  fHitsXYTight->Fill(plane, cell);
223  else
224  fMissXYTight->Fill(plane, cell);
225  }
226 
227  for(std::set<geo::OfflineChan>::iterator it = lefts.begin(); it != lefts.end(); ++it){
228  const int plane = it->Plane();
229  const int cell = it->Cell();
230  if(!rights.count(geo::OfflineChan(plane, cell))) continue;
231 
232  // Expect a hit if both our neighbouring planes are hit here
233  if(hmap.CellExists(plane, cell))
234  fHitsZTight->Fill(plane, cell);
235  else
236  fMissZTight->Fill(plane, cell);
237  }
238  } // end for trkIdx
239  }
240 
241  //......................................................................
243  {
247 
248  // Store what the occupancy cuts think about the channels
249  const unsigned int maxCell = std::max(geom->Plane(0)->Ncells(), geom->Plane(1)->Ncells());
250  TH2* hKnown = tfs->make<TH2F>("known", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
251  // TH2* hOccs = tfs->make<TH2F>("occs", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
252  for(unsigned int p = 0; p < geom->NPlanes(); ++p){
253  for(unsigned int c = 0; c < geom->Plane(p)->Ncells(); ++c){
254  if(bcl->IsBad(p, c)){
255  hKnown->Fill(p, c);
256  }
257  // hOccs->Fill(p, c, bcl->GetOccupancy(p, c));
258  }
259  }
260 
261  // Make sums, ratios, and projections of what we already have:
262 
263  TH2F* tot = tfs->make<TH2F>(*fHits);
264  tot->SetName("tot");
265  tot->Add(fMiss);
266 
267  TH2F* totXYLoose = tfs->make<TH2F>(*fHitsXYLoose);
268  totXYLoose->SetName("totXYLoose");
269  totXYLoose->Add(fMissXYLoose);
270 
271  TH2F* totXYTight = tfs->make<TH2F>(*fHitsXYTight);
272  totXYTight->SetName("totXYTight");
273  totXYTight->Add(fMissXYTight);
274 
275  TH2F* totZLoose = tfs->make<TH2F>(*fHitsZLoose);
276  totZLoose->SetName("totZLoose");
277  totZLoose->Add(fMissZLoose);
278 
279  TH2F* totZTight = tfs->make<TH2F>(*fHitsZTight);
280  totZTight->SetName("totZTight");
281  totZTight->Add(fMissZTight);
282 
283  TH2F* frac = tfs->make<TH2F>(*fMiss);
284  frac->SetName("frac");
285  frac->Divide(tot);
286 
287  TH2F* fracXYLoose = tfs->make<TH2F>(*fMissXYLoose);
288  fracXYLoose->SetName("fracXYLoose");
289  fracXYLoose->Divide(totXYLoose);
290 
291  TH2F* fracXYTight = tfs->make<TH2F>(*fMissXYTight);
292  fracXYTight->SetName("fracXYTight");
293  fracXYTight->Divide(totXYTight);
294 
295  TH2F* fracZLoose = tfs->make<TH2F>(*fMissZLoose);
296  fracZLoose->SetName("fracZLoose");
297  fracZLoose->Divide(totZLoose);
298 
299  TH2F* fracZTight = tfs->make<TH2F>(*fMissZTight);
300  fracZTight->SetName("fracZTight");
301  fracZTight->Divide(totZTight);
302 
303  TH1F* proj = tfs->make<TH1F>("proj", "", 101, 0, 1.01);
304  TH1F* projXYLoose = tfs->make<TH1F>("projXYLoose", "", 101, 0, 1.01);
305  TH1F* projXYTight = tfs->make<TH1F>("projXYTight", "", 101, 0, 1.01);
306  TH1F* projZLoose = tfs->make<TH1F>("projZLoose", "", 101, 0, 1.01);
307  TH1F* projZTight = tfs->make<TH1F>("projZTight", "", 101, 0, 1.01);
308 
309  // Projections don't include channels already known to be bad.
310 
311  const int X = frac->GetNbinsX()+1;
312  const int Y = frac->GetNbinsY()+1;
313  for(int x = 1; x < X; ++x){
314  for(int y = 1; y < Y; ++y){
315  const bool isBad = hKnown->GetBinContent(x, y);
316 
317  if(isBad){
318  frac->SetBinContent(x, y, 0);
319  fracXYLoose->SetBinContent(x, y, 0);
320  fracXYTight->SetBinContent(x, y, 0);
321  fracZLoose->SetBinContent(x, y, 0);
322  fracZTight->SetBinContent(x, y, 0);
323  continue;
324  }
325 
326  const double fr = frac->GetBinContent(x, y);
327  if(tot->GetBinContent(x, y)) proj->Fill(fr);
328 
329  const double frLoXY = fracXYLoose->GetBinContent(x, y);
330  if(totXYLoose->GetBinContent(x, y)) projXYLoose->Fill(frLoXY);
331 
332  const double frTiXY = fracXYTight->GetBinContent(x, y);
333  if(totXYTight->GetBinContent(x, y)) projXYTight->Fill(frTiXY);
334 
335  const double frLoZ = fracZLoose->GetBinContent(x, y);
336  if(totZLoose->GetBinContent(x, y)) projZLoose->Fill(frLoZ);
337 
338  const double frTiZ = fracZTight->GetBinContent(x, y);
339  if(totZTight->GetBinContent(x, y)) projZTight->Fill(frTiZ);
340  } // end for Y
341  } // end for X
342  }
343 
345 
346 } // end namespace chaninfo
347 ///////////////////////////////////////////////////////////////////////////
virtual void analyze(const art::Event &evt)
T max(const caf::Proxy< T > &a, T b)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
virtual void reconfigure(const fhicl::ParameterSet &pset)
set< int >::iterator it
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
CosmicEff(const fhicl::ParameterSet &pset)
void CountCellsOnLineFast(double x1, double y1, double z1, double x2, double y2, double z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
Make list of cells in each view traversed by a line segment.
Definition: Run.h:31
Float_t Y
Definition: plot.C:38
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
unsigned short Cell() const
Definition: CellHit.h:40
T get(std::string const &key) const
Definition: ParameterSet.h:231
double frac(double x)
Fractional part.
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Definition: run.py:1
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
Calculate efficiency for cosmic rays to identify bad channels.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
A (plane, cell) pair.
Definition: OfflineChan.h:17
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
Float_t proj
Definition: plot.C:35
void geom(int which=0)
Definition: geom.C:163
Interface to the run-by-run list of bad channels.
Definition: BadChanList.h:31
virtual void beginRun(const art::Run &run)
unsigned int NPlanes() const
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t X
Definition: plot.C:38
bool IsBad(int plane, int cell)
Float_t track
Definition: plot.C:35
Simple object representing a (plane, cell) pair.
Encapsulate the geometry of one entire detector (near, far, ndos)