MakeAttenuationProfiles_plugin.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file MakeAttenuationProfiles_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 #include <list>
7 #include <map>
8 #include <string>
9 
10 // ROOT includes
11 #include "TH2.h"
12 #include "TProfile.h"
13 #include "TCanvas.h"
14 
15 // Framework includes
23 #include "cetlib_except/exception.h"
24 
25 // NOvASoft includes
29 #include "Geometry/Geometry.h"
34 #include "NovaDAQConventions/DAQConventions.h"
35 #include "RawData/RawDigit.h"
37 
38 namespace calib
39 {
41  {
42  public:
43  explicit MakeAttenuationProfiles(const fhicl::ParameterSet& pset);
45 
46  virtual void writeResults(art::Results & r);
47  virtual void event(art::Event const& evt);
48 
49  virtual void reconfigure(const fhicl::ParameterSet& pset);
50 
51  virtual void beginRun(art::Run const& run);
52  virtual void endJob();
53  virtual void clear();
54 
55  protected:
56  caldp::AttenHists GetChannelHists(bool const& isData,
57  int const& plane,
58  int const& cell);
60  bool const& isData,
61  int const& plane,
62  int const& cell);
63  int GetView(int const& plane);
64  //int GetFiberBrightness(int const& plane,
65  // int const& cell);
66  double PCHitPE(caldp::PCHit const& pchit);
67 
68  void FillPathLengthInfo(std::vector<caldp::PCHit> const& pchits,
69  caldp::PathType_t const& pathType,
70  bool const& isRealData,
71  bool const& useBelowThreshold=false);
72 
74 
76 
77  double fInitEdgeX;
78  double fInitEdgeY;
79  double fTrajectory_pathLow; ///< Condition on trajectory path to be
80  ///< greater than fTrajectory_pathLow (cm).
81  double fTrajectory_pathHigh; ///< Condition on trajectory path to be
82  ///< less than fTrajectory_pathHigh (cm).
85 
92  unsigned int fFiberBrightnessBins;
94  int fPECalculation; ///< method of determining the value of a pe
95  ///< 0 = PCHit::PE(), 1 = DCS ADC value from raw digit
96 
97 
98  // General histograms for investigating variables.
99  // Some of these could be derived from fChannelMap, but it's easier
100  // just to fill them as we go.
101 
103 
104  std::list<std::pair<int, std::string> > fInstanceLabels; ///< labels for the different products
105  std::map<geo::CellUniqueId, rawdata::RawDigit> fCellIdToRawDigit; ///< map cell id to the raw digit
108 
109  TH2* fCellPE[12][4];
110  TH2* fCellRawPE[12][4];
111  TH2* fPlanePE[12][4];
112  TH2* fBrightnessPE[4];
113  TProfile* fBrightnessProf[4];
114  TCanvas* fBrightnessCanv[4];
115 
116  // Direction cosines
117  TH1 *fdx;
118  TH1 *fdy;
119  TH1 *fdz;
120  TH2 *fdxdy;
121  TH2 *fdxdz;
122  TH2 *fdydz;
123 
124  // Track vertices
125  TH2 *fStartXY;
126  TH2 *fStartZX;
127  TH2 *fStartZY;
128  TH2 *fStopXY;
129  TH2 *fStopZX;
130  TH2 *fStopZY;
131 
132  TH1* fPathLength[12][4];
133 
134  //Check on fibre brightness
136 
137  //Check that there are entries for all planes and cells
139 
140  /// This is where \ref GetChannelHist stores its histograms. Map is from block #
142  std::map<int, std::unique_ptr<caldp::AttenProfilesMap>> fChannelMapProf;
143 
144  /// This is what we put in the histogram file
146  std::map<int, caldp::AttenProfilesMap> fChannelMapProfTotal;
147 
148  };
149 
150  //......................................................................
152  : fdx(0)
153  {
154  reconfigure(pset);
155 
156  //make a profile for each block OR for each fibre brightness
158  for(int fb = 0; fb < (int)fFiberBrightnessBins; fb++){
159  std::ostringstream fbLabel;
160  fbLabel << "fiberbrightness" << fb;
161  fInstanceLabels.push_back(std::make_pair(fb,
162  fbLabel.str())
163  );
164  }
165  }
167  std::ostringstream mucLabel;
168  mucLabel << "nofbmuC";
169  fInstanceLabels.push_back(std::make_pair(100,
170  mucLabel.str())
171  );
172  }
174  for(int block = 0; block < 28; ++block){
175  std::ostringstream blockLabel;
176  blockLabel << "block" << block;
177  fInstanceLabels.push_back(std::make_pair(block,
178  blockLabel.str())
179  );
180  }
181  }
182 
183  if(fProfileMode){
184  // TODO: I don't think we can determine the detector in time...
185  for(auto label : fInstanceLabels){
186  produces<caldp::AttenProfilesMap>(label.second);
187  }
188  }
189  else{
190  produces<caldp::AttenHistsMap>();
191  }
192  }
193 
194  //......................................................................
196  {
197  }
198 
199  //......................................................................
201  {
202  fChannelMapProf.clear();
203  fChannelMapProfTotal.clear();
204  }
205 
206  //......................................................................
208  {
209  fPCHitsLabel = pset.get<std::string >("PCHitsLabel" );
210  fIncludeBelowThresholdZeros = pset.get<bool >("IncludeBelowThresholdZeros");
211  fProfileMode = pset.get<bool >("ProfileMode" );
212  fConsolidateViewsMode = pset.get<bool >("ConsolidateViewsMode" );
213  fMuCConsolidateViewsMode = pset.get<bool >("MuCConsolidateViewsMode" );
214  fFiberBrightnessMode = pset.get<bool >("FiberBrightnessMode" );
215  fMuCFiberBrightnessMode = pset.get<bool >("MuCFiberBrightnessMode" );
216  fCutBeamSpills = pset.get<bool >("CutBeamSpills" );
217  fFiberBrightnessBins = pset.get<unsigned int>("FiberBrightnessBins" );
218  fHistoOutput = pset.get<bool >("HistoOutput" );
219  fPECalculation = pset.get<int >("PECalculation", 0 );
220  fRawDataLabel = pset.get<std::string >("RawDataLabel", "daq");
221  fTrajectory_pathLow = pset.get<double >("Trajectory_pathLow" );
222  fTrajectory_pathHigh = pset.get<double >("Trajectory_pathHigh" );
223 
224  //fFiberBrightnessMode = true;
225 
226  if(fFiberBrightnessMode) LOG_VERBATIM("MakeAttenuationProfiles") << "fFiberBrightnessMode is on";
227  else LOG_VERBATIM("MakeAttenuationProfiles") << "fFiberBrightnessMode is off";
228  }
229 
230  //......................................................................
232  {
235 
236  // The rest is making histograms. Don't do it again on a subsequent run
237  if(fdx) return;
238 
239  const TString viewStrArr[4]={"X","Y","muX","muY"};
241  bool isND = fGeom->DetId() == novadaq::cnv::kNEARDET;
242 
244  unsigned int nbrightnesses = 1;
245  if(fFiberBrightnessMode) nbrightnesses = fbhandle->NumberBrightnessBins();
246 
247  if(fFiberBrightnessBins > nbrightnesses){
248  std::cerr << "Asking for too many fiberbrightness bins! FB Service has : " << nbrightnesses << " while you're asking for : " << fFiberBrightnessBins << std::endl;
249  std::abort();
250  }
251  if(fFiberBrightnessBins > 12){
252  std::cerr << "Asking for too many fiberbrightness bins! Maximum is 12" << std::endl;
253  std::abort();
254  }
255 
256  LOG_VERBATIM("MakeAttenuationProfiles") << "Number of fiber brightness bins : " << fFiberBrightnessBins;
257 
258  double maxW;
259  // FD
260  int maxCell = 390;
261  int maxPlane = 900;
263  maxW = 900;
264  }
265  else if(fGeom->DetId() == novadaq::cnv::kNEARDET){
266  maxW = 250;
267  maxPlane = 220;
268  maxCell = 100;
269  }
270  else if(fGeom->DetId() == novadaq::cnv::kTESTBEAM){
271  maxW = 150;
272  }
273  else{
274  // NDOS (hopefully)
275  maxW = 250;
276  }
277 
278  fPlanevsCell = tfs->make<TH2F>("h_Plane_vs_Cell",
279  ";Plane;Cell",
280  maxPlane, 0, maxPlane,
281  maxCell, 0, maxCell);
282 
283  for(unsigned int fiberbrightness = 0; fiberbrightness < fFiberBrightnessBins; fiberbrightness++)
284  {
285  for(int view = geo::kX; view <= (geo::kY+2*isND); ++view){
286  TString viewStr = viewStrArr[view];//was const
287  TString fbStr = "";
288  if(view > geo::kY && fiberbrightness > 0 && !fMuCFiberBrightnessMode) continue; // don't want it for muon catcher unless asked for
290  fbStr.Form("_fb%i", fiberbrightness);
291  }
292 
293 
295  v.WPE = tfs->make<TH2F>("h_WPE_in"+viewStr+fbStr,
296  ";W;PE",
297  2*maxW, -maxW, +maxW,
298  1000, 0, 1000);
299  v.WPE_corr = tfs->make<TH2F>("h_WPE_corr_avg_in"+viewStr+fbStr,
300  ";W;PE/cm",
301  2*maxW, -maxW, +maxW,
302  1000, 0, 200);
303  v.WPE_corr_xy = tfs->make<TH2F>("h_WPE_corr_xy_in"+viewStr+fbStr,
304  ";W;PE/cm",
305  2*maxW, -maxW, +maxW,
306  1000, 0, 200);
307  v.WPE_corr_xy_truth = tfs->make<TH2F>("h_WPE_corr_xy_truth_in"+viewStr+fbStr,
308  ";W;MeV/cm",
309  2*maxW, -maxW, +maxW,
310  1000, 0, 10);
311  v.WPE_corr_z = tfs->make<TH2F>("h_WPE_corr_z_in"+viewStr+fbStr,
312  ";W;PE/cm",
313  2*maxW, -maxW, +maxW,
314  1000, 0, 200);
315  v.WPE_corr_traj = tfs->make<TH2F>("h_WPE_corr_traj_in"+viewStr+fbStr,
316  ";W;PE/cm",
317  2*maxW, -maxW, +maxW,
318  1000, 0, 200);
319  fSumChannels[fiberbrightness][view] = v;
320 
321  fCellPE[fiberbrightness][view] = tfs->make<TH2F>("h_CellPE_in"+viewStr+fbStr, ";Cell;PE/cm",
322  maxCell+5, -5, maxCell, 1000, 0, 200);
323  fCellRawPE[fiberbrightness][view] = tfs->make<TH2F>("h_CellRawPE_in"+viewStr+fbStr, ";Cell;PE",
324  maxCell+5, -5, maxCell, 1000, 0, 1000);
325  fPlanePE[fiberbrightness][view] = tfs->make<TH2F>("h_PlanePE_in"+viewStr+fbStr, ";Plane;PE",
326  maxPlane+5, -5, maxPlane, 1000, 0, 200);
327 
328  fPathLength[fiberbrightness][view] = tfs->make<TH1F>("h_PathLength"+viewStr+fbStr,
329  ";Path length (cm)",
330  80, 0, 20);
331 
332  if(fiberbrightness == 0)
333  {
334  fBrightnessPE[view] = tfs->make<TH2F>("h_BrightnessPE_in"+viewStr+fbStr, ";Fiber Brightness;PE",
335  12, 0, 12, 1000, 0, 1000);
336  fBrightnessProf[view] = tfs->make<TProfile>("p_BrightnessPE_in"+viewStr+fbStr, ";Fiber Brightness;PE",
337  12, 0, 12, 0, 1000);
338  fBrightnessCanv[view] = tfs->make<TCanvas>("c_BrightnessPE_in"+viewStr+fbStr);
339  }
340  } // end for view
341  } // end for fiberbrightness
342 
343  fdx = tfs->make<TH1F>("h_dx", ";|dx|", 100, 0, 1);
344  fdy = tfs->make<TH1F>("h_dy", ";|dy|", 100, 0, 1);
345  fdz = tfs->make<TH1F>("h_dz", ";|dz|", 100, 0, 1);
346 
347  fdxdy = tfs->make<TH2F>("h_dxdy", ";|dx|;|dy|", 100, 0, 1, 100, 0, 1);
348  fdxdz = tfs->make<TH2F>("h_dxdz", ";|dx|;|dz|", 100, 0, 1, 100, 0, 1);
349  fdydz = tfs->make<TH2F>("h_dydz", ";|dy|;|dz|", 100, 0, 1, 100, 0, 1);
350 
351  fStartXY = tfs->make<TH2F>("h_StartXY", ";x;y",
352  100, -200, +200, 100, -250, +250);
353  fStartZX = tfs->make<TH2F>("h_StartZX", ";z;x",
354  100, -50, +900, 100, -200, +200);
355  fStartZY = tfs->make<TH2F>("h_StartZY", ";z;y",
356  100, -50, +900, 100, -250, +250);
357  fStopXY = tfs->make<TH2F>("h_StopXY", ";x;y",
358  100, -200, +200, 100, -250, +250);
359  fStopZX = tfs->make<TH2F>("h_StopZX", ";z;x",
360  100, -50, +900, 100, -200, +200);
361  fStopZY = tfs->make<TH2F>("h_StopZY", ";z;y",
362  100, -50, +900, 100, -250, +250);
363 
364  fBrightnessHist = tfs->make<TH1F>("h_Brightness", ";Brightness;",
365  12, 0, 12);
366 
367  return;
368  }
369 
370  //......................................................................
372  {
373  if(fProfileMode){
374  for(auto label : fInstanceLabels){
375  auto it = fChannelMapProf.find(label.first);
376  if(it == fChannelMapProf.end()) continue;
377 
378  if(fHistoOutput){
379  // Accumulate a running total for the histogram file
380  auto ittot = fChannelMapProfTotal.find(label.first);
381  if(ittot == fChannelMapProfTotal.end())
382  fChannelMapProfTotal.insert(std::make_pair(label.first, *it->second));
383  else
384  ittot->second += *it->second;
385  }
386 
387  r.put(std::move(it->second), label.second);
388  }
389 
390  }
391  else{
392  std::unique_ptr<caldp::AttenHistsMap> cm(new caldp::AttenHistsMap(fChannelMap));
393  r.put(std::move(cm));
394 
395  // Accumulate a running total for the histogram file
397  // Don't double count for the output products
398  }
399 
400  return;
401  }
402 
403  //......................................................................
405  {
406 
407  // handle the case where we want to use the raw ADC values rather than
408  // the PCHit PE value for finding the attenuation curves
409  if(fPECalculation == 1){
410  // get the RawDigit list
411  auto digitcol = evt.getValidHandle< std::vector<rawdata::RawDigit> >(fRawDataLabel);
412 
413  // make a map of the CellUniqueId to rawdigit
414  size_t plane = 0;
415  size_t cell = 0;
416  geo::CellUniqueId cellId = 0;
417  for(auto const& rawdigit : *digitcol){
418  plane = fCMap->GetPlane(&rawdigit);
419  cell = fCMap->GetCell(&rawdigit);
420  cellId = fGeom->Plane(plane)->Cell(cell)->Id();
421  fCellIdToRawDigit[cellId] = rawdigit;
422  }
423  }// end if using ADCs instead of PEs
424 
425  auto pchit_avg = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualAvg" ));
426  auto pchit_xy = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualXY" ));
427  auto pchit_xy_thresh = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualXYThresh"));
428  auto pchit_z = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualZ" ));
429  auto pchit_traj = evt.getValidHandle< std::vector< caldp::PCHit > >(art::InputTag(fPCHitsLabel, "pathQualTraj" ));
430 
431  this->FillPathLengthInfo(*pchit_z, caldp::kCalZ, evt.isRealData());
432  this->FillPathLengthInfo(*pchit_traj, caldp::kCalTraj, evt.isRealData());
433  this->FillPathLengthInfo(*pchit_xy, caldp::kCalXY, evt.isRealData());
434  this->FillPathLengthInfo(*pchit_avg, caldp::kCalAvg, evt.isRealData());
437  }
438 
439  return;
440 
441  }
442 
443  //......................................................................
444  void MakeAttenuationProfiles::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 
453 
454  for(auto const& current_hit : pchits){
455 
456  double w = current_hit.W();
457  double path = current_hit.Path();
458  double pe = this->PCHitPE(current_hit);
459  double mev = current_hit.TrueMeV();
460  int plane = current_hit.Plane();
461  int cell = current_hit.Cell();
462  float time = current_hit.TNS();
463  const int view = this->GetView(plane);
464  int fiberbrightness = 0;
465  if((fFiberBrightnessMode && view <= geo::kY) || (fMuCFiberBrightnessMode && view > geo::kY))
466  fiberbrightness = fbhandle->getBrightnessIndex(plane,cell);//GetFiberBrightness(plane, cell);
467 
468  if(useBelowThreshold && pe != 0.)
469  throw cet::exception("CosmicCalib")
470  << "attempting to use below threshold hit with nonzero pe: " << pe;
471 
473  fBrightnessPE[view]->Fill(fiberbrightness, pe);
474  fBrightnessProf[view]->Fill(fiberbrightness, pe);
475  }
476 
477  if(path > 10) continue;
478 
479  // Beam Spills weren't removed in the miniproduction files, apply Tyler's rough cut instead (docdb : 37413)
480  if(fCutBeamSpills && (time < -240. || time > 1000.)) continue;
481 
482  if(pathType == caldp::kCalAvg){
483  }
484  if(pathType == caldp::kCalTraj && path > fTrajectory_pathLow && path < fTrajectory_pathHigh){
485  fSumChannels[fiberbrightness][view].HistogramByPathType(pathType)->Fill(w, pe/path);
486  }
487  else{
488  if(pathType == caldp::kCalXY){
489  fCellPE[fiberbrightness][view]->Fill(cell, pe/path);
490  fCellRawPE[fiberbrightness][view]->Fill(cell, pe);
491  fPlanePE[fiberbrightness][view]->Fill(plane, pe/path);
492  fPathLength[fiberbrightness][view]->Fill(path);
493  }
494  fSumChannels[fiberbrightness][view].HistogramByPathType(pathType)->Fill(w, pe/path);
495  if(pathType == caldp::kCalXY && !useBelowThreshold)
496  fSumChannels[fiberbrightness][view].HistogramByPathType(caldp::kCalXYTruth)->Fill(w, mev/path);
497  }
498 
499  if(useProfiles){
500  caldp::AttenProfiles& profs = GetChannelProfiles(fGeom->DetId(), isRealData, plane, cell);
501  if(pathType == caldp::kCalTraj && path > fTrajectory_pathLow && path < fTrajectory_pathHigh){
502  profs.LiteProfileByPathType(pathType).Fill(w, pe/path);
503  }//pathlength 4 cm condition
504  else{
505  profs.LiteProfileByPathType(pathType).Fill(w, pe/path);
506  if(pathType == caldp::kCalXY && !useBelowThreshold)
507  profs.LiteProfileByPathType(caldp::kCalXYTruth).Fill(w, mev/path);
508  }
509  }
510  else if(useHists){
511  caldp::AttenHists hists = GetChannelHists(isRealData, plane, cell);
512  hists.HistogramByPathType(pathType)->Fill(w, pe/path);
513  if(pathType == caldp::kCalTraj && path > fTrajectory_pathLow && path < fTrajectory_pathHigh){
514  hists.HistogramByPathType(pathType)->Fill(w, pe/path);
515  }//pathlength 4 cm condition
516  if(pathType == caldp::kCalAvg)
517  hists.HistogramByPathType(caldp::kCalUnknown)->Fill(w, pe/path);
518  else if(pathType == caldp::kCalXY && !useBelowThreshold)
519  hists.HistogramByPathType(caldp::kCalXYTruth)->Fill(w, mev/path);
520  }
521 
522  if(pathType == caldp::kCalXY) fPlanevsCell->Fill(plane, cell);
523 
524  }
525 
526  return;
527  }
528 
529  //......................................................................
531  {
532  // this method provides multiple ways of getting a "PE" value
533  // and putting that into the plots used to make the attenuation
534  // curves. What value is returned is based on a parameter set at
535  // run time
536  double pe = 0.;
537 
538  if(fPECalculation == 0) return pchit.PE();
539  else if(fPECalculation == 1){
540  // get the cellId for the PCHit
541  auto pchitId = fGeom->Plane(pchit.Plane())->Cell(pchit.Cell())->Id();
542 
543  if( fCellIdToRawDigit.count(pchitId) > 0) pe = 1.*fCellIdToRawDigit.find(pchitId)->second.ADC();
544  else
545  throw cet::exception("CosmicCalib")
546  << "Cannot find RawDigit for "
547  << pchit.Plane()
548  << " " << pchit.Cell();
549 
550  }// end if using raw digit ADC values
551 
552  return pe;
553  }
554 
555  //......................................................................
557  int const& plane,
558  int const& cell)
559  {
560  // plane is a const reference so we can't change it.
561  // instead make a new variable called pln that we can change.
562  // if we are consolidating the views, set pln to the view
563  int pln = plane;
564  const int view = GetView(plane);
566  pln = GetView(plane);
567  }
568 
570 
571  // If the histograms are created for this channel, just return them
572  if(hists.WPE) return hists;
573 
574  // Otherwise have to construct them
575  const double maxW = fGeom->Plane(plane)->Cell(0)->HalfL()*1.15;
576 
577  TH2F* wadc = new TH2F(TString::Format("h_WPE_%03d_%03d", pln, cell),
578  ";W;PE", 100, -maxW, +maxW, 100, 0, 1000);
579 
580  TH2F* wadc_corr = new TH2F(TString::Format("h_WPE_corr_avg_%03d_%03d",
581  pln, cell),
582  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
583 
584  TH2F* wadc_corr_xy = new TH2F(TString::Format("h_WPE_corr_xy_%03d_%03d",
585  pln, cell),
586  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
587 
588  TH2F* wadc_corr_xy_truth = new TH2F(TString::Format("h_WPE_corr_xy_truth_%03d_%03d", pln, cell),
589  ";W;MeV/cm", 100, -maxW, +maxW, 100, 0, 10);
590 
591 
592  TH2F* wadc_corr_z = new TH2F(TString::Format("h_WPE_corr_z_%03d_%03d", pln, cell),
593  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
594 
595  TH2F* wadc_corr_traj = new TH2F(TString::Format("h_WPE_corr_traj_%03d_%03d", pln, cell),
596  ";W;PE/cm", 100, -maxW, +maxW, 100, 0, 200);
597 
598  hists.WPE = wadc;
599  hists.WPE_corr = wadc_corr;
600  hists.WPE_corr_xy = wadc_corr_xy;
601  hists.WPE_corr_xy_truth = wadc_corr_xy_truth;
602  hists.WPE_corr_z = wadc_corr_z;
603  hists.WPE_corr_traj = wadc_corr_traj;
604 
605  return hists;
606  }
607 
608  //......................................................................
610  bool const& isData,
611  int const& plane,
612  int const& cell)
613  {
614  // MC only varies for these reasons. Consolidate due to low stats.
615  // Make a new variable, pln to be either the plane or the view if
616  // consolidating
617  int pln = plane;
618  const int view = GetView(plane);
620  pln = GetView(plane);
621  }
622 
623  double maxW;
624  if(det == novadaq::cnv::kFARDET) maxW = 800;
625  else if(det == novadaq::cnv::kNEARDET) maxW = 250;
626  else maxW = 250; // Don't bother to make non-square NDOS profiles
627 
628  if((fFiberBrightnessMode && view <= geo::kY) || (fMuCFiberBrightnessMode && view > geo::kY)){
630  int fb = fbhandle->getBrightnessIndex(plane, cell);
631  if(fChannelMapProf.find(fb) == fChannelMapProf.end())
632  fChannelMapProf.emplace(std::make_pair(fb, std::unique_ptr<caldp::AttenProfilesMap>(new caldp::AttenProfilesMap(-maxW, +maxW))));
633  return fChannelMapProf.find(fb)->second->GetProfiles(pln, cell);
634  } else {
635  int block = pln/32;
636  if(view > geo::kY) block = 100; // special case
637  if(fChannelMapProf.find(block) == fChannelMapProf.end())
638  fChannelMapProf.emplace(std::make_pair(block, std::unique_ptr<caldp::AttenProfilesMap>(new caldp::AttenProfilesMap(-maxW, +maxW))));
639  return fChannelMapProf.find(block)->second->GetProfiles(pln, cell);
640  }
641  }
642 
643  //......................................................................
645  //get the view for the current combination of view and plane (x,y,mux, or muy)
646  int view=0;
648  int firstPlane = fGeom->FirstPlaneInMuonCatcher();
649  int isMuonCatcher = firstPlane <= plane;
650  view = fGeom->Plane(plane)->View()+(2*isMuonCatcher);
651  }
652  else{
653  view = fGeom->Plane(plane)->View();
654  }
655  return view;
656  }
657 
658  //......................................................................
660  {
661 
662  if(!fHistoOutput) return;
663 
665 
666  //LOG_VERBATIM("MakeAttenuationProfiles") << "beginning of the endJob";
667 
668  for(auto& it: fChannelMapTotal.GetAllHistsByPlane()){
669  const int plane = it.first;
670  const std::vector<TH2F*>& hs = it.second;
671 
672  art::TFileDirectory dir = tfs->mkdir(Form("plane_%03d", plane));
673 
674  for(TH2F* h: hs) dir.make<TH2F>(*h);
675  }
676 
677  for(int block = 0; block < 28; ++block){
678  auto itmap = fChannelMapProfTotal.find(block);
679  if(itmap == fChannelMapProfTotal.end()) continue;
680  for(auto& it: itmap->second.GetAllProfilesByPlaneAndCell()){
681  const int plane = it.first;
682  std::cout << plane << std::endl;
683  if(it.second.empty()) continue;
684 
685  art::TFileDirectory dirBlock = tfs->mkdir(Form("block_%02d", plane/32));
686  art::TFileDirectory dirPlane = dirBlock.mkdir(Form("plane_%03d", plane));
687 
688  for(auto& it2: it.second){
689  const int cell = it2.first;
690  const std::vector<TH1F*>& hs = it2.second;
691  if(hs.empty()) continue;
692 
693  art::TFileDirectory dirMod = dirPlane.mkdir(Form("module_%02d", cell/32));
694  art::TFileDirectory dirCell = dirMod.mkdir(Form("cell_%03d", cell));
695 
696  for(TH1F* h: hs) dirCell.make<TH1F>(*h);
697  }
698  }
699  }
700  }
701 
703 
704 } // end namespace
705 ////////////////////////////////////////////////////////////////////////
std::map< int, caldp::AttenProfilesMap > fChannelMapProfTotal
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
virtual void reconfigure(const fhicl::ParameterSet &pset)
caldp::AttenHistsMap fChannelMapTotal
This is what we put in the histogram file.
set< int >::iterator it
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
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
art::ProductID put(std::unique_ptr< PROD > &&product)
Definition: Results.h:78
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
bool isRealData() const
Definition: Event.h:83
TString hists[nhists]
Definition: bdt_com.C:3
TH2F * WPE_corr_xy_truth
Definition: AttenHists.h:67
unsigned int NumberBrightnessBins() const
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 char * label
art::ServiceHandle< geo::Geometry > fGeom
const CellUniqueId & Id() const
Definition: CellGeo.h:55
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
art::ServiceHandle< cmap::CMap > fCMap
connection map
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
CDPStorage service.
block
print "ROW IS " print row
Definition: elec2geo.py:31
virtual void writeResults(art::Results &r)
float PE() const
Return PE value.
Definition: PCHit.h:38
enum caldp::_path_type PathType_t
#define DEFINE_ART_RESULTS_PLUGIN(klass)
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
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
Histograms used by attenuation calibration.
Definition: AttenHists.h:27
double PCHitPE(caldp::PCHit const &pchit)
virtual void beginRun(art::Run const &run)
double DetHalfHeight() const
TH2F * HistogramByPathType(caldp::PathType_t const &pt)
Definition: AttenHists.h:52
int Plane() const
Return plane value.
Definition: PCHit.h:24
int fPECalculation
0 = PCHit::PE(), 1 = DCS ADC value from raw digit
AttenHists for many channels.
Definition: AttenHists.h:77
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
Profiles used by attenuation calibration.
Definition: AttenProfiles.h:61
const std::string path
Definition: plot_BEN.C:43
virtual void event(art::Event const &evt)
caldp::AttenHists GetChannelHists(bool const &isData, int const &plane, int const &cell)
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
caldp::AttenHistsMap fChannelMap
This is where GetChannelHist stores its histograms. Map is from block #.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
LiteProfile & LiteProfileByPathType(caldp::PathType_t const &pt)
Definition: AttenProfiles.h:71
double DetHalfWidth() const
void FillPathLengthInfo(std::vector< caldp::PCHit > const &pchits, caldp::PathType_t const &pathType, bool const &isRealData, bool const &useBelowThreshold=false)
TDirectory * dir
Definition: macro.C:5
std::list< std::pair< int, std::string > > fInstanceLabels
labels for the different products
caldp::AttenProfiles & GetChannelProfiles(novadaq::cnv::DetId const &det, bool const &isData, int const &plane, int const &cell)
double fTrajectory_pathLow
greater than fTrajectory_pathLow (cm).
double fTrajectory_pathHigh
less than fTrajectory_pathHigh (cm).
std::map< int, std::unique_ptr< caldp::AttenProfilesMap > > fChannelMapProf
std::map< geo::CellUniqueId, rawdata::RawDigit > fCellIdToRawDigit
map cell id to the raw digit
TRandom3 r(0)
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
#define LOG_VERBATIM(category)
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
MakeAttenuationProfiles(const fhicl::ParameterSet &pset)
Float_t w
Definition: plot.C:20
TH2F * WPE_corr_traj
Definition: AttenHists.h:67
static constexpr Double_t cm
Definition: Munits.h:140
Encapsulate the geometry of one entire detector (near, far, ndos)
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.