ThresholdAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ThresholdAna_module.cc
3 // \brief Make plots of threshold effect from PCHits files
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 //C++ and C:
8 #include <cmath>
9 #include <vector>
10 
11 //Framework:
16 
17 //NovaSoft:
20 #include "Geometry/Geometry.h"
22 
23 // ROOT
24 #include "TH2.h"
25 
26 namespace calib
27 {
28  /// Make plots of threshold and shielding effects from PCHits files
30  {
31  public:
32  explicit ThresholdAna(const fhicl::ParameterSet& pset);
33  virtual ~ThresholdAna();
34 
35  virtual void beginRun(const art::Run& run) override;
36 
37  virtual void analyze(const art::Event& evt) override;
38 
39  virtual void reconfigure(const fhicl::ParameterSet& pset);
40 
41  protected:
43  //int GetFiberBrightness(int const& plane, int const& cell);
44  TH1* fNCalibVsW[4];//one per view (4)
45  TH1* fNZerosVsW[4];
46 
47  TH1* fNCalibVsCell[4];
48  TH1* fNZerosVsCell[4];
49 
52 
53  TH2* fPEpercmVsWTot[4];
54  TH2* fRecoVsWTot[4];
55  TH2* fTrueVsWTot[4];
56  TH2* fThreshVsWTot[4];
57  TH2* fShadowVsWTot[4];
58  TH2* fRatioVsWTot[4];
59 
60  TH2* fPEpercmVsW[4][32*12];
61  TH2* fRecoVsW[4][32*12];//one per view, per cell
62  TH2* fTrueVsW[4][32*12];
63  TH2* fThreshVsW[4][32*12];
64  TH2* fShadowVsW[4][32*12];
65  TH2* fRatioVsW[4][32*12];
66 
67  TH2* fPEpercmVsMeVpercmTot[4];//one per view
68  TH2* fPEpercmVsMeVpercm[4][32*12];//one per view, per cell
69 
71 
76  };
77 
78  //---------------------------------------------------------------------
80  : EDAnalyzer(pset)
81  {
82  reconfigure(pset);
83 
84  // Flag that beginRun hasn't run yet
85  fNZerosVsW[0] = 0;
86  }
87 
88  //---------------------------------------------------------------------
90  {
91  }
92 
93  //---------------------------------------------------------------------
95  {
96  fPathMode = pset.get<std::string>("PathMode", "pathQualXY");
97  fFiberBrightnessMode = pset.get<bool> ("FiberBrightnessMode", false);
98  fMuCFiberBrightnessMode = pset.get<bool> ("MuCFiberBrightnessMode", false);
99  fFiberBrightness = pset.get<int> ("FiberBrightness", 0);
100 
101  if(fFiberBrightnessMode)LOG_VERBATIM("ThresholdAna") << "Looking at FB " << fFiberBrightness;
102  else LOG_VERBATIM("ThresholdAna") << "FiberBrightnessMode is OFF, fFiberBrightness = " << fFiberBrightness;
103 
104  }
105 
106  //---------------------------------------------------------------------
108  {
109  //LOG_VERBATIM("ThresholdAna") << "beginRun\n";
110  //LOG_VERBATIM("ThresholdAna") << fNZerosVsW[0] << " should be 0\n";
111 
112  // We've already been called for a previous run, no need to make histograms
113  // again.
114  if(fNZerosVsW[0]){
115  //LOG_VERBATIM("ThresholdAna") << "Already made histograms (fNZerosVsW[0])\n";
116  return;
117  }
118 
119  //LOG_VERBATIM("ThresholdAna") << "Going to make histograms (!fNZerosVsW[0])\n";
120 
122  //LOG_VERBATIM("ThresholdAna") << "Got geom\n";
123 
125  int nbrightnesses = 1;
126  if(fFiberBrightnessMode || fMuCFiberBrightnessMode) nbrightnesses = fbhandle->NumberBrightnessBins();
127  if(fFiberBrightness >= nbrightnesses){
128  std::cerr << "Asking for too many fiberbrightness bins! FB Service has : " << nbrightnesses << " while you're asking for : " << fFiberBrightness << std::endl;
129  std::abort();
130  }
131  if(fFiberBrightness > 11){
132  std::cerr << "Asking for too many fiberbrightness bins! Maximum is 12" << std::endl;
133  std::abort();
134  }
136  std::cout << "WARNING!!!" << std::endl;
137  std::cout << "Asking for fiber brightness only in muon catcher but not the active detector, Why??" << std::endl;
138  }
139 
140  double maxW = -1;
141  int nCell = -1;
142  if(geom->DetId() == novadaq::cnv::kFARDET){
143  maxW = 900;
144  nCell = 32*12;
145  //LOG_VERBATIM("ThresholdAna") << "FD: maxW = 900, nCell = 32*12";
146  }
147  else if(geom->DetId() == novadaq::cnv::kNEARDET ||
148  geom->DetId() == novadaq::cnv::kNDOS){
149  maxW = 250;
150  nCell = 32*3;
151  //LOG_VERBATIM("ThresholdAna") << "ND: maxW = 250, nCell = 32*3";
152  }
153  else if(geom->DetId() == novadaq::cnv::kTESTBEAM){
154  maxW = 150;
155  nCell = 32*2;
156  //LOG_VERBATIM("ThresholdAna") << "TB: maxW = 150, nCell = 32*3";
157  }
158  else{
159  assert(0 && "Unknown detector");
160  LOG_VERBATIM("ThresholdAna") << "What detector??";
161  }
162 
164  //LOG_VERBATIM("ThresholdAna") << "Started TFileService\n";
165 
166  const int nW = 100;
167 
168  const int nY = 1000;
169 
170  const TString viewStrArr[4]={"_xview","_yview","_muxview","_muyview"};
171  bool isND=geom->DetId() == novadaq::cnv::kNEARDET;
172 
173  //LOG_VERBATIM("ThresholdAna") << "Making histograms\n";
174  for(int view = geo::kX; view <= (geo::kY+2*isND); ++view){
175 
176  TString viewStr = viewStrArr[view];
177 
178  if(view > geo::kY && fFiberBrightness > 0 && !fMuCFiberBrightnessMode) continue; // don't want it for muon catcher unless asked for
179 
180  //LOG_VERBATIM("ThresholdAna") << "View: " << viewStrArr[view];
181 
182  fNCalibVsBrightness[view] = tfs->make<TH1F>("nCalibVsBrightness"+viewStr, ";Fiber Brightness;Number of PCHits", nbrightnesses, 0, nbrightnesses);
183  fNZerosVsBrightness[view] = tfs->make<TH1F>("nZerosVsBrightness"+viewStr, ";Fiber Brightness;Number of PCHits below threshold", nbrightnesses, 0, nbrightnesses);
184 
186  viewStr = TString::Format("%s_fb%i", viewStr.Data(), fFiberBrightness);
187  //LOG_VERBATIM("ThresholdAna") << "View = " << viewStrArr[view] << "\t fibre brightness = " << fFiberBrightness;
188  }
189 
190  fNZerosVsW[view] = tfs->make<TH1F>("nZerosVsW"+viewStr, ";W;Number of PCHits below threshold", nW, -maxW, +maxW);//";W;"
191  fNCalibVsW[view] = tfs->make<TH1F>("nCalibVsW"+viewStr, ";W;Number of PCHits", nW, -maxW, +maxW);
192 
193  fNZerosVsCell[view] = tfs->make<TH1F>("nZerosVsCell"+viewStr, ";Cell;Number of PCHits below threshold", nCell, 0, nCell);//";Cell"
194  fNCalibVsCell[view] = tfs->make<TH1F>("nCalibVsCell"+viewStr, ";Cell;Number of PCHits", nCell, 0, nCell);
195 
196  fPEpercmVsWTot[view] = tfs->make<TH2F>(TString::Format("fPEpercmVsW%s", viewStr.Data()),"PE/cm over W;W (cm);PE/cm", nW, -maxW, +maxW, 200, 0 , 200);
197 
198  fRatioVsWTot[view] = tfs->make<TH2F>(TString::Format("ratioVsW%s_tot", viewStr.Data()), "Total threshold and shadow correction;W (cm);PE / PE poisson #lambda * E_{true} (MeV) / PE_{MIP}^{exp}", nW, -maxW, +maxW, nY, 0, 10);//";W (cm);Reco / Ideal"
199 
200  fRecoVsWTot[view] = tfs->make<TH2F>(TString::Format("recoVsW%s_tot", viewStr.Data()), ";W (cm);PE / PE_{MIP}^{exp}", nW, -maxW, +maxW, nY, 0, 100);//alt titles: ";W (cm);PE / MIP expectation for path" ";W (cm);Reco = Number of PE deposited / Number of PE expected for path for 1 MIP"
201  fTrueVsWTot[view] = tfs->make<TH2F>(TString::Format("trueVsW%s_tot", viewStr.Data()), ";W (cm);PE poisson #lambda / E_{true} (MeV)", nW, -maxW, +maxW, nY, 0, 50);//alt title: ";W (cm);PE poisson lambda / true E deposit (MeV)"
202 
203  fThreshVsWTot[view] = tfs->make<TH2F>(TString::Format("threshVsW%s_tot", viewStr.Data()), "Threshold correction;W (cm);PE / PE poisson #lambda", nW, -maxW, +maxW, nY, 0, 10);//";W (cm);PE / poisson lambda"
204  fShadowVsWTot[view] = tfs->make<TH2F>(TString::Format("shadowVsW%s_tot", viewStr.Data()), "Shadow correction;W (cm);E_{true} (MeV) / PE_{MIP}^{exp}", nW, -maxW, +maxW, nY, 0, 10);//";W (cm);True E / MIP expectation"
205 
206  fPEpercmVsMeVpercmTot[view] = tfs->make<TH2F>(TString::Format("PEpercmVsMeVpercm%s_tot", viewStr.Data()), ";PE/cm;MeV/cm", 100, 0, 10, 100, 0, 100);
207 
208  for(int cell = 0; cell < nCell; ++cell){
209  //LOG_VERBATIM("ThresholdAna") << "\t cell = " << cell;
210  fPEpercmVsW[view][cell] = tfs->make<TH2F>(TString::Format("fPEpercmVsW%s_cell%d", viewStr.Data(), cell),"PE/cm over W;W (cm);PE/cm", nW, -maxW, +maxW, 200, 0 ,200);
211 
212  fRatioVsW[view][cell] = tfs->make<TH2F>(TString::Format("ratioVsW%s_cell%d", viewStr.Data(), cell), "Total threshold and shadow correction;W (cm);PE / PE poisson #lambda * E_{true} (MeV) / PE_{MIP}^{exp}", nW, -maxW, +maxW, nY, 0, 10);
213 
214  fRecoVsW[view][cell] = tfs->make<TH2F>(TString::Format("recoVsW%s_cell%d", viewStr.Data(), cell), ";W (cm);PE / PE_{MIP}^{exp}", nW, -maxW, +maxW, nY, 0, 100);
215  fTrueVsW[view][cell] = tfs->make<TH2F>(TString::Format("trueVsW%s_cell%d", viewStr.Data(), cell), ";W (cm);PE poisson #lambda / E_{true} (MeV)", nW, -maxW, +maxW, nY, 0, 50);
216 
217  fThreshVsW[view][cell] = tfs->make<TH2F>(TString::Format("threshVsW%s_cell%d", viewStr.Data(), cell), "Threshold correction;W (cm);PE / PE poisson #lambda", nW, -maxW, +maxW, nY, 0, 10);
218  fShadowVsW[view][cell] = tfs->make<TH2F>(TString::Format("shadowVsW%s_cell%d", viewStr.Data(), cell), "Shadow correction;W (cm);E_{true} (MeV) / PE_{MIP}^{exp}", nW, -maxW, +maxW, nY, 0, 10);
219 
220  fPEpercmVsMeVpercm[view][cell] = tfs->make<TH2F>(TString::Format("PEpercmVsMeVpercm%s_cell%d", viewStr.Data(), cell), ";PE/cm;MeV/cm", 100, 0, 10, 100, 0, 100);
221 
222  } // end for cell
223 
224  } // end for view
225  LOG_VERBATIM("ThresholdAna") << "Ending beginRun\n";
226  }
227 
228  //......................................................................
230  //get the view for the current combination of view and plane (x,y,mux, or muy)
231  //LOG_VERBATIM("ThresholdAna") << "Running GetView...";
232  int view=0;
233  if(geom->DetId() == novadaq::cnv::kNEARDET){
234  int firstPlane=geom->FirstPlaneInMuonCatcher();
235  int isMuonCatcher=firstPlane <= plane;
236  view = geom->Plane(plane)->View()+(2*isMuonCatcher);
237  }
238  else{
239  view = geom->Plane(plane)->View();
240  }
241  //LOG_VERBATIM("ThresholdAna") << "Successful\n";
242  return view;
243  }
244 
245  //---------------------------------------------------------------------
247  {
248  //LOG_VERBATIM("ThresholdAna") << "Running analyze...\n";
252 
253  // LOG_VERBATIM("ThresholdAna") << "PathType = " << fPathMode << std::endl;
254 
255  evt.getByLabel("pclist", fPathMode, hits);
256 
257  //LOG_DEBUG("ThresholdAna") << hits->size() << " hits";
258 
259  for(unsigned int n = 0; n < hits->size(); ++n){
260  const caldp::PCHit& hit = (*hits)[n];
261  int plane = hit.Plane();
262  int view = GetView(plane,geom);
263  int brightness;
264  if((fFiberBrightnessMode && view <= geo::kY) || (fMuCFiberBrightnessMode && view > geo::kY))
265  brightness = fbhandle->getBrightnessIndex(plane, hit.Cell());
266  else brightness = 0;
267 
268  // Only make the n calib vs fb plot for brightness=0
270 
271  // Only fill if it's the brightness that we're looking at
272  if((fFiberBrightnessMode || fMuCFiberBrightnessMode) && brightness != fFiberBrightness) continue;
273 
274  //LOG_VERBATIM("ThresholdAna") << "index = " << index;
275 
276  // Saturated, skip
277  if(hit.PE() > 2000){
278  LOG_VERBATIM("ThresholdAna") << "Saturated, skip";
279  continue;
280  }
281 
282  fNCalibVsW[view]->Fill(hit.W());
283  fNCalibVsCell[view]->Fill(hit.Cell());
284 
285  // This is the expected number of photons produced by 1 MIP. No need for
286  // high number of decimal places, the absolute calibration actually sorts
287  // out the MIP scale.
288  const double mipExp = 1.78*hit.Path();
289 
290  // This is essentially what the actual attenuation fit uses (except the
291  // factor to convert PEs to MIPs).
292  const double recoPart = hit.PE()/mipExp;
293 
294  // This is what we actually want to measure, the true transmission
295  // function.
296  const double truePart = hit.PoissonLambda()/hit.TrueMeV();
297 
298  // Correction due to threshold effects: that we get a biased-high number
299  // of PEs because of fluctuations sometimes needed to get over threshold.
300  const double threshPart = hit.PE()/hit.PoissonLambda();
301 
302  // Correction due to "shadowing" (that the true energy deposit isn't
303  // equally MIPish everywhere).
304  const double shadowPart = hit.TrueMeV()/mipExp;
305 
306  // Written this way, it's clear that applying this ratio to the
307  // observations will correct them to be the attenuation function we
308  // actually want to measure. Could also have written as
309  // threshPart*shadowPart, which makes the physical origins of the
310  // corrections clearer.
311  const double ratio = recoPart/truePart;
312 
313  // We fill all the different combinations, because when the distributions
314  // are broad, different ways of taking the mean, eg <reco/true> vs
315  // <reco>/<true> make a difference.
316 
317  fPEpercmVsWTot[view]->Fill(hit.W(), hit.PE() / hit.Path());
318  fRecoVsWTot [view]->Fill(hit.W(), recoPart);
319  fTrueVsWTot [view]->Fill(hit.W(), truePart);
320  fThreshVsWTot [view]->Fill(hit.W(), threshPart);
321  fShadowVsWTot [view]->Fill(hit.W(), shadowPart);
322  fRatioVsWTot [view]->Fill(hit.W(), ratio);
323 
324  fPEpercmVsW[view][hit.Cell()]->Fill(hit.W(), hit.PE() / hit.Path());
325  fRecoVsW [view][hit.Cell()]->Fill(hit.W(), recoPart);
326  fTrueVsW [view][hit.Cell()]->Fill(hit.W(), truePart);
327  fThreshVsW [view][hit.Cell()]->Fill(hit.W(), threshPart);
328  fShadowVsW [view][hit.Cell()]->Fill(hit.W(), shadowPart);
329  fRatioVsW [view][hit.Cell()]->Fill(hit.W(), ratio);
330 
331  fPEpercmVsMeVpercmTot[view] ->Fill(hit.TrueMeV() / hit.Path() , hit.PE() / hit.Path());
332  fPEpercmVsMeVpercm[view][hit.Cell()]->Fill(hit.TrueMeV() / hit.Path() , hit.PE() / hit.Path());
333  }
334  //LOG_VERBATIM("ThresholdAna") << "Finished running over PCHitlist\n";
336  evt.getByLabel("pclist", "pathQualXYThresh", hitsXYThresh);
337 
338  //LOG_VERBATIM("ThresholdAna") << "Filling hists\n";
339  for(unsigned int n = 0; n < hitsXYThresh->size(); ++n){
340  const caldp::PCHit& hit = (*hitsXYThresh)[n];
341  int plane = hit.Plane();
342  int view = GetView(plane,geom);
343  int brightness = fbhandle->getBrightnessIndex(plane, hit.Cell());
345  fNZerosVsBrightness[view]->Fill(brightness);
346  }
347  if(fFiberBrightnessMode && brightness != fFiberBrightness){
348  if(view > geo::kY && !fMuCFiberBrightnessMode && fFiberBrightness > 0) continue;
349  fNZerosVsW[view]->Fill(hit.W());
350  fNZerosVsCell[view]->Fill(hit.Cell());
351  }
352  }
353  //LOG_VERBATIM("ThresholdAna") << "Finished analyze\n";
354  }
355 
357 
358 } //end namespace
virtual void reconfigure(const fhicl::ParameterSet &pset)
float TrueMeV() const
Return True energy (MeV) value.
Definition: PCHit.h:44
TH2 * fRecoVsW[4][32 *12]
TH1 * ratio(TH1 *h1, TH1 *h2)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
int Cell() const
Return cell value.
Definition: PCHit.h:26
OStream cerr
Definition: OStream.cxx:7
TH2 * fPEpercmVsMeVpercm[4][32 *12]
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
TH2 * fRatioVsW[4][32 *12]
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
float Path() const
Return Path value.
Definition: PCHit.h:40
Definition: Run.h:31
virtual void beginRun(const art::Run &run) override
float W() const
Return W value.
Definition: PCHit.h:42
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
float PoissonLambda() const
Return number of simulated photons at readout before fluctuations.
Definition: PCHit.h:60
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
CDPStorage service.
void hits()
Definition: readHits.C:15
float PE() const
Return PE value.
Definition: PCHit.h:38
TH2 * fShadowVsW[4][32 *12]
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
Prototype Near Detector on the surface at FNAL.
TH2 * fTrueVsW[4][32 *12]
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
virtual void analyze(const art::Event &evt) override
Near Detector in the NuMI cavern.
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
int Plane() const
Return plane value.
Definition: PCHit.h:24
TH2 * fPEpercmVsW[4][32 *12]
TH2 * fThreshVsW[4][32 *12]
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
Make plots of threshold and shielding effects from PCHits files.
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
int GetView(int plane, art::ServiceHandle< geo::Geometry > geom)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
#define LOG_VERBATIM(category)
Encapsulate the geometry of one entire detector (near, far, ndos)
ThresholdAna(const fhicl::ParameterSet &pset)
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string