CosmicCalib_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CosmicCalib_module.cc
3 // \brief Write out PE vs W plots from the input PCHits files
4 // \author Christopher Backhouse, Dan Pershey - bckhouse@caltech.edu, pershey@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 // ROOT includes
8 #include "TH2.h"
9 
10 // Framework includes
13 #include "art_root_io/TFileService.h"
18 #include "cetlib_except/exception.h"
19 
20 // NOvASoft includes
24 #include "Geometry/Geometry.h"
29 #include "NovaDAQConventions/DAQConventions.h"
30 #include "RawData/RawDigit.h"
31 
32 namespace calib
33 {
34  class CosmicCalib: public art::EDFilter
35  {
36  public:
37  explicit CosmicCalib(const fhicl::ParameterSet& pset);
38  ~CosmicCalib();
39 
40  virtual bool filter(art::Event& evt);
41 
42  void reconfigure(const fhicl::ParameterSet& pset);
43 
44  virtual void beginJob();
45  virtual bool beginRun(art::Run& run);
46  virtual bool endSubRun(art::SubRun& sr);
47  virtual void endJob();
48 
49  protected:
50  caldp::AttenHists GetChannelHists(bool const& isData,
51  int const& plane,
52  int const& cell);
54  bool const& isData,
55  int const& plane,
56  int const& cell);
57  int GetView(int const& plane);
58  double PCHitPE(caldp::PCHit const& pchit);
59 
60  void FillPathLengthInfo(std::vector<caldp::PCHit> const& pchits,
61  caldp::PathType_t const& pathType,
62  bool const& isRealData,
63  bool const& useBelowThreshold=false);
64 
66 
67  // General histograms for investigating variables.
68  // Some of these could be derived from fChannelMap, but it's easier
69  // just to fill them as we go.
70 
72 
76 
80  int fPECalculation; ///< method of determining the value of a pe
81  ///< 0 = PCHit::PE(), 1 = DCS ADC value from raw digit
82  std::map<geo::CellUniqueId, rawdata::RawDigit> fCellIdToRawDigit; ///< map cell id to the raw digit
84 
85 
86  TH2* fCellPE[4];
87  TH2* fPlanePE[4];
88 
89  // Direction cosines
90  TH1 *fdx, *fdy, *fdz;
91  TH2 *fdxdy, *fdxdz, *fdydz;
92 
93  // Track vertices
96 
97  TH1* fPathLength[4];
98 
99  /// This is where \ref GetChannelHist stores its histograms. Map is from block #
101  std::map<int, std::unique_ptr<caldp::AttenProfilesMap>> fChannelMapProf;
102  /// This is what we put in the histogram file
104  std::map<int, caldp::AttenProfilesMap> fChannelMapProfTotal;
105 
107  };
108 
109  //......................................................................
111  : EDFilter(pset)
112  , fdx(0) // This one is the signal that no histograms are created yet
113  {
114  reconfigure(pset);
115 
116  if(fProfileMode){
117  // TODO: I don't think we can determine the detector in time...
118  for(int block = 0; block < 28; ++block){
119  produces<caldp::AttenProfilesMap, art::InSubRun>(TString::Format("block%d", block).Data());
120  }
121  }
122  else{
123  produces<caldp::AttenHistsMap, art::InRun>();
124  }
125  }
126 
127  //......................................................................
129  {
130  }
131 
132  //......................................................................
133 
135  {
136  fPCHitsLabel = pset.get<std::string>("PCHitsLabel");
137  fIncludeBelowThresholdZeros = pset.get<bool>("IncludeBelowThresholdZeros");
138  fProfileMode = pset.get<bool>("ProfileMode");
139  fConsolidateViewsMode = pset.get<bool>("ConsolidateViewsMode");
140  fHistoOutput = pset.get<bool>("HistoOutput");
141  fPECalculation = pset.get<int >("PECalculation", 0);
142  fRawDataLabel = pset.get<std::string>("RawDataLabel", "daq");
143  }
144 
145  //......................................................................
147  {
148  }
149 
150  //......................................................................
152  {
155 
156  // The rest is making histograms. Don't do it again on a subsequent run
157  if(fdx) return true;
158 
159  const TString viewStrArr[4]={"X","Y","muX","muY"};
161  bool isND = fGeom->DetId() == novadaq::cnv::kNEARDET;
162  for(int view = geo::kX; view <= (geo::kY+2*isND); ++view){
163  const TString viewStr = viewStrArr[view];
164 
165  double maxW;
167  maxW = 900;
168  }
169  else if(fGeom->DetId() == novadaq::cnv::kNEARDET){
170  maxW = 250;
171  }
172  else{
173  // NDOS (hopefully)
174  maxW = (view == geo::kX) ? 250 : 150;
175  }
176 
178 
179  v.WPE = tfs->make<TH2F>("h_WPE_in"+viewStr,
180  ";W;PE",
181  2*maxW, -maxW, +maxW,
182  1000, 0, 1000);
183 
184  v.WPE_corr = tfs->make<TH2F>("h_WPE_corr_avg_in"+viewStr,
185  ";W;PE/cm",
186  2*maxW, -maxW, +maxW,
187  1000, 0, 200);
188 
189  v.WPE_corr_xy = tfs->make<TH2F>("h_WPE_corr_xy_in"+viewStr,
190  ";W;PE/cm",
191  2*maxW, -maxW, +maxW,
192  1000, 0, 200);
193 
194  v.WPE_corr_xy_truth = tfs->make<TH2F>("h_WPE_corr_xy_truth_in"+viewStr,
195  ";W;MeV/cm",
196  2*maxW, -maxW, +maxW,
197  1000, 0, 10);
198 
199  v.WPE_corr_z = tfs->make<TH2F>("h_WPE_corr_z_in"+viewStr,
200  ";W;PE/cm",
201  2*maxW, -maxW, +maxW,
202  1000, 0, 200);
203  v.WPE_corr_traj = tfs->make<TH2F>("h_WPE_corr_traj_in"+viewStr,
204  ";W;PE/cm",
205  2*maxW, -maxW, +maxW,
206  1000, 0, 200);
207 
208  fSumChannels[view] = v;
209 
210  fCellPE[view] = tfs->make<TH2F>("h_CellPE_in"+viewStr, ";Cell;PE",
211  105, -5, 100, 1000, 0, 1000);
212 
213  fPlanePE[view] = tfs->make<TH2F>("h_PlanePE_in"+viewStr, ";Plane;PE",
214  103, -6, 200, 1000, 0, 1000);
215 
216  fPathLength[view] = tfs->make<TH1F>("h_PathLength"+viewStr,
217  ";Path length (cm)",
218  80, 0, 20);
219  } // end for view
220 
221  fdx = tfs->make<TH1F>("h_dx", ";|dx|", 100, 0, 1);
222  fdy = tfs->make<TH1F>("h_dy", ";|dy|", 100, 0, 1);
223  fdz = tfs->make<TH1F>("h_dz", ";|dz|", 100, 0, 1);
224 
225  fdxdy = tfs->make<TH2F>("h_dxdy", ";|dx|;|dy|", 100, 0, 1, 100, 0, 1);
226  fdxdz = tfs->make<TH2F>("h_dxdz", ";|dx|;|dz|", 100, 0, 1, 100, 0, 1);
227  fdydz = tfs->make<TH2F>("h_dydz", ";|dy|;|dz|", 100, 0, 1, 100, 0, 1);
228 
229  fStartXY = tfs->make<TH2F>("h_StartXY", ";x;y",
230  100, -200, +200, 100, -250, +250);
231  fStartZX = tfs->make<TH2F>("h_StartZX", ";z;x",
232  100, -50, +900, 100, -200, +200);
233  fStartZY = tfs->make<TH2F>("h_StartZY", ";z;y",
234  100, -50, +900, 100, -250, +250);
235 
236  fStopXY = tfs->make<TH2F>("h_StopXY", ";x;y",
237  100, -200, +200, 100, -250, +250);
238  fStopZX = tfs->make<TH2F>("h_StopZX", ";z;x",
239  100, -50, +900, 100, -200, +200);
240  fStopZY = tfs->make<TH2F>("h_StopZY", ";z;y",
241  100, -50, +900, 100, -250, +250);
242 
243  return true;
244  }
245 
246  //......................................................................
248  {
249  if(fProfileMode){
250  for(int block = 0; block < 28; ++block){
251  auto it = fChannelMapProf.find(block);
252  if(it == fChannelMapProf.end()) continue;
253 
254  if(fHistoOutput){
255  // Accumulate a running total for the histogram file
256  auto ittot = fChannelMapProfTotal.find(block);
257  if(ittot == fChannelMapProfTotal.end())
258  fChannelMapProfTotal.insert(std::make_pair(block, *it->second));
259  else
260  ittot->second += *it->second;
261  }
262 
263  sr.put(std::move(it->second), TString::Format("block%d", block).Data());
264  }
265 
266  fChannelMapProf.clear();
267  }
268  else{
269  std::unique_ptr<caldp::AttenHistsMap> cm(new caldp::AttenHistsMap(fChannelMap));
270  sr.put(std::move(cm));
271 
272  // Accumulate a running total for the histogram file
274  // Don't double count for the output products
275  fChannelMap.Clear();
276  }
277 
278  return true;
279  }
280 
281  //......................................................................
283  {
284 
285  // handle the case where we want to use the raw ADC values rather than
286  // the PCHit PE value for finding the attenuation curves
287  if(fPECalculation == 1){
288  // get the RawDigit list
289  auto digitcol = evt.getValidHandle< std::vector<rawdata::RawDigit> >(fRawDataLabel);
290 
291  // make a map of the CellUniqueId to rawdigit
292  size_t plane = 0;
293  size_t cell = 0;
294  geo::CellUniqueId cellId = 0;
295  for(auto const& rawdigit : *digitcol){
296  plane = fCMap->GetPlane(&rawdigit);
297  cell = fCMap->GetCell(&rawdigit);
298  cellId = fGeom->Plane(plane)->Cell(cell)->Id();
299  fCellIdToRawDigit[cellId] = rawdigit;
300  }
301  }// end if using ADCs instead of PEs
302 
303  auto pchit_avg = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualAvg" ));
304  auto pchit_xy = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualXY" ));
305  auto pchit_xy_thresh = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualXYThresh"));
306  auto pchit_z = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualZ" ));
307  auto pchit_traj = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualTraj" ));
308 
309  this->FillPathLengthInfo(*pchit_z, caldp::kCalZ, evt.isRealData());
310  this->FillPathLengthInfo(*pchit_traj, caldp::kCalTraj, evt.isRealData());
311  this->FillPathLengthInfo(*pchit_xy, caldp::kCalXY, evt.isRealData());
312  this->FillPathLengthInfo(*pchit_avg, caldp::kCalAvg, evt.isRealData());
315  }
316 
317  // For automatic output file naming to work, need to allow at least one
318  // event per subrun for now. ART issue 5427
319  static unsigned int prevR = -1, prevSR = -1;
320  if(evt.run() != prevR || evt.subRun() != prevSR){
321  prevR = evt.run();
322  prevSR = evt.subRun();
323  return true;
324  }
325  return false;
326 
327 
328 // for (size_t idx = 0; idx < pchit_avg->size(); ++idx){
329 //
330 // caldp::PCHit current_hit = pchit_avg->at(idx);
331 // double w = current_hit.W();
332 // double path = current_hit.Path();
333 // double pe = this->PCHitPE(current_hit);
334 // int plane = current_hit.Plane();
335 // int cell = current_hit.Cell();
336 // const int view = this->GetView(plane);
337 //
338 // if(path > 10) continue;
339 //
340 // fCellPE[view]->Fill(cell, pe);
341 // fPlanePE[view]->Fill(plane, pe);
342 //
343 // // fSumChannels[view].WPE->Fill(w, pe);
344 // // fSumChannels[view].WPE_corr->Fill(w, pe/path);
345 //
346 // if(fProfileMode){
347 // // caldp::AttenProfiles& profs = GetChannelProfiles(evt.isRealData(), plane, cell);
348 // // profs.WPE.Fill(w, pe);
349 // // profs.WPE_corr.Fill(w, pe/path);
350 // }
351 // else{
352 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
353 // hists.WPE->Fill(w, pe);
354 // hists.WPE_corr->Fill(w, pe/path);
355 // }
356 //
357 // fPathLength[view]->Fill(path);
358 //
359 // } // loop over pc_hit_avg
360 //
361 // for(size_t idx = 0; idx < pchit_xy->size(); ++idx){
362 //
363 // caldp::PCHit current_hit = pchit_xy->at(idx);
364 // double w = current_hit.W();
365 // double path = current_hit.Path();
366 // double pe = this->PCHitPE(current_hit);
367 // double mev = current_hit.TrueMeV();
368 // int plane = current_hit.Plane();
369 // int cell = current_hit.Cell();
370 // const int view = GetView(plane);
371 // if(path > 10) continue;
372 //
373 // fSumChannels[view].WPE_corr_xy->Fill(w, pe/path);
374 // fSumChannels[view].WPE_corr_xy_truth->Fill(w, mev/path);
375 //
376 // if(fProfileMode){
377 // caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), evt.isRealData(), plane, cell);
378 // profs.WPE_corr_xy.Fill(w, pe/path);
379 // profs.WPE_corr_xy_truth.Fill(w, mev/path);
380 // }
381 // else{
382 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
383 // hists.WPE_corr_xy->Fill(w, pe/path);
384 // hists.WPE_corr_xy_truth->Fill(w, mev/path);
385 // }
386 // }// loop over pc_hit_xy
387 //
388 // if(fIncludeBelowThresholdZeros){
389 // for (size_t idx = 0; idx < pchit_xy_thresh->size(); ++idx){
390 //
391 // caldp::PCHit current_hit = pchit_xy_thresh->at(idx);
392 // double w = current_hit.W();
393 // double path = current_hit.Path();
394 // double pe = this->PCHitPE(current_hit);
395 // assert(pe == 0);
396 // int plane = current_hit.Plane();
397 // int cell = current_hit.Cell();
398 // const int view = GetView(plane);
399 //
400 // if(path > 10) continue;
401 //
402 // fSumChannels[view].WPE_corr_xy->Fill(w, 0);
403 //
404 // if(fProfileMode){
405 // caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), evt.isRealData(), plane, cell);
406 // profs.WPE_corr_xy.Fill(w, 0);
407 // }
408 // else{
409 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
410 // hists.WPE_corr_xy->Fill(w, 0);
411 // }
412 // // TODO: think about the other hists and profs
413 // }// loop over pc_hit_xy
414 // }
415 //
416 // for (size_t idx = 0; idx < pchit_z->size(); ++idx){
417 //
418 // caldp::PCHit current_hit = pchit_z->at(idx);
419 // double w = current_hit.W();
420 // double path = current_hit.Path();
421 // double pe = this->PCHitPE(current_hit);
422 // int plane = current_hit.Plane();
423 // int cell = current_hit.Cell();
424 // const int view = GetView(plane);
425 //
426 // if(path > 10) continue;
427 //
428 // fSumChannels[view].WPE_corr_z->Fill(w, pe/path);
429 //
430 // if(fProfileMode){
431 // caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), evt.isRealData(), plane, cell);
432 // profs.WPE_corr_z.Fill(w, pe/path);
433 // }
434 // else{
435 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
436 // hists.WPE_corr_z->Fill(w, pe/path);
437 // }
438 //
439 // }// loop over pc_hit_z
440 
441  }
442 
443  //......................................................................
444  void CosmicCalib::FillPathLengthInfo(std::vector<caldp::PCHit> const& pchits,
445  caldp::PathType_t const& pathType,
446  bool const& isRealData,
447  bool const& useBelowThreshold)
448  {
449  bool useProfiles = (fProfileMode && pathType != caldp::kCalAvg);
450  bool useHists = (!fProfileMode);
451 
452  for(auto const& current_hit : pchits){
453 
454  double w = current_hit.W();
455  double path = current_hit.Path();
456  double pe = this->PCHitPE(current_hit);
457  double mev = current_hit.TrueMeV();
458  int plane = current_hit.Plane();
459  int cell = current_hit.Cell();
460  const int view = this->GetView(plane);
461 
462  if(useBelowThreshold && pe != 0.)
463  throw cet::exception("CosmicCalib")
464  << "attempting to use below threshold hit with nonzero pe: " << pe;
465 
466  if(path > 10) continue;
467 
468  if(pathType == caldp::kCalAvg){
469  fCellPE[view]->Fill(cell, pe);
470  fPlanePE[view]->Fill(plane, pe);
471  fPathLength[view]->Fill(path);
472  }
473  else{
474  fSumChannels[view].HistogramByPathType(pathType)->Fill(w, pe/path);
475  if(pathType == caldp::kCalXY && !useBelowThreshold){
477  }
478  }
479 
480  if(useProfiles){
481  caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), isRealData, plane, cell);
482  profs.LiteProfileByPathType(pathType).Fill(w, pe/path);
483  if(pathType == caldp::kCalXY && !useBelowThreshold)
484  profs.LiteProfileByPathType(caldp::kCalXYTruth).Fill(w, mev/path);
485  }
486  else if(useHists){
487  caldp::AttenHists hists = GetChannelHists(isRealData, plane, cell);
488  hists.HistogramByPathType(pathType)->Fill(w, pe/path);
489  if(pathType == caldp::kCalAvg)
490  hists.HistogramByPathType(caldp::kCalUnknown)->Fill(w, pe/path);
491  else if(pathType == caldp::kCalXY && !useBelowThreshold)
492  hists.HistogramByPathType(caldp::kCalXYTruth)->Fill(w, mev/path);
493  }
494 
495  }
496 
497  return;
498  }
499 
500  //......................................................................
501  double CosmicCalib::PCHitPE(caldp::PCHit const& pchit)
502  {
503  // this method provides multiple ways of getting a "PE" value
504  // and putting that into the plots used to make the attenuation
505  // curves. What value is returned is based on a parameter set at
506  // run time
507  double pe = 0.;
508 
509  if(fPECalculation == 0) return pchit.PE();
510  else if(fPECalculation == 1){
511  // get the cellId for the PCHit
512  auto pchitId = fGeom->Plane(pchit.Plane())->Cell(pchit.Cell())->Id();
513 
514  if( fCellIdToRawDigit.count(pchitId) > 0) pe = 1.*fCellIdToRawDigit.find(pchitId)->second.ADC();
515  else
516  throw cet::exception("CosmicCalib")
517  << "Cannot find RawDigit for " << pchit.Plane() << " " << pchit.Cell();
518 
519  }// end if using raw digit ADC values
520 
521  return pe;
522  }
523 
524  //......................................................................
526  int const& plane,
527  int const& cell)
528  {
529  // plane is a const reference so we can't change it.
530  // instead make a new variable called pln that we can change.
531  // if we are consolidating the views, set pln to the view
532  int pln = plane;
534  pln = GetView(plane);
535  }
536 
538 
539  // If the histograms are created for this channel, just return them
540  if(hists.WPE) return hists;
541 
542  // Otherwise have to construct them
543  const double maxW = fGeom->Plane(plane)->Cell(0)->HalfL()*1.15;
544 
545  TH2F* wadc = new TH2F(TString::Format("h_WPE_%03d_%03d", pln, cell),
546  ";W;PE", 100, -maxW, +maxW, 100, 0, 1000);
547 
548  TH2F* wadc_corr = new TH2F(TString::Format("h_WPE_corr_avg_%03d_%03d",
549  pln, cell),
550  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
551 
552  TH2F* wadc_corr_xy = new TH2F(TString::Format("h_WPE_corr_xy_%03d_%03d",
553  pln, cell),
554  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
555 
556  TH2F* wadc_corr_xy_truth = new TH2F(TString::Format("h_WPE_corr_xy_truth_%03d_%03d", pln, cell),
557  ";W;MeV/cm", 100, -maxW, +maxW, 100, 0, 10);
558 
559 
560  TH2F* wadc_corr_z = new TH2F(TString::Format("h_WPE_corr_z_%03d_%03d", pln, cell),
561  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
562 
563  TH2F* wadc_corr_traj = new TH2F(TString::Format("h_WPE_corr_traj_%03d_%03d", pln, cell),
564  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
565 
566  hists.WPE = wadc;
567  hists.WPE_corr = wadc_corr;
568  hists.WPE_corr_xy = wadc_corr_xy;
569  hists.WPE_corr_xy_truth = wadc_corr_xy_truth;
570  hists.WPE_corr_z = wadc_corr_z;
571  hists.WPE_corr_traj = wadc_corr_traj;
572 
573  return hists;
574  }
575 
576  //......................................................................
579  bool const& isData,
580  int const& plane,
581  int const& cell)
582  {
583  // MC only varies for these reasons. Consolidate due to low stats.
584  // Make a new variable, pln to be either the plane or the view if
585  // consolidating
586  int pln = plane;
588  pln = GetView(plane);
589  }
590 
591  double maxW;
592  if(det == novadaq::cnv::kFARDET) maxW = 800;
593  else if(det == novadaq::cnv::kNEARDET) maxW = 250;
594  else maxW = 250; // Don't bother to make non-square NDOS profiles
595 
596  const int block = pln/32;
597  if(fChannelMapProf.find(block) == fChannelMapProf.end())
598  fChannelMapProf.emplace(std::make_pair(block, std::unique_ptr<caldp::AttenProfilesMap>(new caldp::AttenProfilesMap(-maxW, +maxW))));
599  return fChannelMapProf.find(block)->second->GetProfiles(pln, cell);
600  }
601 
602  //......................................................................
603  int CosmicCalib::GetView(int const& plane){
604  //get the view for the current combination of view and plane (x,y,mux, or muy)
605  int view=0;
607  int firstPlane = fGeom->FirstPlaneInMuonCatcher();
608  int isMuonCatcher = firstPlane <= plane;
609  view = fGeom->Plane(plane)->View()+(2*isMuonCatcher);
610  }
611  else{
612  view = fGeom->Plane(plane)->View();
613  }
614  return view;
615  }
616 
617  //......................................................................
619  {
620  if(!fHistoOutput) return;
621 
623 
624  for(auto& it: fChannelMapTotal.GetAllHistsByPlane()){
625  const int plane = it.first;
626  const std::vector<TH2F*>& hs = it.second;
627 
628  art::TFileDirectory dir = tfs->mkdir(Form("plane_%03d", plane));
629 
630  for(TH2F* h: hs) dir.make<TH2F>(*h);
631  }
632 
633  for(int block = 0; block < 28; ++block){
634  auto itmap = fChannelMapProfTotal.find(block);
635  if(itmap == fChannelMapProfTotal.end()) continue;
636  for(auto& it: itmap->second.GetAllProfilesByPlaneAndCell()){
637  const int plane = it.first;
638  std::cout << plane << std::endl;
639  if(it.second.empty()) continue;
640 
641  art::TFileDirectory dirBlock = tfs->mkdir(Form("block_%02d", plane/32));
642  art::TFileDirectory dirPlane = dirBlock.mkdir(Form("plane_%03d", plane));
643 
644  for(auto& it2: it.second){
645  const int cell = it2.first;
646  const std::vector<TH1F*>& hs = it2.second;
647  if(hs.empty()) continue;
648 
649  art::TFileDirectory dirMod = dirPlane.mkdir(Form("module_%02d", cell/32));
650  art::TFileDirectory dirCell = dirMod.mkdir(Form("cell_%03d", cell));
651 
652  for(TH1F* h: hs) dirCell.make<TH1F>(*h);
653  }
654  }
655  }
656  }
657 
659 
660 } // end namespace
661 ////////////////////////////////////////////////////////////////////////
std::map< int, std::vector< TH2F * > > GetAllHistsByPlane()
Intended for writing out/drawing all the histograms.
Definition: AttenHists.h:92
double HalfL() const
Definition: CellGeo.cxx:198
caldp::AttenHists GetChannelHists(bool const &isData, int const &plane, int const &cell)
virtual bool filter(art::Event &evt)
set< int >::iterator it
caldp::AttenHistsMap fChannelMapTotal
This is what we put in the histogram file.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
Vertical planes which measure X.
Definition: PlaneGeo.h:28
void Fill(float x, float y)
int Cell() const
Return cell value.
Definition: PCHit.h:26
AttenProfiles for many channels.
Definition: AttenProfiles.h:89
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
art::ServiceHandle< cmap::CMap > fCMap
connection map
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const PlaneGeo * Plane(unsigned int i) const
void reconfigure(const fhicl::ParameterSet &pset)
TString hists[nhists]
Definition: bdt_com.C:3
DEFINE_ART_MODULE(TestTMapFile)
caldp::AttenHistsMap fChannelMap
This is where GetChannelHist stores its histograms. Map is from block #.
TH2F * WPE_corr_xy_truth
Definition: AttenHists.h:67
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
"Pre-calibration hit". Common input to calibration procedures
Definition: PCHit.h:16
Definition: Run.h:21
Give every cell in the geometry a unique ID number based on the TGeo path to the node.
const CellUniqueId & Id() const
Definition: CellGeo.h:55
art::ServiceHandle< geo::Geometry > fGeom
bool isRealData() const
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
std::map< geo::CellUniqueId, rawdata::RawDigit > fCellIdToRawDigit
map cell id to the raw digit
Far Detector at Ash River, MN.
CDPStorage service.
block
print "ROW IS " print row
Definition: elec2geo.py:31
float PE() const
Return PE value.
Definition: PCHit.h:38
virtual bool endSubRun(art::SubRun &sr)
enum caldp::_path_type PathType_t
std::map< int, caldp::AttenProfilesMap > fChannelMapProfTotal
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
int GetView(int const &plane)
AttenHists & GetHists(const geo::OfflineChan &c)
Definition: AttenHists.h:81
Near Detector in the NuMI cavern.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:491
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
double PCHitPE(caldp::PCHit const &pchit)
Histograms used by attenuation calibration.
Definition: AttenHists.h:27
caf::StandardRecord * sr
SubRunNumber_t subRun() const
double DetHalfHeight() const
RunNumber_t run() const
TH2F * HistogramByPathType(caldp::PathType_t const &pt)
Definition: AttenHists.h:52
int Plane() const
Return plane value.
Definition: PCHit.h:24
AttenHists for many channels.
Definition: AttenHists.h:77
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
caldp::AttenHists fSumChannels[4]
Profiles used by attenuation calibration.
Definition: AttenProfiles.h:61
const std::string path
Definition: plot_BEN.C:43
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
LiteProfile & LiteProfileByPathType(caldp::PathType_t const &pt)
Definition: AttenProfiles.h:71
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.h:22
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
void FillPathLengthInfo(std::vector< caldp::PCHit > const &pchits, caldp::PathType_t const &pathType, bool const &isRealData, bool const &useBelowThreshold=false)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
CosmicCalib(const fhicl::ParameterSet &pset)
virtual bool beginRun(art::Run &run)
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
Float_t w
Definition: plot.C:20
TH2F * WPE_corr_traj
Definition: AttenHists.h:67
static constexpr Double_t cm
Definition: Munits.h:140
caldp::AttenProfiles & GetChannelProfiles(novadaq::cnv::DetId const &det, bool const &isData, int const &plane, int const &cell)
Encapsulate the geometry of one entire detector (near, far, ndos)
std::map< int, std::unique_ptr< caldp::AttenProfilesMap > > fChannelMapProf
TH2F * WPE_corr_xy
Definition: AttenHists.h:67
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string