test_ExponentialChannelNoiseService.cxx
Go to the documentation of this file.
1 // test_ExponentialChannelNoiseService.cxx
2 
3 // David Adams
4 // February 2015
5 //
6 // Test ExponentialChannelNoiseService.
7 
8 #include "../ExponentialChannelNoiseService.h"
9 
10 #include <fstream>
11 #include <iomanip>
12 #include <iostream>
13 #include <map>
14 #include <sstream>
15 #include <string>
16 
18 #include "art_root_io/TFileService.h"
20 #include "fhiclcpp/ParameterSet.h"
25 #include "lardataobj/RawData/raw.h"
26 
27 using std::string;
28 using std::cout;
29 using std::endl;
30 using std::istringstream;
31 using std::ostringstream;
32 using std::ofstream;
33 using std::setw;
34 using art::ServiceHandle;
35 using util::LArFFT;
36 
37 typedef std::vector<short> AdcVector;
38 
39 #undef NDEBUG
40 #include <cassert>
41 
42 // service configuration strings (definitions are at the bottom)
45 
46 namespace {
47 
48  std::string map_to_string(std::map<std::string, std::string> const& configs)
49  {
50  std::string full_config;
51  for (auto const& [service, config] : configs) {
52  full_config.append(service)
53  .append(": {")
54  .append(config)
55  .append("}\n");
56  }
57  return full_config;
58  }
59 
60  void load_services(int fftsize =0) {
61  const string myname = "load_services: ";
62  const string line = "-----------------------------";
63 
64  cout << myname << line << endl;
65  cout << myname << "Fetch art service helper." << endl;
66 
67  std::map<std::string, std::string> service_configs;
68  service_configs.emplace("TFileService", "fileName: 'test.root'");
69  string gname = "dune35t4apa_v6";
70  service_configs.emplace("Geometry",
71  "DisableWiresInG4: true GDML: \"dune35t4apa_v6.gdml\" Name: \"" + gname +
72  "\" ROOT: \"" + gname + "\""
73  " SortingParameters: { DetectorVersion: \"" + gname + "\""
74  " ChannelsPerOpDet: 12" +
75  "} SurfaceY: 0");
76  service_configs.emplace("ExptGeoHelperInterface", "service_provider: \"DUNEGeometryHelper\"");
77  service_configs.emplace("LArPropertiesService", LArPropertiesServiceConfigurationString);
78  service_configs.emplace("DetectorPropertiesService", DetectorPropertiesServiceConfigurationString);
79  service_configs.emplace("DetectorClocksService", "ClockSpeedExternal: 3.125e1 ClockSpeedOptical: 128 ClockSpeedTPC: 2 ClockSpeedTrigger: 16 DefaultBeamTime: 0 DefaultTrigTime: 0 FramePeriod: 1600 G4RefTime: 0 InheritClockConfig: false TrigModuleName: \"\" TriggerOffsetTPC: 0 service_provider: \"DetectorClocksServiceStandard\" ");
80  ostringstream sscfg;
81  sscfg << "FFTOption: \"\" FFTSize: " << fftsize << " FitBins: 5";
82 
83  service_configs.emplace("LArFFT", sscfg.str());
84  ArtServiceHelper::load_services(map_to_string(service_configs));
85 
86  // Work around defect in larsoft: https://cdcvs.fnal.gov/redmine/issues/10618.
87  TH1::AddDirectory(kFALSE);
88  }
89 }
90 
91 int test_ExponentialChannelNoiseService(unsigned int ntick, unsigned int maxchan =100) {
92  const string myname = "test_ExponentialChannelNoiseService: ";
93 #ifdef NDEBUG
94  cout << myname << "NDEBUG must be off." << endl;
95  abort();
96 #endif
97  const string line = "-----------------------------";
98 
99  load_services();
100 
101  cout << myname << line << endl;
102  cout << myname << "Create noise service." << endl;
103  fhicl::ParameterSet pset;
104  pset.put("NoiseNormZ", 3.16);
105  pset.put("NoiseNormU", 3.16);
106  pset.put("NoiseNormV", 3.16);
107  pset.put("NoiseWidthZ", 2000.0);
108  pset.put("NoiseWidthU", 2000.0);
109  pset.put("NoiseWidthV", 2000.0);
110  pset.put("LowCutoffZ", 7.5);
111  pset.put("LowCutoffU", 7.5);
112  pset.put("LowCutoffV", 7.5);
113  pset.put("WhiteNoiseZ", 0.0);
114  pset.put("WhiteNoiseU", 0.0);
115  pset.put("WhiteNoiseV", 0.0);
116  pset.put("NoiseArrayPoints", 1000);
117  pset.put("OldNoiseIndex", true);
118  pset.put("RandomSeed", 54321);
119  ExponentialChannelNoiseService noisvc(pset);
120 
121  ServiceHandle<LArFFT> hfftsvc;
123 
124  AdcSignal ped = 1000;
125  cout << myname << line << endl;
126  cout << myname << "Create signal with ped=" << ped << " and " << ntick << " ticks." << endl;
127  cout << myname << "FFT size: " << hfftsvc->FFTSize() << endl;
128  cout << myname << "Detector: " << hgeosvc->DetectorName() << endl;
129  AdcSignalVector sigsin(ntick, ped);
130  for ( unsigned int isig=0; isig<20; ++isig ) {
131  }
132  double sum = 0.0;
133  double sumsq = 0.0;
134  double sumdsq = 0.0;
135  double sumcol = 0.0;
136  double sumsqcol = 0.0;
137  double sumdsqcol = 0.0;
138  double sumind = 0.0;
139  double sumsqind = 0.0;
140  double sumdsqind = 0.0;
141  unsigned int nchan = 0;
142  unsigned int nchancol = 0;
143  unsigned int nchanind = 0;
144  unsigned int count = 0;
145  unsigned int countcol = 0;
146  unsigned int countind = 0;
147  // Loop over channels.
148  cout << myname << "Looping over channels." << endl;
149  for ( unsigned int icha=0; icha<1000000; ++icha ) {
150  geo::SigType_t sigtyp = hgeosvc->SignalType(icha);
151  if ( sigtyp == geo::kMysteryType ) break;
152  if ( icha >= maxchan ) break;
153  ++nchan;
154  bool firstcol = false;
155  bool firstind = false;
156  string labtyp = "unknown";
157  if ( sigtyp == geo::kCollection ) {
158  firstcol = nchancol == 0;
159  labtyp = "collection";
160  ++nchancol;
161  } else if ( sigtyp == geo::kInduction ) {
162  firstind = nchanind == 0;
163  labtyp = "induction";
164  ++nchanind;
165  } else abort();
166  AdcSignalVector sigs = sigsin;
167  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
168  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
169  assert( noisvc.addNoise(clockData, detProp, icha, sigs) == 0 );
170  if ( firstcol || firstind ) {
171  cout << myname << "First " << labtyp << " channel: " << icha << endl;
172  }
173  for ( unsigned int isig=0; isig<ntick; ++isig ) {
174  if ( firstcol || firstind ) {
175  if ( isig < 20 || ntick-isig < 20 ) {
176  cout << myname << setw(5) << isig << ": " << setw(8) << sigs[isig] << endl;
177  }
178  if ( isig == 20 && ntick-isig > 20 ) {
179  cout << myname << setw(5) << "..." << endl;
180  }
181  }
182  float sig = sigs[isig];
183  double dif = sig - ped;
184  sum += sig;
185  sumsq += sig*sig;
186  sumdsq += dif*dif;
187  ++count;
188  if ( sigtyp == geo::kCollection ) {
189  sumcol += sig;
190  sumsqcol += sig*sig;
191  sumdsqcol += dif*dif;
192  ++countcol;
193  } else if ( sigtyp == geo::kInduction ) {
194  sumind += sig;
195  sumsqind += sig*sig;
196  sumdsqind += dif*dif;
197  ++countind;
198  }
199  }
200  }
201  cout << myname << " # channels: " << nchan << endl;
202  cout << myname << "# collection channels: " << nchancol << endl;
203  cout << myname << " # induction channels: " << nchanind << endl;
204  double nticktot = double(nchan)*ntick;
205  double nticktotcol = double(nchancol)*ntick;
206  double nticktotind = double(nchanind)*ntick;
207  cout << myname << " # tick total: " << nticktot << endl;
208  cout << myname << " Count: " << count << endl;
209  cout << myname << " Collection # tick total: " << nticktotcol << endl;
210  cout << myname << " Collection Count: " << countcol << endl;
211  cout << myname << " Induction # tick total: " << nticktotind << endl;
212  cout << myname << " Induction Count: " << countind << endl;
213  assert( nchan > 0 );
214  float mean = sum/nticktot;
215  float rms = sqrt(sumsq/nticktot - mean*mean);
216  float rms2 = sqrt(sumdsq/nticktot);
217  float meancol = sumcol/nticktotcol;
218  float rmscol = sqrt(sumsqcol/nticktotcol - meancol*meancol);
219  float rms2col = sqrt(sumdsqcol/nticktotcol);
220  float meanind = sumind/nticktotind;
221  float rmsind = sqrt(sumsqind/nticktotind - meanind*meanind);
222  float rms2ind = sqrt(sumdsqind/nticktotind);
223  cout << myname << " Mean: " << mean << endl;
224  cout << myname << " RMS: " << rms << endl;
225  cout << myname << " RMS2: " << rms2 << endl;
226  cout << myname << " Collection Mean: " << meancol << endl;
227  cout << myname << " Collection RMS: " << rmscol << endl;
228  cout << myname << " Collection RMS2: " << rms2col << endl;
229  cout << myname << " Induction Mean: " << meanind << endl;
230  cout << myname << " Induction RMS: " << rmsind << endl;
231  cout << myname << " Induction RMS2: " << rms2ind << endl;
232 
233  cout << myname << "Done." << endl;
234  return 0;
235 }
236 
237 int main(int argc, char* argv[]) {
238  int fftsize = 0;
239  int maxchan = 2000;
240  int ntick = 1000;
241  if ( argc > 1 ) {
242  istringstream ssarg(argv[1]);
243  ssarg >> fftsize;
244  }
245  if ( argc > 2 ) {
246  istringstream ssarg(argv[2]);
247  ssarg >> maxchan;
248  }
249  load_services(fftsize);
250  test_ExponentialChannelNoiseService(ntick, maxchan);
251  // We must close TFileService to write out histograms.
252  cout << "Close services." << endl;
253  return 0;
254 }
255 
256 
257 //******************************************************************************
258 //*** Big configuration chunks
259 //***
260 
261 // this is a copy of the "standard" LArProperties configuration
263  service_provider: "LArPropertiesServiceStandard"
264 
265  # For following parameters, see http://pdg.lbl.gov/AtomicNuclearProperties/
266  RadiationLength: 19.55 # g/cm^2
267  AtomicNumber: 18 # Ar atomic number.
268  AtomicMass: 39.948 # Ar atomic mass (g/mol).
269  ExcitationEnergy: 188.0 # Ar mean excitation energy (eV).
270 
271 # realistic Argon 39 decays
272 # Argon39DecayRate: 0.00141 # decays per cm^3 per second. Assumes 1.01 Bq/kg and a density of 1.396 g/cc
273 # switch them off for faster simulation
274  Argon39DecayRate: 0.0
275 
276  # Optical properties
277  # Fast and slow scintillation emission spectra, from [J Chem Phys vol 91 (1989) 1469]
278  FastScintEnergies: [ 6.0, 6.7, 7.1, 7.4, 7.7, 7.9, 8.1, 8.4, 8.5, 8.6, 8.8, 9.0, 9.1, 9.4, 9.8, 10.4, 10.7]
279  SlowScintEnergies: [ 6.0, 6.7, 7.1, 7.4, 7.7, 7.9, 8.1, 8.4, 8.5, 8.6, 8.8, 9.0, 9.1, 9.4, 9.8, 10.4, 10.7]
280  FastScintSpectrum: [ 0.0, 0.04, 0.12, 0.27, 0.44, 0.62, 0.80, 0.91, 0.92, 0.85, 0.70, 0.50, 0.31, 0.13, 0.04, 0.01, 0.0]
281  SlowScintSpectrum: [ 0.0, 0.04, 0.12, 0.27, 0.44, 0.62, 0.80, 0.91, 0.92, 0.85, 0.70, 0.50, 0.31, 0.13, 0.04, 0.01, 0.0]
282  ScintResolutionScale: 1. # resolution factor used by G4 scintillation
283  ScintFastTimeConst: 6. # fast scintillation time constant (ns)
284  ScintSlowTimeConst: 1590. # slow scintillation time constant (ns)
285  ScintBirksConstant: 0.069 # birks constant for LAr (1/MeV cm)
286  ScintYield: 24000. # total scintillation yield (ph/Mev)
287  ScintPreScale: 0.03 # Scale factor to reduce number of photons simulated
288  # Later QE should be corrected for this scale
289  ScintYieldRatio: 0.3 # fast / slow scint ratio (needs revisitting)
290  ScintByParticleType: false # whether to use different yields and
291  # quenching per particle in fast op sim
292  EnableCerenkovLight: true # whether to switch on cerenkov light (slow)
293 
294 
295 
296 
297  # Scintillation yields and fast/slow ratios per particle type
298  MuonScintYield: 24000
299  MuonScintYieldRatio: 0.23
300  PionScintYield: 24000
301  PionScintYieldRatio: 0.23
302  ElectronScintYield: 20000
303  ElectronScintYieldRatio: 0.27
304  KaonScintYield: 24000
305  KaonScintYieldRatio: 0.23
306  ProtonScintYield: 19200
307  ProtonScintYieldRatio: 0.29
308  AlphaScintYield: 16800
309  AlphaScintYieldRatio: 0.56
310 
311 
312  # Refractive index as a function of energy (eV) from arXiv:1502.04213v1
313  RIndexEnergies: [ 1.900, 2.934, 3.592, 5.566, 6.694, 7.540, 8.574, 9.044, 9.232, 9.420, 9.514, 9.608, 9.702, 9.796, 9.890, 9.984, 10.08, 10.27, 10.45, 10.74, 10.92 ]
314  RIndexSpectrum: [ 1.232, 1.236, 1.240, 1.261, 1.282, 1.306, 1.353, 1.387, 1.404, 1.423, 1.434, 1.446, 1.459, 1.473, 1.488, 1.505, 1.524, 1.569, 1.627, 1.751, 1.879 ]
315 
316  # absorption length as function of energy
317  AbsLengthEnergies: [ 4, 5, 6, 7, 8, 9, 10, 11 ]
318  AbsLengthSpectrum: [ 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.]
319 
320  # Rayleigh scattering length (cm) @ 90K as a function of energy (eV) from arXiv:1502.04213
321  RayleighEnergies: [ 2.80, 3.00, 3.50, 4.00, 5.00, 6.00, 7.00, 8.00, 8.50, 9.00, 9.20, 9.40, 9.50, 9.60, 9.70, 9.80, 9.90, 10.0, 10.2, 10.4, 10.6, 10.8 ]
322  RayleighSpectrum: [ 47923., 35981., 18825., 10653., 3972., 1681., 750.9, 334.7, 216.8, 135.0, 109.7, 88.06, 78.32, 69.34, 61.06, 53.46, 46.50, 40.13, 28.91, 19.81, 12.61, 7.20 ]
323 
324  # Surface reflectivity data - vector of energy spectrum per
325  # surface type
326  ReflectiveSurfaceEnergies: [ 7, 9, 10 ]
327  ReflectiveSurfaceNames: [ "STEEL_STAINLESS_Fe7Cr2Ni" ]
328  ReflectiveSurfaceReflectances: [ [ 0.25, 0.25, 0.25 ] ]
329  ReflectiveSurfaceDiffuseFractions: [ [ 0.5, 0.5, 0.5 ] ]
330 
331  # Information related with the simulation of the Wavelength Shifter (TPB)
332  LoadExtraMatProperties: false
333 
334  # TPB - WLS
335  TpbTimeConstant: 2.5 #wls time constant in s J. Lumin 81(1999) 285
336 
337  # WLS - TPB properties original tpb [0.0, 0.0, 0.0, 0.0588,0.235, 0.853, 1.0,1.0,0.9259,0.704,0.0296,0.011, 0.0,0.0, 0.]
338  TpbEmmisionEnergies: [0.05,1.0,1.5, 2.25, 2.481, 2.819, 2.952,2.988,3.024, 3.1, 3.14,3.1807, 3.54, 5.5, 50.39]
339  TpbEmmisionSpectrum: [0.0, 0.0, 0.0, 0.0588,0.235, 0.853, 1.0,1.0,0.9259,0.704,0.0296,0.011, 0.0,0.0, 0.]
340  TpbAbsorptionEnergies: [0.05,1.77,2.0675, 7.42, 7.75, 8.16, 8.73, 9.78,10.69, 50.39]
341  TpbAbsorptionSpectrum: [100000.0,100000.0, 100000.0,0.001,0.00000000001,0.00000000001, 0.00000000001, 0.00000000001, 0.00000000001, 0.00000000001]
342 )cfg"};
343 
344 
345 // this is a copy of the "standard" DetectorProperties configuration
347 service_provider: "DetectorPropertiesServiceStandard"
348 
349 # Drift properties
350 SternheimerA: 0.1956 # Ar Sternheimer parameter a.
351 SternheimerK: 3.0000 # Ar Sternheimer parameter k.
352 SternheimerX0: 0.2000 # Ar Sternheimer parameter x0.
353 SternheimerX1: 3.0000 # Ar Sternheimer parameter x0.
354 SternheimerCbar: 5.2146 # Ar Sternheimer parameter Cbar.
355 
356 Temperature: 87
357 Electronlifetime: 3.0e3
358 Efield: [0.5,0.666,0.8] #(predicted for microBooNE)
359 ElectronsToADC: 6.8906513e-3 # 1fC = 43.008 ADC counts for DUNE fd
360 NumberTimeSamples: 4492 # drift length/drift velocity*sampling rate = (359.4 cm)/(0.16 cm/us)*(2 MHz)
361 ReadOutWindowSize: 4492 # drift length/drift velocity*sampling rate = (359.4 cm)/(0.16 cm/us)*(2 MHz)
362 TimeOffsetU: 0.
363 TimeOffsetV: 0.
364 TimeOffsetZ: 0.
365 DriftVelFudgeFactor: 1.
366 
367 SimpleBoundaryProcess: true #enable opticalBoundaryProcessSimple instead of G4 default
368 
369 )cfg"};
int main(int argc, char *argv[])
Who knows?
Definition: geo_types.h:147
double rms(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:40
std::string string
Definition: nybbler.cc:12
float AdcSignal
Definition: AdcTypes.h:21
const std::string DetectorPropertiesServiceConfigurationString
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
static void load_services(std::string const &config)
art framework interface to geometry description
int addNoise(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, Channel chan, AdcSignalVector &sigs) const
int test_ExponentialChannelNoiseService(unsigned int ntick, unsigned int maxchan=100)
int FFTSize() const
Definition: LArFFT.h:69
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Signal from induction planes.
Definition: geo_types.h:145
static Config * config
Definition: config.cpp:1054
enum geo::_plane_sigtype SigType_t
std::vector< short > AdcVector
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
const std::string LArPropertiesServiceConfigurationString
void line(double t, double *p, double &x, double &y, double &z)
std::vector< AdcSignal > AdcSignalVector
Definition: AdcTypes.h:22
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
void put(std::string const &key)
QTextStream & endl(QTextStream &s)
Signal from collection planes.
Definition: geo_types.h:146