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
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  : fdx(0) // This one is the signal that no histograms are created yet
112  {
113  reconfigure(pset);
114 
115  if(fProfileMode){
116  // TODO: I don't think we can determine the detector in time...
117  for(int block = 0; block < 28; ++block){
118  produces<caldp::AttenProfilesMap, art::InSubRun>(TString::Format("block%d", block).Data());
119  }
120  }
121  else{
122  produces<caldp::AttenHistsMap, art::InRun>();
123  }
124  }
125 
126  //......................................................................
128  {
129  }
130 
131  //......................................................................
132 
134  {
135  fPCHitsLabel = pset.get<std::string>("PCHitsLabel");
136  fIncludeBelowThresholdZeros = pset.get<bool>("IncludeBelowThresholdZeros");
137  fProfileMode = pset.get<bool>("ProfileMode");
138  fConsolidateViewsMode = pset.get<bool>("ConsolidateViewsMode");
139  fHistoOutput = pset.get<bool>("HistoOutput");
140  fPECalculation = pset.get<int >("PECalculation", 0);
141  fRawDataLabel = pset.get<std::string>("RawDataLabel", "daq");
142  }
143 
144  //......................................................................
146  {
147  }
148 
149  //......................................................................
151  {
154 
155  // The rest is making histograms. Don't do it again on a subsequent run
156  if(fdx) return true;
157 
158  const TString viewStrArr[4]={"X","Y","muX","muY"};
160  bool isND = fGeom->DetId() == novadaq::cnv::kNEARDET;
161  for(int view = geo::kX; view <= (geo::kY+2*isND); ++view){
162  const TString viewStr = viewStrArr[view];
163 
164  double maxW;
166  maxW = 900;
167  }
168  else if(fGeom->DetId() == novadaq::cnv::kNEARDET){
169  maxW = 250;
170  }
171  else{
172  // NDOS (hopefully)
173  maxW = (view == geo::kX) ? 250 : 150;
174  }
175 
177 
178  v.WPE = tfs->make<TH2F>("h_WPE_in"+viewStr,
179  ";W;PE",
180  2*maxW, -maxW, +maxW,
181  1000, 0, 1000);
182 
183  v.WPE_corr = tfs->make<TH2F>("h_WPE_corr_avg_in"+viewStr,
184  ";W;PE/cm",
185  2*maxW, -maxW, +maxW,
186  1000, 0, 200);
187 
188  v.WPE_corr_xy = tfs->make<TH2F>("h_WPE_corr_xy_in"+viewStr,
189  ";W;PE/cm",
190  2*maxW, -maxW, +maxW,
191  1000, 0, 200);
192 
193  v.WPE_corr_xy_truth = tfs->make<TH2F>("h_WPE_corr_xy_truth_in"+viewStr,
194  ";W;MeV/cm",
195  2*maxW, -maxW, +maxW,
196  1000, 0, 10);
197 
198  v.WPE_corr_z = tfs->make<TH2F>("h_WPE_corr_z_in"+viewStr,
199  ";W;PE/cm",
200  2*maxW, -maxW, +maxW,
201  1000, 0, 200);
202  v.WPE_corr_traj = tfs->make<TH2F>("h_WPE_corr_traj_in"+viewStr,
203  ";W;PE/cm",
204  2*maxW, -maxW, +maxW,
205  1000, 0, 200);
206 
207  fSumChannels[view] = v;
208 
209  fCellPE[view] = tfs->make<TH2F>("h_CellPE_in"+viewStr, ";Cell;PE",
210  105, -5, 100, 1000, 0, 1000);
211 
212  fPlanePE[view] = tfs->make<TH2F>("h_PlanePE_in"+viewStr, ";Plane;PE",
213  103, -6, 200, 1000, 0, 1000);
214 
215  fPathLength[view] = tfs->make<TH1F>("h_PathLength"+viewStr,
216  ";Path length (cm)",
217  80, 0, 20);
218  } // end for view
219 
220  fdx = tfs->make<TH1F>("h_dx", ";|dx|", 100, 0, 1);
221  fdy = tfs->make<TH1F>("h_dy", ";|dy|", 100, 0, 1);
222  fdz = tfs->make<TH1F>("h_dz", ";|dz|", 100, 0, 1);
223 
224  fdxdy = tfs->make<TH2F>("h_dxdy", ";|dx|;|dy|", 100, 0, 1, 100, 0, 1);
225  fdxdz = tfs->make<TH2F>("h_dxdz", ";|dx|;|dz|", 100, 0, 1, 100, 0, 1);
226  fdydz = tfs->make<TH2F>("h_dydz", ";|dy|;|dz|", 100, 0, 1, 100, 0, 1);
227 
228  fStartXY = tfs->make<TH2F>("h_StartXY", ";x;y",
229  100, -200, +200, 100, -250, +250);
230  fStartZX = tfs->make<TH2F>("h_StartZX", ";z;x",
231  100, -50, +900, 100, -200, +200);
232  fStartZY = tfs->make<TH2F>("h_StartZY", ";z;y",
233  100, -50, +900, 100, -250, +250);
234 
235  fStopXY = tfs->make<TH2F>("h_StopXY", ";x;y",
236  100, -200, +200, 100, -250, +250);
237  fStopZX = tfs->make<TH2F>("h_StopZX", ";z;x",
238  100, -50, +900, 100, -200, +200);
239  fStopZY = tfs->make<TH2F>("h_StopZY", ";z;y",
240  100, -50, +900, 100, -250, +250);
241 
242  return true;
243  }
244 
245  //......................................................................
247  {
248  if(fProfileMode){
249  for(int block = 0; block < 28; ++block){
250  auto it = fChannelMapProf.find(block);
251  if(it == fChannelMapProf.end()) continue;
252 
253  if(fHistoOutput){
254  // Accumulate a running total for the histogram file
255  auto ittot = fChannelMapProfTotal.find(block);
256  if(ittot == fChannelMapProfTotal.end())
257  fChannelMapProfTotal.insert(std::make_pair(block, *it->second));
258  else
259  ittot->second += *it->second;
260  }
261 
262  sr.put(std::move(it->second), TString::Format("block%d", block).Data());
263  }
264 
265  fChannelMapProf.clear();
266  }
267  else{
268  std::unique_ptr<caldp::AttenHistsMap> cm(new caldp::AttenHistsMap(fChannelMap));
269  sr.put(std::move(cm));
270 
271  // Accumulate a running total for the histogram file
273  // Don't double count for the output products
274  fChannelMap.Clear();
275  }
276 
277  return true;
278  }
279 
280  //......................................................................
282  {
283 
284  // handle the case where we want to use the raw ADC values rather than
285  // the PCHit PE value for finding the attenuation curves
286  if(fPECalculation == 1){
287  // get the RawDigit list
288  auto digitcol = evt.getValidHandle< std::vector<rawdata::RawDigit> >(fRawDataLabel);
289 
290  // make a map of the CellUniqueId to rawdigit
291  size_t plane = 0;
292  size_t cell = 0;
293  geo::CellUniqueId cellId = 0;
294  for(auto const& rawdigit : *digitcol){
295  plane = fCMap->GetPlane(&rawdigit);
296  cell = fCMap->GetCell(&rawdigit);
297  cellId = fGeom->Plane(plane)->Cell(cell)->Id();
298  fCellIdToRawDigit[cellId] = rawdigit;
299  }
300  }// end if using ADCs instead of PEs
301 
302  auto pchit_avg = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualAvg" ));
303  auto pchit_xy = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualXY" ));
304  auto pchit_xy_thresh = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualXYThresh"));
305  auto pchit_z = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualZ" ));
306  auto pchit_traj = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualTraj" ));
307 
308  this->FillPathLengthInfo(*pchit_z, caldp::kCalZ, evt.isRealData());
309  this->FillPathLengthInfo(*pchit_traj, caldp::kCalTraj, evt.isRealData());
310  this->FillPathLengthInfo(*pchit_xy, caldp::kCalXY, evt.isRealData());
311  this->FillPathLengthInfo(*pchit_avg, caldp::kCalAvg, evt.isRealData());
314  }
315 
316  // For automatic output file naming to work, need to allow at least one
317  // event per subrun for now. ART issue 5427
318  static unsigned int prevR = -1, prevSR = -1;
319  if(evt.run() != prevR || evt.subRun() != prevSR){
320  prevR = evt.run();
321  prevSR = evt.subRun();
322  return true;
323  }
324  return false;
325 
326 
327 // for (size_t idx = 0; idx < pchit_avg->size(); ++idx){
328 //
329 // caldp::PCHit current_hit = pchit_avg->at(idx);
330 // double w = current_hit.W();
331 // double path = current_hit.Path();
332 // double pe = this->PCHitPE(current_hit);
333 // int plane = current_hit.Plane();
334 // int cell = current_hit.Cell();
335 // const int view = this->GetView(plane);
336 //
337 // if(path > 10) continue;
338 //
339 // fCellPE[view]->Fill(cell, pe);
340 // fPlanePE[view]->Fill(plane, pe);
341 //
342 // // fSumChannels[view].WPE->Fill(w, pe);
343 // // fSumChannels[view].WPE_corr->Fill(w, pe/path);
344 //
345 // if(fProfileMode){
346 // // caldp::AttenProfiles& profs = GetChannelProfiles(evt.isRealData(), plane, cell);
347 // // profs.WPE.Fill(w, pe);
348 // // profs.WPE_corr.Fill(w, pe/path);
349 // }
350 // else{
351 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
352 // hists.WPE->Fill(w, pe);
353 // hists.WPE_corr->Fill(w, pe/path);
354 // }
355 //
356 // fPathLength[view]->Fill(path);
357 //
358 // } // loop over pc_hit_avg
359 //
360 // for(size_t idx = 0; idx < pchit_xy->size(); ++idx){
361 //
362 // caldp::PCHit current_hit = pchit_xy->at(idx);
363 // double w = current_hit.W();
364 // double path = current_hit.Path();
365 // double pe = this->PCHitPE(current_hit);
366 // double mev = current_hit.TrueMeV();
367 // int plane = current_hit.Plane();
368 // int cell = current_hit.Cell();
369 // const int view = GetView(plane);
370 // if(path > 10) continue;
371 //
372 // fSumChannels[view].WPE_corr_xy->Fill(w, pe/path);
373 // fSumChannels[view].WPE_corr_xy_truth->Fill(w, mev/path);
374 //
375 // if(fProfileMode){
376 // caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), evt.isRealData(), plane, cell);
377 // profs.WPE_corr_xy.Fill(w, pe/path);
378 // profs.WPE_corr_xy_truth.Fill(w, mev/path);
379 // }
380 // else{
381 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
382 // hists.WPE_corr_xy->Fill(w, pe/path);
383 // hists.WPE_corr_xy_truth->Fill(w, mev/path);
384 // }
385 // }// loop over pc_hit_xy
386 //
387 // if(fIncludeBelowThresholdZeros){
388 // for (size_t idx = 0; idx < pchit_xy_thresh->size(); ++idx){
389 //
390 // caldp::PCHit current_hit = pchit_xy_thresh->at(idx);
391 // double w = current_hit.W();
392 // double path = current_hit.Path();
393 // double pe = this->PCHitPE(current_hit);
394 // assert(pe == 0);
395 // int plane = current_hit.Plane();
396 // int cell = current_hit.Cell();
397 // const int view = GetView(plane);
398 //
399 // if(path > 10) continue;
400 //
401 // fSumChannels[view].WPE_corr_xy->Fill(w, 0);
402 //
403 // if(fProfileMode){
404 // caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), evt.isRealData(), plane, cell);
405 // profs.WPE_corr_xy.Fill(w, 0);
406 // }
407 // else{
408 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
409 // hists.WPE_corr_xy->Fill(w, 0);
410 // }
411 // // TODO: think about the other hists and profs
412 // }// loop over pc_hit_xy
413 // }
414 //
415 // for (size_t idx = 0; idx < pchit_z->size(); ++idx){
416 //
417 // caldp::PCHit current_hit = pchit_z->at(idx);
418 // double w = current_hit.W();
419 // double path = current_hit.Path();
420 // double pe = this->PCHitPE(current_hit);
421 // int plane = current_hit.Plane();
422 // int cell = current_hit.Cell();
423 // const int view = GetView(plane);
424 //
425 // if(path > 10) continue;
426 //
427 // fSumChannels[view].WPE_corr_z->Fill(w, pe/path);
428 //
429 // if(fProfileMode){
430 // caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), evt.isRealData(), plane, cell);
431 // profs.WPE_corr_z.Fill(w, pe/path);
432 // }
433 // else{
434 // caldp::AttenHists hists = GetChannelHists(evt.isRealData(), plane, cell);
435 // hists.WPE_corr_z->Fill(w, pe/path);
436 // }
437 //
438 // }// loop over pc_hit_z
439 
440  }
441 
442  //......................................................................
443  void CosmicCalib::FillPathLengthInfo(std::vector<caldp::PCHit> const& pchits,
444  caldp::PathType_t const& pathType,
445  bool const& isRealData,
446  bool const& useBelowThreshold)
447  {
448  bool useProfiles = (fProfileMode && pathType != caldp::kCalAvg);
449  bool useHists = (!fProfileMode);
450 
451  for(auto const& current_hit : pchits){
452 
453  double w = current_hit.W();
454  double path = current_hit.Path();
455  double pe = this->PCHitPE(current_hit);
456  double mev = current_hit.TrueMeV();
457  int plane = current_hit.Plane();
458  int cell = current_hit.Cell();
459  const int view = this->GetView(plane);
460 
461  if(useBelowThreshold && pe != 0.)
462  throw cet::exception("CosmicCalib")
463  << "attempting to use below threshold hit with nonzero pe: " << pe;
464 
465  if(path > 10) continue;
466 
467  if(pathType == caldp::kCalAvg){
468  fCellPE[view]->Fill(cell, pe);
469  fPlanePE[view]->Fill(plane, pe);
470  fPathLength[view]->Fill(path);
471  }
472  else{
473  fSumChannels[view].HistogramByPathType(pathType)->Fill(w, pe/path);
474  if(pathType == caldp::kCalXY && !useBelowThreshold){
476  }
477  }
478 
479  if(useProfiles){
480  caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), isRealData, plane, cell);
481  profs.LiteProfileByPathType(pathType).Fill(w, pe/path);
482  if(pathType == caldp::kCalXY && !useBelowThreshold)
483  profs.LiteProfileByPathType(caldp::kCalXYTruth).Fill(w, mev/path);
484  }
485  else if(useHists){
486  caldp::AttenHists hists = GetChannelHists(isRealData, plane, cell);
487  hists.HistogramByPathType(pathType)->Fill(w, pe/path);
488  if(pathType == caldp::kCalAvg)
489  hists.HistogramByPathType(caldp::kCalUnknown)->Fill(w, pe/path);
490  else if(pathType == caldp::kCalXY && !useBelowThreshold)
491  hists.HistogramByPathType(caldp::kCalXYTruth)->Fill(w, mev/path);
492  }
493 
494  }
495 
496  return;
497  }
498 
499  //......................................................................
500  double CosmicCalib::PCHitPE(caldp::PCHit const& pchit)
501  {
502  // this method provides multiple ways of getting a "PE" value
503  // and putting that into the plots used to make the attenuation
504  // curves. What value is returned is based on a parameter set at
505  // run time
506  double pe = 0.;
507 
508  if(fPECalculation == 0) return pchit.PE();
509  else if(fPECalculation == 1){
510  // get the cellId for the PCHit
511  auto pchitId = fGeom->Plane(pchit.Plane())->Cell(pchit.Cell())->Id();
512 
513  if( fCellIdToRawDigit.count(pchitId) > 0) pe = 1.*fCellIdToRawDigit.find(pchitId)->second.ADC();
514  else
515  throw cet::exception("CosmicCalib")
516  << "Cannot find RawDigit for " << pchit.Plane() << " " << pchit.Cell();
517 
518  }// end if using raw digit ADC values
519 
520  return pe;
521  }
522 
523  //......................................................................
525  int const& plane,
526  int const& cell)
527  {
528  // plane is a const reference so we can't change it.
529  // instead make a new variable called pln that we can change.
530  // if we are consolidating the views, set pln to the view
531  int pln = plane;
533  pln = GetView(plane);
534  }
535 
537 
538  // If the histograms are created for this channel, just return them
539  if(hists.WPE) return hists;
540 
541  // Otherwise have to construct them
542  const double maxW = fGeom->Plane(plane)->Cell(0)->HalfL()*1.15;
543 
544  TH2F* wadc = new TH2F(TString::Format("h_WPE_%03d_%03d", pln, cell),
545  ";W;PE", 100, -maxW, +maxW, 100, 0, 1000);
546 
547  TH2F* wadc_corr = new TH2F(TString::Format("h_WPE_corr_avg_%03d_%03d",
548  pln, cell),
549  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
550 
551  TH2F* wadc_corr_xy = new TH2F(TString::Format("h_WPE_corr_xy_%03d_%03d",
552  pln, cell),
553  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
554 
555  TH2F* wadc_corr_xy_truth = new TH2F(TString::Format("h_WPE_corr_xy_truth_%03d_%03d", pln, cell),
556  ";W;MeV/cm", 100, -maxW, +maxW, 100, 0, 10);
557 
558 
559  TH2F* wadc_corr_z = new TH2F(TString::Format("h_WPE_corr_z_%03d_%03d", pln, cell),
560  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
561 
562  TH2F* wadc_corr_traj = new TH2F(TString::Format("h_WPE_corr_traj_%03d_%03d", pln, cell),
563  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
564 
565  hists.WPE = wadc;
566  hists.WPE_corr = wadc_corr;
567  hists.WPE_corr_xy = wadc_corr_xy;
568  hists.WPE_corr_xy_truth = wadc_corr_xy_truth;
569  hists.WPE_corr_z = wadc_corr_z;
570  hists.WPE_corr_traj = wadc_corr_traj;
571 
572  return hists;
573  }
574 
575  //......................................................................
578  bool const& isData,
579  int const& plane,
580  int const& cell)
581  {
582  // MC only varies for these reasons. Consolidate due to low stats.
583  // Make a new variable, pln to be either the plane or the view if
584  // consolidating
585  int pln = plane;
587  pln = GetView(plane);
588  }
589 
590  double maxW;
591  if(det == novadaq::cnv::kFARDET) maxW = 800;
592  else if(det == novadaq::cnv::kNEARDET) maxW = 250;
593  else maxW = 250; // Don't bother to make non-square NDOS profiles
594 
595  const int block = pln/32;
596  if(fChannelMapProf.find(block) == fChannelMapProf.end())
597  fChannelMapProf.emplace(std::make_pair(block, std::unique_ptr<caldp::AttenProfilesMap>(new caldp::AttenProfilesMap(-maxW, +maxW))));
598  return fChannelMapProf.find(block)->second->GetProfiles(pln, cell);
599  }
600 
601  //......................................................................
602  int CosmicCalib::GetView(int const& plane){
603  //get the view for the current combination of view and plane (x,y,mux, or muy)
604  int view=0;
606  int firstPlane = fGeom->FirstPlaneInMuonCatcher();
607  int isMuonCatcher = firstPlane <= plane;
608  view = fGeom->Plane(plane)->View()+(2*isMuonCatcher);
609  }
610  else{
611  view = fGeom->Plane(plane)->View();
612  }
613  return view;
614  }
615 
616  //......................................................................
618  {
619  if(!fHistoOutput) return;
620 
622 
623  for(auto& it: fChannelMapTotal.GetAllHistsByPlane()){
624  const int plane = it.first;
625  const std::vector<TH2F*>& hs = it.second;
626 
627  art::TFileDirectory dir = tfs->mkdir(Form("plane_%03d", plane));
628 
629  for(TH2F* h: hs) dir.make<TH2F>(*h);
630  }
631 
632  for(int block = 0; block < 28; ++block){
633  auto itmap = fChannelMapProfTotal.find(block);
634  if(itmap == fChannelMapProfTotal.end()) continue;
635  for(auto& it: itmap->second.GetAllProfilesByPlaneAndCell()){
636  const int plane = it.first;
637  std::cout << plane << std::endl;
638  if(it.second.empty()) continue;
639 
640  art::TFileDirectory dirBlock = tfs->mkdir(Form("block_%02d", plane/32));
641  art::TFileDirectory dirPlane = dirBlock.mkdir(Form("plane_%03d", plane));
642 
643  for(auto& it2: it.second){
644  const int cell = it2.first;
645  const std::vector<TH1F*>& hs = it2.second;
646  if(hs.empty()) continue;
647 
648  art::TFileDirectory dirMod = dirPlane.mkdir(Form("module_%02d", cell/32));
649  art::TFileDirectory dirCell = dirMod.mkdir(Form("cell_%03d", cell));
650 
651  for(TH1F* h: hs) dirCell.make<TH1F>(*h);
652  }
653  }
654  }
655  }
656 
658 
659 } // end namespace
660 ////////////////////////////////////////////////////////////////////////
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)
SubRunNumber_t subRun() const
Definition: Event.h:72
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
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
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)
bool isRealData() const
Definition: Event.h:83
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:31
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
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.
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
double DetHalfHeight() 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
T * make(ARGS...args) const
LiteProfile & LiteProfileByPathType(caldp::PathType_t const &pt)
Definition: AttenProfiles.h:71
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
ProductID put(std::unique_ptr< PROD > &&)
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)
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
virtual bool beginRun(art::Run &run)
RunNumber_t run() const
Definition: Event.h:77
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