Calibrator.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file Calibrator.h
3 /// \brief Service that handles calibrations for a particular
4 /// validity context (eg: detector+timestamp).
5 /// \author Christopher Backhouse - bckhouse@caltech.edu
6 ////////////////////////////////////////////////////////////////////////
7 #ifndef CALIBRATOR_H
8 #define CALIBRATOR_H
9 
10 #include "TTimeStamp.h"
11 
12 // ART includes
17 #include "fhiclcpp/ParameterSet.h"
19 
20 
22 #include "Geometry/Geometry.h"
24 #include "RawData/RawDigit.h"
27 
28 #include "TF1.h"
29 
30 
31 namespace calib
32 {
33 
35  {
36  template<class T> using Atom = fhicl::Atom<T>;
38  using Name = fhicl::Name;
40  Name("ShapeTableFilenameFDMC"),
41  Comment("FDMC ADC chape fit table")
42  };
44  Name("ShapeTableFilenameFDData"),
45  Comment("FD Data ADC chape fit table")
46  };
48  Name("ShapeTableFilenameNDMC"),
49  Comment("NDMC ADC chape fit table")
50  };
52  Name("ShapeTableFilenameNDData"),
53  Comment("ND Data ADC chape fit table")
54  };
55  };
56 
58  {
59  template<class T> using Atom = fhicl::Atom<T>;
60  template<class T> using Sequence = fhicl::Sequence<T>;
61  template<class T> using Table = fhicl::Table<T>;
63  using Name = fhicl::Name;
64 
65 
66  // ADC fit file paths. Config could be cleaner
67  Table<calib::ShapeTableParams> DefaultShapeTable{
68  Name("default"),
69  Comment("Default ADC shape fit tables")
70  };
71  Table<calib::ShapeTableParams> NewFuncShapeTable{
72  Name("newfunc"),
73  Comment("New ADC shape fit tables")
74  };
75  Atom<std::string> CalibrationMode{
76  Name("CalibrationMode"),
77  Comment("Name of pset holding the desired ADC fits")
78  };
79 
80 
82  Name("Tag"),
83  Comment("Tag dictating which set of UPS CSV constants to calibrate with")
84  };
85  Atom<bool> UseAttenEpochs{
86  Name("UseAttenEpochs"),
87  Comment("Use run-boundary epoch strings from DB for retrieving attenuation CSVs")
88  };
89  Atom<std::string> AttenEpochTag{
90  Name("AttenEpochTag"),
91  Comment("Set DB tag for retrieving attenuation CSV epoch string")
92  };
93  Atom<bool> UseCSVsFromUPS{
94  Name("UseCSVsFromUPS"),
95  Comment("Use official calibrator-tagged CSVs; override this to use local constants")
96  };
97  Atom<std::string> AttenCSVPath{
98  Name("AttenCSVPath"),
99  Comment("Override path to use unofficial attenuation CSV files, when UseCSVsFromUPS=false ")
100  };
101  Atom<bool> UseAbsEpochs{
102  Name("UseAbsEpochs"),
103  Comment("Use run-boundary epoch strings from DB for retrieving absolute scale CSVs")
104  };
105  Atom<bool> ReadEpochsFromCSV{
106  Name("ReadEpochsFromCSV"),
107  Comment("Read the run-boundary epoch strings directly from the names of the CSVs")
108  };
109  Atom<std::string> AbsEpochTag{
110  Name("AbsEpochTag"),
111  Comment("Set DB tag for retrieving absolute CSV epoch string")
112  };
113  Atom<std::string> AbsConstsCSVPath{
114  Name("AbsConstsCSVPath"),
115  Comment("Override path to use unofficial absolute CSV files, when UseCSVsFromUPS=false ")
116  };
117  Atom<bool> UseTimingEpochs{
118  Name("UseTimingEpochs"),
119  Comment("Use run-boundary epoch strings from DB for retrieving timing CSVs")
120  };
121  Atom<std::string> TimingEpochTag{
122  Name("TimingEpochTag"),
123  Comment("Set DB tag for retrieving timing CSV epoch string")
124  };
125  Atom<std::string> TimingTag{
126  Name("TimingTag"),
127  Comment("Which timing calibration to use, tagged seperately from rest of Calibrator")
128  };
129  Atom<std::string> TimingConstsCSVPath{
130  Name("TimingConstsCSVPath"),
131  Comment("Override path to use unofficial timing CSV files, when UseCSVsFromUPS=false ")
132  };
133  Atom<bool> UseDrift{
134  Name("UseDrift"),
135  Comment("Use the drift calibration or not")
136  };
137  Atom<std::string> DriftConstsCSVPath{
138  Name("DriftConstsCSVPath"),
139  Comment("")
140  };
142  Name("DriftCSV"),
143  Comment("")
144  };
145  Atom<bool> UseDriftEpochs{
146  Name("UseDriftEpochs"),
147  Comment("")
148  };
149  Atom<std::string> DriftEpochTag{
150  Name("DriftEpochTag"),
151  Comment("")
152  };
153  Atom<bool> UseGainSetting4Abs{
154  Name("UseGainSetting4Abs"),
155  Comment("Reinitialize gain setting from RunHistoryService in postBeginRun")
156  };
157  Atom<bool> UseGainSetting4Atten{
158  Name("UseGainSetting4Atten"),
159  Comment("Reinitialize gain setting from RunHistoryService in postBeginRun")
160  };
162  Name("Gain"),
163  Comment("-1 to use Database, positive number to override")
164  };
165  Atom<bool> UseXFunction{
166  Name("UseXFunction"),
167  Comment("Apply W-dependent systematic shift on PECorr returned by Calibrator in X view")
168  };
169  Atom<bool> UseYFunction{
170  Name("UseYFunction"),
171  Comment("Apply W-dependent systematic shift on PECorr returned by Calibrator in Y view")
172  };
173  Atom<std::string> Xg100FunctionForm{
174  Name("Xg100FunctionForm"),
175  Comment("TF1 equation of W; shifts PECorr in X view, gain 100 mode")
176  };
177  Atom<std::string> XMuCg100FunctionForm{
178  Name("XMuCg100FunctionForm"),
179  Comment("TF1 equation of W; shifts PECorr in X view of ND Muon Catcher, Xg100FunctionForm used if not specified")
180  };
181  Atom<std::string> Yg100FunctionForm{
182  Name("Yg100FunctionForm"),
183  Comment("TF1 equation of W; shifts PECorr in Y view, gain 100 mode")
184  };
185  Atom<std::string> Xg140FunctionForm{
186  Name("Xg140FunctionForm"),
187  Comment("TF1 equation of W; shifts PECorr in X view, gain 140 mode (FD only)")
188  };
189  Atom<std::string> Yg140FunctionForm{
190  Name("Yg140FunctionForm"),
191  Comment("TF1 equation of W; shifts PECorr in Y view, gain 140 mode (FD only)")
192  };
193  Sequence<double> Xg100FunctionParams{
194  Name("Xg100FunctionParams"),
195  Comment("TF1 equation constants for Xg100FunctionForm")
196  };
197  Sequence<double> XMuCg100FunctionParams{
198  Name("XMuCg100FunctionParams"),
199  Comment("TF1 equation constants for XMuCg100FunctionForm")
200  };
201  Sequence<double> Yg100FunctionParams{
202  Name("Yg100FunctionParams"),
203  Comment("TF1 equation constants for Yg100FunctionForm")
204  };
205  Sequence<double> Xg140FunctionParams{
206  Name("Xg140FunctionParams"),
207  Comment("TF1 equation constants for Xg140FunctionForm")
208  };
209  Sequence<double> Yg140FunctionParams{
210  Name("Yg140FunctionParams"),
211  Comment("TF1 equation constants for Yg140FunctionForm")
212  };
213  Atom<bool> UseCellByCellCalibRandomOffset{
214  Name("UseCellByCellCalibRandomOffset"),
215  Comment("Smear attenuation scale cell-by-cell for systematic studies.")
216  };
217  Atom<bool> CellByCellCalibRandomOffsetHighResolution{
218  Name("CellByCellCalibRandomOffsetHighResolution"),
219  Comment("Smear attenuation scale within cell for systematic studiess.")
220  };
221  Atom<double> CellByCellCalibRandomOffset{
222  Name("CellByCellCalibRandomOffset"),
223  Comment("One-sigma value in TRandom3 smear for systematic studies")
224  };
225  Atom<int> CellByCellCalibRandomOffsetSeed{
226  Name("CellByCellCalibRandomOffsetSeed"),
227  Comment("TRandom3 seed for systematic studies, smearing either per-cell or cell-to-cell")
228  };
229  Atom<bool> MaskUncalibratedChannelsInMC{
230  Name("MaskUncalibratedChannelsInMC"),
231  Comment("Force channels uncalibrated in Data to also be uncalibrated in MC")
232  };
233  Atom<bool> MakeSinglePointCellHits{
234  Name("MakeSinglePointCellHits"),
235  Comment("When CalHit transforms RawDigits->CellHits, produce CellHits that look like v0 data, whatever the input")
236  };
237  Atom<bool> UseTimingOffsets{
238  Name("UseTimingOffsets"),
239  Comment("Turn off use of timing offsets for purpose of running the timing calibration")
240  };
241  Atom<bool> SimulateDrift{
242  Name("SimulateDrift"),
243  Comment("Whether to throw a linear time-based scale factor to absolute calibration")
244  };
245  Atom<int> DriftReference{
246  Name("DriftReference"),
247  Comment("Unix time specifying nominal time to simulate detector drift around")
248  };
249  Atom<double> DriftGradient{
250  Name("DriftGradient"),
251  Comment("Gradient of detector drift, with units of fraction/year")
252  };
253  };
254 }
255 
256 
257 namespace rb {
258  class CellHit;
259  class RecoHit;
260 }
261 
262 ///Calibration service
263 namespace calib {
264 
265  class AttenCurve;
266  class AbsCache;
267  class AttenCache;
268  class ADCShapeFitTable;
269  class DriftCache;
270  class TimingCache;
271 
272  class Calibrator {
273 
274  public:
276  Calibrator(const Parameters& params,
277  art::ActivityRegistry & reg);
278  ~Calibrator();
279 
280  void preEvent(art::Event const& evt);
281  void postBeginRun(art::Run const& run);
282  void postBeginSubRun(art::SubRun const& subrun);
283 
284  // Get the timestamp currently assumed by the Calibrator (unused for now)
285  TTimeStamp GetTimeStamp() const { return fCurrTimeStamp; }
286 
287  //---------------------------------------------
288  // RawDigit to CellHit methods
289  //---------------------------------------------
290 
291  // Main thing you'll want:
292  rb::CellHit MakeCellHit(const rawdata::RawDigit *rawdigit); // Do everything
293 
294  float GetPE (rb::CellHit const&);
295  float GetPE (const rawdata::RawDigit *rawdigit);
296  /// \param goodTime Did the timing fit go well?
297  /// \param maxadc Internal hook for use of \ref GetPE
298  float GetTNS(const rawdata::RawDigit *rawdigit,
299  bool & goodTime,
300  double* maxadc = 0);
301 
302  //---------------------------------------------
303  // CellHit to RecoHit methods
304  //---------------------------------------------
305  // (NOTE: w != dist_to_readout. Those are
306  // different arguments.)
307  //---------------------------------------------
308 
309  // Main thing you'll want:
310 
311  // Do everything
312  rb::RecoHit MakeRecoHit(rb::CellHit const& cellhit,
313  double w);
314  rb::RecoHit MakeRecoHit(rb::CellHit const& cellhit,
315  double* xyz);
316 
317  double GetT (rb::CellHit const& cellhit,
318  double const& dist_to_readout);
319  double GetPECorr (rb::CellHit const& cellhit,
320  double w);
321  double GetTimeRes(rb::CellHit const& cellhit);
322 
323  double GetAttenScale(rb::CellHit const& cellhit,
324  double w); ///< for PE->PECorr conversion
325  double GetAttenScale(const calib::AttenCurve* curve,
326  double const& minW,
327  double const& maxW,
328  double & w);
329 
330  /// Part of the PE->PECorr conversion
331  double GetDriftScale(rb::CellHit const& cellhit);
332 
333  // For PECorr->GeV conversion
334  double GetPECorrToGeVScale(rb::CellHit const& cellhit);
335 
336  /// For GeV->MIP conversion
337  double GetGeVToMIPScale(rb::CellHit const& cellhit);
338 
339  /// \brief For use by RecoBase classes
340  ///
341  /// Conversion from TotalGeV to CalorimetricEnergy. Takes account of dead
342  /// material and threshold effects in an averaged fashion. Makes no sense
343  /// to apply to individual hits. You probably don't want to use this, you
344  /// probably want to look at functions of your RecoBase object.
345  double GetGeVToCalorimetricScale() const;
346 
347  /// \brief Estimate of true shower energy, adapted from RecoJMShower
348  ///
349  /// Accounts for dead material and threshold effects that can't be done at
350  /// the individual cell level.
351  /// Pass in the TVectors and energy explicitly rather than the rb::Track
352  /// object to avoid a circular dependency with RecoBase
353  double GetShowerEnergy(TVector3 const& start,
354  TVector3 const& stop,
355  double const& totalGeV);
356 
357  /// Given plane, cell and whether the data is real,
358  /// get attenuation curve and Max w for the plane
359  bool GetAttenCurve(int const& plane,
360  int const& cell,
361  bool const& is_realdata,
362  const calib::AttenCurve*&,
363  double & minW,
364  double & maxW) const;
365 
366  /// Get the timing offset for a given plane, cell. Useful downstream to check calibration
367  double GetTimingOffset(unsigned int const& plane,
368  unsigned int const& cell,
369  bool const& isData);
370 
371  /// Method to determine the systematic uncertainty scale factor to apply
372  /// when those options are turned on.
373  /// This function should only be called by code wishing to check the
374  /// behavior of the systematic scaling, it is not useful outside of the
375  /// Calibrator service otherwise
376  double SystematicUncertaintyScale(geo::View_t const& view,
377  unsigned short const& plane,
378  double const& w);
379 
380 
382 
383  // Provide exact calibration configuration so
384  // there need not be confusion over what was done
385  std::string GetAttenCSVPath() {return (fParams.UseCSVsFromUPS()) ? getenv("CALIBCSVS_CSV_PATH") : fParams.AttenCSVPath();}
386  std::string GetAbsConstsCSVPath() {return (fParams.UseCSVsFromUPS()) ? getenv("CALIBCSVS_CSV_PATH") : fParams.AbsConstsCSVPath();}
387  std::string GetTimingConstsCSVPath(){return (fParams.UseCSVsFromUPS()) ? getenv("CALIBCSVS_CSV_PATH") : fParams.TimingConstsCSVPath();}
388  std::string GetDriftConstsCSVPath() {return (fParams.UseCSVsFromUPS()) ? getenv("CALIBCSVS_CSV_PATH") : fParams.DriftConstsCSVPath();}
390 
391  protected:
392  CalibratorParams fParams; // the fcl parameter table
393 
394  private:
395  TTimeStamp fCurrTimeStamp; ///< current time stamp
396 
398  /// Root file \ref ADCShapeFitTable will load from
403 
404  // Not initialized by default, must check and create if required
406  ADCShapeFitTable* fShapeTable5TB; ///< FEBv5 table for TestBeam (fShapeTable will be FEBv4)
407 
408  AttenCache* fAttenCacheData; ///< Data attenutation corrections
409  std::vector< AttenCache* > fAttenCacheMC; ///< MC attenutation corrections indexed by fiber brightness
413  bool fTimingCacheVldTimeSet; ///< Is fTimingCache set to the right time?
414  int fCurrentRun; ///< get the current run to make sure
415  ///< that we know if we have perfect MC
416 
417  /// Helper function for the two public MakeRecoHit methods
418  rb::RecoHit MakeRecoHit(rb::CellHit const& cellhit,
419  double & w,
420  double* xyz);
421 
422  void EnsureTimingCacheVldTime();
423 
424  /// Helper function for GetPE()
425  double GetAdcPerPE(const rawdata::RawDigit* dig);
426 
427  /// Get PE when we have already calculated the peakadc through a call to GetTNS
428  float GetPE (rb::CellHit const&, const double peakadc);
429  float GetPE (const rawdata::RawDigit *rawdigit, const double peakadc);
430 
431  // save some fcl parameters for configuring
432  // AttenCache vector in postBeginRun
434  double fGain;
435 
439 
440 
441  // Cache service handles for speed (we're in the inner loop of lots of
442  // algorithms)
447 
448  /// Number of calibrated hits returned
449  unsigned int fStatsNumCalibrated;
450  /// Number of uncalibrated hits returned
451  unsigned int fStatsNumUncalibrated;
452  /// Number of hits from uncalibrated diblocks using averaged calibrations
453  unsigned int fStatsNumAveraged;
454 
455  };
456 } // end namespace calib
458 #endif
459 ///////////////////////////////////////////////////////////////////////////
Atom< std::string > ShapeTableFilenameFDMC
Definition: Calibrator.h:39
Atom< std::string > ShapeTableFilenameFDData
Definition: Calibrator.h:43
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
TimingCache * fTimingCache
Definition: Calibrator.h:412
unsigned int fStatsNumUncalibrated
Number of uncalibrated hits returned.
Definition: Calibrator.h:451
TTimeStamp GetTimeStamp() const
Definition: Calibrator.h:285
art::ServiceHandle< geo::Geometry > fGeom
Definition: Calibrator.h:444
art::ServiceHandle< nova::dbi::RunHistoryService > fRH
Definition: Calibrator.h:446
#define DECLARE_ART_SERVICE(svc, scope)
Definition: ServiceMacros.h:91
DriftCache * fDriftCache
Definition: Calibrator.h:410
Definition: Run.h:31
art::ServiceHandle< cmap::CMap > fCMap
Definition: Calibrator.h:443
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
CDPStorage service.
std::string getenv(std::string const &name)
EFEBType
Enumeration for the FEB 4.2 and 5.1 which have different clocks.
Definition: CalibUtil.h:25
std::string fShapeTableFilenameFDMC
Root file ADCShapeFitTable will load from.
Definition: Calibrator.h:399
int evt
ADCShapeFitTable * fShapeTable
Definition: Calibrator.h:405
Atom< std::string > ShapeTableFilenameNDData
Definition: Calibrator.h:51
ADCShapeFitTable * fShapeTable5TB
FEBv5 table for TestBeam (fShapeTable will be FEBv4)
Definition: Calibrator.h:406
Atom< std::string > ShapeTableFilenameNDMC
Definition: Calibrator.h:47
std::string GetDriftConstsCSVPath()
Definition: Calibrator.h:388
Look up absolute attenuation constants.
Definition: AbsCache.h:16
std::string GetTimingConstsCSVPath()
Definition: Calibrator.h:387
Perform a "2 point" Hough transform on a collection of hits.
CalibratorParams fParams
Definition: Calibrator.h:392
bool fTimingCacheVldTimeSet
Is fTimingCache set to the right time?
Definition: Calibrator.h:413
std::string GetAbsConstsCSVPath()
Definition: Calibrator.h:386
Definition: run.py:1
std::string fShapeTableFilenameNDMC
Definition: Calibrator.h:401
TTimeStamp fCurrTimeStamp
current time stamp
Definition: Calibrator.h:395
AbsCache * fAbsCache
Definition: Calibrator.h:411
unsigned int fStatsNumAveraged
Number of hits from uncalibrated diblocks using averaged calibrations.
Definition: Calibrator.h:453
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
std::string fTag
Definition: Calibrator.h:433
fhicl::Comment Comment
Definition: Calibrator.h:37
std::string fShapeTableFilenameFDData
Definition: Calibrator.h:400
std::string fShapeTableFilenameNDData
Definition: Calibrator.h:402
std::vector< AttenCache * > fAttenCacheMC
MC attenutation corrections indexed by fiber brightness.
Definition: Calibrator.h:409
std::string GetAttenCSVPath()
Definition: Calibrator.h:385
art::ServiceHandle< photrans::FiberBrightness > fFibBrightness
Definition: Calibrator.h:445
AttenCache * fAttenCacheData
Data attenutation corrections.
Definition: Calibrator.h:408
unsigned int fStatsNumCalibrated
Number of calibrated hits returned.
Definition: Calibrator.h:449
CalibratorParams GetCalibratorParams()
Definition: Calibrator.h:389
Float_t w
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)