LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
ShowerCalibrationGaloreFromPID.cxx
Go to the documentation of this file.
1 
15 
16 
18 #include "cetlib_except/exception.h"
19 #include "cetlib/search_path.h"
20 
21 // ROOT libraries
22 #include "TGraphErrors.h"
23 #include "TFile.h"
24 #include "TDirectory.h"
25 
26 
27 // C/C++ standard libraries
28 #include <algorithm> // std::lower_bound(), std::sort(), std::min(), std::max()
29 #include <memory> // std::make_unique()
30 #include <utility> // std::pair<>
31 #include <tuple> // std::tie()
32 
33 
34 
35 //------------------------------------------------------------------------------
37  (recob::Shower const& shower, PDGID_t id /* = unknownID */) const
38 {
39  CalibrationInfo_t const& corr = selectCorrection(id);
40 
41  double const E = shower.Energy().at(shower.best_plane());
42 
43  return (float) corr.evalFactor(E);
44 } // lar::example::ShowerCalibrationGaloreFromPID::correctionFactor()
45 
46 
47 //------------------------------------------------------------------------------
50  (recob::Shower const& shower, PDGID_t id /* = unknownID */) const
51 {
52  CalibrationInfo_t const& corr = selectCorrection(id);
53 
54  double const E = shower.Energy().at(shower.best_plane());
55 
56  return { (float) corr.evalFactor(E), (float) corr.evalError(E) };
57 } // lar::example::ShowerCalibrationGaloreFromPID::correction()
58 
59 
60 //------------------------------------------------------------------------------
62  (std::string path)
63 {
64 
65  //
66  // open the input file
67  //
68  TDirectory* CalibDir = nullptr;
69  try {
70  CalibDir = OpenROOTdirectory(path);
71  }
72  catch (cet::exception const& e) {
73  // wrap the exception with context information
74  throw cet::exception
75  ("ShowerCalibrationGaloreFromPID", "readCalibration()", e)
76  << "Reading calibration from: '" << path << "'";
77  }
78  if (!CalibDir) { // this would be actually a bug
79  throw cet::exception
80  ("ShowerCalibrationGaloreFromPID", "readCalibration()")
81  << "Null directory reading calibration from: '" << path << "'";
82  }
83 
84  // make sure that when this is over we are done with the ROOT file
85  std::unique_ptr<TFile> InputFile(CalibDir->GetFile());
86 
87  //
88  // read each calibration object and associate it with its particle category
89  //
90 
91  Calibration_pi0 = readParticleCalibration(CalibDir, "Pi0", 111);
92 
93  Calibration_photon = readParticleCalibration(CalibDir, "Photon", 22);
94 
95  Calibration_electron
96  = readParticleCalibration(CalibDir, "Electron", { -11, 11 });
97 
98  Calibration_muon = readParticleCalibration(CalibDir, "Muon", { -13, 13 });
99 
100  Calibration_other = readParticleCalibration(CalibDir, "Default", unknownID);
101 
102  //
103  // release resources
104  //
105 
106  // TFile freed by its unique pointer
107 
108 } // lar::example::ShowerCalibrationGaloreFromPID::readCalibration()
109 
110 
111 //------------------------------------------------------------------------------
114 {
115  switch (id) {
116  case 111: // pi0
117  return Calibration_pi0;
118  case 22: // photon
119  return Calibration_photon;
120  case -11: // electron
121  case +11: // positron
122  return Calibration_electron;
123  case -13: // muon
124  case +13: // antimuon
125  return Calibration_muon;
126  case unknownID:
127  default:
128  return Calibration_other;
129  } // switch
130 } // lar::example::ShowerCalibrationGaloreFromPID::selectCorrection()
131 
132 
133 //------------------------------------------------------------------------------
136  (TDirectory* SourceDir, std::string GraphName) const
137 {
138  CalibrationInfo_t info;
139 
140  // apply list is left empty
141 
142  //
143  // retrieve the object
144  //
145  auto graph = details::readROOTobject<TGraphErrors>(SourceDir, GraphName);
146 
147  verifyOrder(graph.get());
148 
149  size_t const N = (size_t) graph->GetN();
150  if (N == 0) {
151  throw cet::exception("ShowerCalibrationGaloreFromPID")
152  << "No point in graph " << SourceDir->GetPath() << "/" << GraphName
153  << "\n";
154  }
155 
156  // include the "error" on x in the full range
157  info.minE = graph->GetX()[0];
158  info.maxE = graph->GetX()[N - 1];
159 
160  // a spline; if there are at least 5 points, use AKIMA spline, that is
161  // "stable" for outliers (reducing over/undershoot)
162  // set to zero
163  info.factor = createInterpolator(N, graph->GetX(), graph->GetY());
164 
165  // compute the error in the same way; kind of an approximation here
166  info.error = createInterpolator(N, graph->GetX(), graph->GetEY());
167 
168  return info;
169 } // lar::example::ShowerCalibrationGaloreFromPID::readParticleCalibration()
170 
171 
172 //------------------------------------------------------------------------------
174  (TGraph const* graph)
175 {
176 
177  if (!graph) {
178  throw cet::exception("ShowerCalibrationGaloreFromPID")
179  << "VerifyOrder(): invalid graph specified\n";
180  }
181 
182  size_t const N = graph->GetN();
183  if (N < 2) return;
184 
185  Double_t const* x = graph->GetX();
186 
187  for (size_t i = 1; i < N; ++i) {
188  if (x[i-1] > x[i]) {
189  throw cet::exception("ShowerCalibrationGaloreFromPID")
190  << "VerifyOrder(): points in graph '" << graph->GetName()
191  << "' are not sorted in abscissa\n";
192  }
193  } // while
194 } // lar::example::ShowerCalibrationGaloreFromPID::verifyOrder()
195 
196 
197 //------------------------------------------------------------------------------
200  (TDirectory* SourceDir, std::string GraphName, PDGID_t id) const
201 {
202  CalibrationInfo_t info = readParticleCalibration(SourceDir, GraphName);
203  info.applyTo(id);
204  return info;
205 } // lar::example::ShowerCalibrationGaloreFromPID::readParticleCalibration(ID)
206 
207 
210  TDirectory* SourceDir, std::string GraphName,
211  std::initializer_list<PDGID_t> ids
212 ) const {
213  CalibrationInfo_t info = readParticleCalibration(SourceDir, GraphName);
214  info.applyTo(ids);
215  return info;
216 } // lar::example::ShowerCalibrationGaloreFromPID::readParticleCalibration(IDs)
217 
218 
219 //------------------------------------------------------------------------------
221  (std::string path)
222 {
223  //
224  // split the data file path
225  //
226  std::string filePath, ROOTdirPath;
227  std::tie(filePath, ROOTdirPath) = splitROOTpath(path);
228 
229  //
230  // find the ROOT file in the search path
231  //
232  std::string fullFilePath;
233  cet::search_path sp("FW_SEARCH_PATH");
234  // if we can't find the file in FW_SEARCH_PATH, we try in current directory
235  if (!sp.find_file(filePath, fullFilePath)) fullFilePath = filePath;
236 
237  //
238  // open the ROOT file (created new)
239  //
240  auto inputFile = std::make_unique<TFile>(fullFilePath.c_str(), "READ");
241  if (!(inputFile->IsOpen())) {
242  throw cet::exception("ShowerCalibrationGaloreFromPID")
243  << "ShowerCalibrationGaloreFromPID::OpenROOTdirectory() can't read '"
244  << fullFilePath << "' (from '" << filePath << "' specification)\n";
245  }
246 
247  //
248  // get the ROOT directory
249  //
250  TDirectory* dir = ROOTdirPath.empty()?
251  inputFile.get(): inputFile->GetDirectory(ROOTdirPath.c_str());
252  if (!dir) {
253  throw cet::exception("ShowerCalibrationGaloreFromPID")
254  << "ShowerCalibrationGaloreFromPID::OpenROOTdirectory() can't find '"
255  << ROOTdirPath << "' in ROOT file '" << inputFile->GetPath() << "'\n";
256  }
257 
258  //
259  // return the directory
260  //
261  inputFile.release(); // do not delete the file any more
262  return dir;
263 } // lar::example::ShowerCalibrationGaloreFromPID::OpenROOTdirectory()
264 
265 
266 //------------------------------------------------------------------------------
267 // lar::example::ShowerCalibrationGaloreFromPID::CalibrationInfo_t
268 //
271  (PDGID_t id)
272 {
273  auto it = std::lower_bound(appliesTo.begin(), appliesTo.end(), id);
274  if ((it == appliesTo.end()) || (*it != id))
275  appliesTo.insert(it, id);
276  return *this;
277 } // ShowerCalibrationGaloreFromPID::CalibrationInfo_t::applyTo(ID)
278 
279 
282  (std::initializer_list<PDGID_t> ids)
283 {
284  appliesTo.insert(appliesTo.begin(), ids.begin(), ids.end());
285  std::sort(appliesTo.begin(), appliesTo.end());
286  return *this;
287 } // ShowerCalibrationGaloreFromPID::CalibrationInfo_t::applyTo(IDs)
288 
289 
290 //--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ---
291 double
293  (double E) const
294 {
295  double const boundE = std::min(maxE, std::max(minE, E));
296  return factor->Eval(boundE);
297 }
298 
299 double
301  (double E) const
302 {
303  double const boundE = std::min(maxE, std::max(minE, E));
304  return error->Eval(boundE);
305 }
306 
307 
308 //------------------------------------------------------------------------------
309 std::unique_ptr<ROOT::Math::Interpolator>
311  (unsigned int N, double const* x, double const* y)
312 {
313 
314  // decide the type of interpolation based on the available number of points
315  ROOT::Math::Interpolation::Type type;
316  if (N >= 5) type = ROOT::Math::Interpolation::kAKIMA;
317  else if (N >= 3) type = ROOT::Math::Interpolation::kCSPLINE;
318  else type = ROOT::Math::Interpolation::kLINEAR;
319 
320  auto interp
321  = std::make_unique<ROOT::Math::Interpolator>(std::max(N, 2U), type);
322  if (N > 1) interp->SetData(N, x, y);
323  else { // we need to make up the second point
324  double const x_p[2] = { *x, *x + 1. };
325  double const y_p[2] = { *y, *y };
326  interp->SetData(2, x_p, y_p);
327  }
328  return interp;
329 } // lar::example::ShowerCalibrationGaloreFromPID::createInterpolator()
330 
331 
332 //------------------------------------------------------------------------------
333 
334 std::pair<std::string, std::string>
335 lar::example::splitROOTpath(std::string path) {
336  const std::string suffix = ".root";
337 
338  // find the ROOT file name
339  std::string::size_type iSuffix = std::string::npos;
340  do {
341  iSuffix = path.rfind(suffix, iSuffix);
342  if (iSuffix == std::string::npos) return {}; // failure: no suffix
343 
344  // if it's not "suffix:" or "suffix/" or at end of string, it's not it
345  auto iAfter = iSuffix + suffix.length();
346  if ((iAfter < path.length())
347  && (path[iAfter] != '/')
348  && (path[iAfter] != ':')
349  )
350  {
351  // oops... this suffix is invalid; keep searching
352  if (iSuffix == 0) return {}; // failure: no suffix
353  --iSuffix;
354  continue;
355  }
356 
357  // we found a proper suffix
358  std::pair<std::string, std::string> result;
359 
360  result.first = path.substr(0U, iAfter); // file path
361  if (iAfter < path.length())
362  result.second = path.substr(iAfter + 1, path.length()); // ROOT dir
363  return result;
364 
365  } while (true);
366 } // lar::example::splitROOTpath()
367 
368 
369 //------------------------------------------------------------------------------
CalibrationInfo_t Calibration_pi0
neutral pion calibration
Float_t x
Definition: compare.C:6
int best_plane() const
Definition: Shower.h:200
virtual Correction_t correction(recob::Shower const &, PDGID_t=unknownID) const override
Returns the correction for a given reconstructed shower.
CalibrationInfo_t Calibration_electron
electron/positron calibration
Internal structure containing the calibration information.
Float_t E
Definition: plot.C:23
std::unique_ptr< ROOT::Math::Interpolator > error
parametrisation of the correction uncertainty
Float_t y
Definition: compare.C:6
const std::vector< double > & Energy() const
Definition: Shower.h:195
CalibrationInfo_t const & selectCorrection(PDGID_t id) const
Returns the correct CalibrationInfo_t for specified id.
static std::unique_ptr< ROOT::Math::Interpolator > createInterpolator(unsigned int N, double const *x, double const *y)
Creates a ROOT interpolator from a set of N points.
A correction factor with global uncertainty.
CalibrationInfo_t Calibration_photon
photon calibration
Int_t max
Definition: plot.C:27
CalibrationInfo_t Calibration_other
default calibration
CalibrationInfo_t readParticleCalibration(TDirectory *SourceDir, std::string GraphName) const
Reads and returns calibration information from the specified graph.
CalibrationInfo_t Calibration_muon
muon/antimuon calibration
std::unique_ptr< ROOT::Math::Interpolator > factor
parametrisation of the correction factor
virtual float correctionFactor(recob::Shower const &, PDGID_t=unknownID) const override
Returns a correction factor for a given reconstructed shower.
static constexpr PDGID_t unknownID
A mnemonic constant for unknown particle ID.
int PDGID_t
A type representing a particle ID in Particle Data Group convention.
TDirectory * dir
Definition: macro.C:5
Int_t min
Definition: plot.C:26
double minE
Definition: plot_hist.C:7
Shower energy calibration according to particle type.
double maxE
Definition: plot_hist.C:8
std::pair< std::string, std::string > splitROOTpath(std::string path)
Splits path into ROOT file name and directory path.
Float_t e
Definition: plot.C:34
void readCalibration(std::string path)
Reads the calibration information from the specified file.
static TDirectory * OpenROOTdirectory(std::string path)
Opens the specified ROOT directory, as in path/to/file.root:dir/dir.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33