LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
OpDigiProperties.cc
Go to the documentation of this file.
1 //
3 // \file OpDigiProperties_service.cc
4 //
6 
7 // LArSoft includes
9 
10 // Framework includes
12 
13 // CLHEP includes
14 #include "CLHEP/Random/RandGauss.h"
15 
16 // ROOT includes
17 #include "TF1.h"
18 
19 // C++ includes
20 #include <fstream>
21 
22 namespace opdet {
23 
24  //--------------------------------------------------------------------
26  {
27  fSampleFreq = p.get<double>("SampleFreq");
28  fTimeBegin = p.get<double>("TimeBegin");
29  fTimeEnd = p.get<double>("TimeEnd");
30  fWaveformFile = p.get<std::string>("WaveformFile");
31  fPERescale = p.get<double>("PERescale");
32 
33  // PMT properties
34  fQE = p.get<double>("QE");
35  fDarkRate = p.get<double>("DarkRate");
36  fSaturationScale = p.get<optdata::ADC_Count_t>("SaturationScale");
37 
38  // Shaper properties
39  fUseEmpiricalGain = p.get<bool>("UseEmpiricalGain");
40  fHighGainFile = p.get<std::string>("HighGainFile");
41  fLowGainFile = p.get<std::string>("LowGainFile");
42  fGainSpreadFile = p.get<std::string>("GainSpreadFile");
43 
44  fHighGainMean = p.get<double>("HighGainMean");
45  fLowGainMean = p.get<double>("LowGainMean");
46  fGainSpread = p.get<double>("GainSpread");
47  fGainSpread_PMT2PMT = p.get<double>("GainSpread_PMT2PMT");
48 
49  // Digitization ped fluc
50  fPedFlucRate = p.get<double>("PedFlucRate");
51  fPedFlucAmp = p.get<optdata::ADC_Count_t>("PedFlucAmp");
52  fADCBaseline = p.get<optdata::ADC_Count_t>("ADCBaseline");
53  fADCBaseSpread = p.get<double>("ADCBaseSpread");
54 
55  // WF related stuff
56  fUseEmpiricalShape = p.get<bool>("UseEmpiricalShape");
57  fWFLength = p.get<double>("WFLength");
58 
59  fWaveformFile = p.get<std::string>("WaveformFile");
60  fChargeNormalized = p.get<bool>("WaveformChargeNormalized", false);
61 
62  // Option 2: WF from analytical function
63  fWFPowerFactor = p.get<double>("WFPowerFactor");
64  fWFTimeConstant = p.get<double>("WFTimeConstant");
65  fVoltageAmpForSPE = p.get<double>("VoltageAmpForSPE");
66 
67  // Generate the SPE waveform (i.e. fWaveform)
69 
70  // Fill gain array
71  FillGainArray();
72 
73  // Fill baseline mean
75 
76  // Report
77  std::string msg(Form(
78  "%-10s ... %-10s ... %-10s ... %-10s\n", "Ch. Number", "Pedestal", "High Gain", "Low Gain"));
79  for (unsigned int i = 0; i < fGeometry->NOpChannels(); ++i) {
80  msg += Form("%-10d ... %-10d ... %-10g ... %-10g\n",
81  i,
82  fPedMeanArray[i],
83  fHighGainArray[i],
84  fLowGainArray[i]);
85  }
86  mf::LogInfo(__FUNCTION__) << msg.c_str();
87  }
88 
89  //--------------------------------------------------------------------
91  {
92  double SPEArea = 0;
93  for (size_t i = 0; i != fWaveform.size(); ++i)
94  SPEArea += fWaveform.at(i);
95  return SPEArea;
96  }
97 
98  //--------------------------------------------------------------------
100  {
101  double AmpSoFar = 0;
102  for (size_t i = 0; i != fWaveform.size(); ++i)
103  if (fWaveform.at(i) > AmpSoFar) AmpSoFar = fWaveform.at(i);
104  return AmpSoFar;
105  }
106 
107  //--------------------------------------------------------------------
109  {
110  double Cumulative = 0, SPEArea = 0;
111  for (size_t i = 0; i != fWaveform.size(); ++i) {
112  Cumulative += fWaveform.at(i);
113  SPEArea += Cumulative;
114  }
115  return SPEArea;
116  }
117 
118  //--------------------------------------------------------------------
120  {
121  double AmpSoFar = 0, Cumulative = 0;
122  for (size_t i = 0; i != fWaveform.size(); ++i) {
123  Cumulative += fWaveform.at(i);
124  if (Cumulative > AmpSoFar) AmpSoFar = Cumulative;
125  }
126  return AmpSoFar;
127  }
128 
129  //--------------------------------------------------------------------
130 
132  {
133  return fLowGainArray[ch];
134  }
135 
136  //--------------------------------------------------------------------
138  {
139  return fHighGainArray[ch];
140  }
141  //--------------------------------------------------------------------
143  {
144  return CLHEP::RandGauss::shoot(fLowGainArray[ch], fGainSpreadArray[ch] * fLowGainArray[ch]);
145  }
146  //--------------------------------------------------------------------
148  {
149  return CLHEP::RandGauss::shoot(fHighGainArray[ch], fGainSpreadArray[ch] * fHighGainArray[ch]);
150  }
151  //--------------------------------------------------------------------
153  {
154  if (time_ns / 1.e3 > (fTimeEnd - fTimeBegin))
155  return std::numeric_limits<optdata::TimeSlice_t>::max();
156 
157  else
158  return optdata::TimeSlice_t((time_ns / 1.e3 - fTimeBegin) * fSampleFreq);
159  }
160 
161  //
162  // As far as Kazu is concerned, this function is deprecated.
163  // Any comment? --Kazu 08/05/2013
164  //
165  std::vector<double> OpDigiProperties::WaveformInit(std::string fWaveformFile)
166  {
167  // Read in waveform vector from text file
168  std::ifstream WaveformFile(fWaveformFile.c_str());
169  std::string line;
170 
171  mf::LogInfo("OpDigiProperties")
172  << "OpDigiProperties opening OpDet waveform at " << fWaveformFile.c_str();
173 
174  // Read in each line and place into waveform vector
175  std::vector<double> PEWaveform;
176  if (WaveformFile.is_open()) {
177  while (WaveformFile.good()) {
178  getline(WaveformFile, line);
179  PEWaveform.push_back(fPERescale * strtod(line.c_str(), NULL));
180  }
181  }
182  else
183  throw cet::exception("OpDigiProperties") << "No Waveform File: Unable to open file\n";
184 
185  WaveformFile.close();
186  return (PEWaveform);
187  }
188 
189  // Fill the array of pedestal mean
191  {
192  for (unsigned int i = 0; i < fGeometry->NOpChannels(); ++i)
194  CLHEP::RandGauss::shoot((double)fADCBaseline, fADCBaseSpread) + 0.5));
195  }
196 
197  // Fill arrays (std::vector<double>) for PMT gain mean & spread information.
199  {
200  if (fUseEmpiricalGain) {
201  // Fill fron user's text files.
202  mf::LogWarning("OpDigiProperties") << "Using empirical table of gain for each PMT...";
203  std::string FullPath;
204  cet::search_path sp("FW_SEARCH_PATH");
205 
206  if (!sp.find_file(fHighGainFile, FullPath))
207  throw cet::exception("OpDigiProperties")
208  << "Unable to find high gain spread file in " << sp.to_string() << "\n";
209 
210  mf::LogWarning("OpDigiProperties")
211  << "OpDigiProperties opening high gain spread file at " << FullPath.c_str();
212  std::ifstream HighGainFile(FullPath.c_str());
213  if (HighGainFile.is_open()) {
214  std::string line;
215  while (HighGainFile.good()) {
216  getline(HighGainFile, line);
217  fHighGainArray.push_back(strtod(line.c_str(), NULL));
218  }
219  }
220  else
221  throw cet::exception("OpDigiProperties") << "Unable to open file!\n";
222 
223  FullPath = "";
224  if (!sp.find_file(fLowGainFile, FullPath))
225  throw cet::exception("OpDigiProperties")
226  << "Unable to find low gain spread file in " << sp.to_string() << "\n";
227 
228  mf::LogWarning("OpDigiProperties")
229  << "OpDigiProperties opening low gain spread file at " << FullPath.c_str();
230  std::ifstream LowGainFile(FullPath.c_str());
231  if (LowGainFile.is_open()) {
232  std::string line;
233  while (LowGainFile.good()) {
234  getline(LowGainFile, line);
235  fLowGainArray.push_back(strtod(line.c_str(), NULL));
236  }
237  }
238  else
239  throw cet::exception("OpDigiProperties") << "Unable to open file!\n";
240 
241  FullPath = "";
242  if (!sp.find_file(fGainSpreadFile, FullPath))
243  throw cet::exception("OpDigiProperties")
244  << "Unable to find low gain spread file in " << sp.to_string() << "\n";
245 
246  mf::LogWarning("OpDigiProperties")
247  << "OpDigiProperties opening low gain spread file at " << FullPath.c_str();
248  std::ifstream GainSpreadFile(FullPath.c_str());
249  if (GainSpreadFile.is_open()) {
250  std::string line;
251  while (GainSpreadFile.good()) {
252  getline(GainSpreadFile, line);
253  fGainSpreadArray.push_back(strtod(line.c_str(), NULL));
254  }
255  }
256  else
257  throw cet::exception("OpDigiProperties") << "Unable to open file!\n";
258  }
259  else {
260  // Generate a set of gake gain mean & spread.
261  std::string txt;
262  txt += " Generating gain for each pmt.\n";
263  txt += Form(" High gain mean: %g ADC/p.e.\n", fHighGainMean);
264  txt += Form(" Low gain mean: %g ADC/p.e.\n", fLowGainMean);
265  txt += Form(" PMT-to-PMT gain spread : %g \n", fGainSpread_PMT2PMT);
266  txt += Form(" Intrinsic gain spread : %g \n", fGainSpread);
267  mf::LogWarning("OpDigiProperties") << txt.c_str();
268  for (unsigned int i = 0; i < fGeometry->NOpChannels(); ++i) {
269  fLowGainArray.push_back(
270  CLHEP::RandGauss::shoot(fLowGainMean, fLowGainMean * fGainSpread_PMT2PMT));
271  fHighGainArray.push_back(
272  CLHEP::RandGauss::shoot(fHighGainMean, fHighGainMean * fGainSpread_PMT2PMT));
273  fGainSpreadArray.push_back(fGainSpread);
274  }
275  }
276 
277  //
278  // Note:
279  // Check for # entries. These exceptions ensures we have enough # of elements.
280  // When a user access the elements by a channel number using LowGainMean(Channel_t ch) etc.,
281  // those functions do not check if the given channel number is valid or not to keep a high
282  // performance of the code. If there's an invalid memory access error in run-time, then
283  // it must mean the user provided an invalid channel number and not due to insufficient
284  // vector elements filled in this function.
285  //
286  if (fLowGainArray.size() < fGeometry->NOpChannels())
287  throw cet::exception("OpDigiProperties") << "Low gain missing for some channels!\n";
288  if (fHighGainArray.size() < fGeometry->NOpChannels())
289  throw cet::exception("OpDigiProperties") << "High gain missing for some channels!\n";
290  if (fGainSpreadArray.size() < fGeometry->NOpChannels())
291  throw cet::exception("OpDigiProperties") << "Gain spread missing for some channels!\n";
292  }
293 
295  {
296  std::string FullPath;
297 
298  if (fUseEmpiricalShape) {
299  cet::search_path sp("FW_SEARCH_PATH");
300  if (!sp.find_file(fWaveformFile, FullPath))
301 
302  throw cet::exception("OpDigiProperties")
303  << "Unable to find PMT waveform file in " << sp.to_string() << "\n";
304 
305  fWaveform = GenEmpiricalWF(FullPath);
306  }
307  else
309  }
310 
311  std::vector<double> OpDigiProperties::GenEmpiricalWF(std::string fWaveformFile)
312  {
313  // Read in waveform vector from text file
314  std::ifstream WaveformFile(fWaveformFile.c_str());
315  std::string line;
316 
317  mf::LogWarning("OpDigiProperties")
318  << "OpDigiProperties opening OpDet waveform at " << fWaveformFile.c_str();
319 
320  // Read in each line and place into waveform vector
321  std::vector<double> PEWaveform;
322  if (WaveformFile.is_open()) {
323  double MaxAmp = 0;
324  double Charge = 0;
325  int NSample = 0;
326  while (WaveformFile.good() && NSample < int(fWFLength * fSampleFreq)) {
327  getline(WaveformFile, line);
328  double Amp = strtod(line.c_str(), NULL);
329  PEWaveform.push_back(Amp);
330  if (Amp > MaxAmp) MaxAmp = Amp;
331  NSample++;
332  Charge += Amp;
333  }
334  // rescale
335  if (MaxAmp <= 0)
336  throw cet::exception("OpDigiProperties_module") << "Waveform amplitude <=0!\n";
337  if (!fChargeNormalized)
338  for (unsigned short i = 0; i < PEWaveform.size(); i++) {
339  PEWaveform[i] = PEWaveform[i] / MaxAmp;
340  }
341  else
342  for (unsigned short i = 0; i < PEWaveform.size(); i++) {
343  PEWaveform[i] = PEWaveform[i] / Charge;
344  }
345  }
346  else
347  throw cet::exception("No Waveform File") << "Unable to open file\n";
348 
349  WaveformFile.close();
350  return (PEWaveform);
351  }
352 
353  std::vector<double> OpDigiProperties::GenAnalyticalWF()
354  {
355  mf::LogWarning("OpDigiProperties")
356  << " OpDigiProperties using analytical function for WF generation.";
357  //
358  // Generate waveform from analytical form
359  //
360  if (!fAnalyticalSPE) { // Create analytical function if not yet created
361  fAnalyticalSPE = new TF1("_analyticalSPE",
362  "10^(22)*x^[1]*[0]*exp(-x/[2])/TMath::Factorial([1])",
363  0,
364  fWFLength * 2); // x is time in micro-seconds
365  fAnalyticalSPE->SetParameter(0, fVoltageAmpForSPE); // voltage amplitude for s.p.e.
366  fAnalyticalSPE->SetParameter(1, fWFPowerFactor); // power factor (no unit)
367  fAnalyticalSPE->SetParameter(2, fWFTimeConstant); // time constant in us
368  }
369  //
370  // Define a waveform vector
371  //
372  // Size of WF = time width [us] * frequency [MHz]
373  std::vector<double> PEWaveform(int(fWFLength * fSampleFreq), 0.0);
374  double SamplingDuration = 1. / fSampleFreq; // in micro seconds
375  double MaxAmp = 0;
376  double Charge = 0;
377  for (unsigned short i = 0; i < PEWaveform.size(); ++i) {
378  double Value = fAnalyticalSPE->Integral(i * SamplingDuration, (i + 1) * SamplingDuration) /
379  SamplingDuration;
380  PEWaveform[i] = Value;
381  if (PEWaveform[i] > MaxAmp) MaxAmp = PEWaveform[i];
382  Charge += Value;
383  }
384 
385  // rescale
386  if (MaxAmp <= 0) throw cet::exception("OpDigiProperties_module") << "Waveform amplitude <=0!\n";
387 
388  if (!fChargeNormalized) {
389  for (unsigned short i = 0; i < PEWaveform.size(); i++) {
390  PEWaveform[i] = PEWaveform[i] / MaxAmp;
391  if (PEWaveform[i] < 1.e-4) PEWaveform[i] = 0;
392  }
393  }
394  else
395  for (unsigned short i = 0; i < PEWaveform.size(); i++) {
396  PEWaveform[i] = PEWaveform[i] / Charge;
397  }
398 
399  return PEWaveform;
400  }
401 
402 } // namespace
std::vector< double > fHighGainArray
art::ServiceHandle< geo::Geometry const > fGeometry
OpDigiProperties(fhicl::ParameterSet const &pset)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< double > GenAnalyticalWF()
double GetSPEAmplitude()
Utility function ... To be verified (Kazu 08/05/13)
std::vector< double > fLowGainArray
unsigned int NOpChannels() const
Number of electronics channels for all the optical detectors.
double LowGainMean() const noexcept
Returns set mean gain value for LOW gain.
double GetSPECumulativeArea()
Utility function ... To be verified (Kazu 08/05/13)
uint16_t ADC_Count_t
Definition: OpticalTypes.h:16
std::vector< double > WaveformInit(std::string WaveformFile)
double GetSPECumulativeAmplitude()
Utility function ... To be verified (Kazu 08/05/13)
std::vector< double > GenEmpiricalWF(std::string WaveformFile)
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::vector< double > fWaveform
optdata::ADC_Count_t fPedFlucAmp
unsigned int TimeSlice_t
Definition: OpticalTypes.h:20
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double LowGain(optdata::Channel_t ch) const
Generate & return LOW gain value for an input channel using mean & spread for this channel...
double HighGain(optdata::Channel_t ch) const
Generate & return HIGH gain value for an input channel using mean & spread for this channel...
std::vector< optdata::ADC_Count_t > fPedMeanArray
optdata::TimeSlice_t GetTimeSlice(double time_ns)
std::vector< double > fGainSpreadArray
Float_t e
Definition: plot.C:35
optdata::ADC_Count_t fADCBaseline
double GetSPEArea()
Utility function ... To be verified (Kazu 08/05/13)
unsigned int Channel_t
Definition: OpticalTypes.h:19
optdata::ADC_Count_t fSaturationScale
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double HighGainMean() const noexcept
Returns set mean gain value for HIGH gain.