Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
Sensor.cc
Go to the documentation of this file.
1#include <iostream>
2#include <fstream>
3#include <cmath>
4
5#include "Sensor.hh"
7#include "Plotting.hh"
8#include "Numerics.hh"
10
11namespace Garfield {
12
13double Sensor::m_signalConversion = ElementaryCharge;
14
16 : m_nComponents(0),
17 m_lastComponent(-1),
18 m_nElectrodes(0),
19 m_nTimeBins(200),
20 m_tStart(0.),
21 m_tStep(10.),
22 m_nEvents(0),
23 m_hasTransferFunction(false),
24 m_fTransfer(0),
25 m_hasNoiseFunction(false),
26 m_fNoise(0),
27 m_nThresholdCrossings(0),
28 m_hasUserArea(false),
29 m_xMinUser(0.),
30 m_yMinUser(0.),
31 m_zMinUser(0.),
32 m_xMaxUser(0.),
33 m_yMaxUser(0.),
34 m_zMaxUser(0.),
35 m_debug(false) {
36
37 m_className = "Sensor";
38
39 m_components.clear();
40 m_electrodes.clear();
41 m_thresholdCrossings.clear();
42}
43
44void Sensor::ElectricField(const double x, const double y, const double z,
45 double& ex, double& ey, double& ez, double& v,
46 Medium*& medium, int& status) {
47
48 ex = ey = ez = v = 0.;
49 status = -10;
50 medium = 0;
51 double fx, fy, fz, p;
52 Medium* med = 0;
53 int stat;
54 // Add up electric field contributions from all components.
55 for (int i = m_nComponents; i--;) {
56 m_components[i].comp->ElectricField(x, y, z, fx, fy, fz, p, med, stat);
57 if (status != 0) {
58 status = stat;
59 medium = med;
60 }
61 if (stat == 0) {
62 ex += fx;
63 ey += fy;
64 ez += fz;
65 v += p;
66 }
67 }
68}
69
70void Sensor::ElectricField(const double x, const double y, const double z,
71 double& ex, double& ey, double& ez, Medium*& medium,
72 int& status) {
73
74 ex = ey = ez = 0.;
75 status = -10;
76 medium = 0;
77 double fx, fy, fz;
78 Medium* med = 0;
79 int stat;
80 // Add up electric field contributions from all components.
81 for (int i = m_nComponents; i--;) {
82 m_components[i].comp->ElectricField(x, y, z, fx, fy, fz, med, stat);
83 if (status != 0) {
84 status = stat;
85 medium = med;
86 }
87 if (stat == 0) {
88 ex += fx;
89 ey += fy;
90 ez += fz;
91 }
92 }
93}
94
95void Sensor::MagneticField(const double x, const double y, const double z,
96 double& bx, double& by, double& bz, int& status) {
97
98 bx = by = bz = 0.;
99 double fx, fy, fz;
100 // Add up contributions.
101 for (int i = m_nComponents; i--;) {
102 m_components[i].comp->MagneticField(x, y, z, fx, fy, fz, status);
103 if (status != 0) continue;
104 bx += fx;
105 by += fy;
106 bz += fz;
107 }
108}
109
110void Sensor::WeightingField(const double x, const double y, const double z,
111 double& wx, double& wy, double& wz,
112 const std::string label) {
113
114 wx = wy = wz = 0.;
115 double fx = 0., fy = 0., fz = 0.;
116 // Add up field contributions from all components.
117 for (int i = m_nElectrodes; i--;) {
118 if (m_electrodes[i].label == label) {
119 fx = fy = fz = 0.;
120 m_electrodes[i].comp->WeightingField(x, y, z, fx, fy, fz, label);
121 wx += fx;
122 wy += fy;
123 wz += fz;
124 }
125 }
126}
127
128double Sensor::WeightingPotential(const double x, const double y,
129 const double z, const std::string label) {
130
131 double v = 0.;
132 // Add up contributions from all components.
133 for (int i = m_nElectrodes; i--;) {
134 if (m_electrodes[i].label == label) {
135 v += m_electrodes[i].comp->WeightingPotential(x, y, z, label);
136 }
137 }
138 return v;
139}
140
141bool Sensor::GetMedium(const double x, const double y, const double z,
142 Medium*& m) {
143
144 m = NULL;
145
146 // Make sure there is at least one component.
147 if (m_lastComponent < 0) return false;
148
149 // Check if we are still in the same component as in the previous call.
150 m = m_components[m_lastComponent].comp->GetMedium(x, y, z);
151 // Cross-check that the medium is defined.
152 if (m) return true;
153
154 for (int i = m_nComponents; i--;) {
155 m = m_components[i].comp->GetMedium(x, y, z);
156 // Cross-check that the medium is defined.
157 if (m) {
158 m_lastComponent = i;
159 return true;
160 }
161 }
162 return false;
163}
164
166
167 if (!GetBoundingBox(m_xMinUser, m_yMinUser, m_zMinUser, m_xMaxUser,
168 m_yMaxUser, m_zMaxUser)) {
169 std::cerr << m_className << "::SetArea:\n";
170 std::cerr << " Bounding box is not known.\n";
171 return false;
172 }
173
174 std::cout << m_className << "::SetArea:\n";
175 std::cout << " " << m_xMinUser << " < x [cm] < " << m_xMaxUser << "\n";
176 std::cout << " " << m_yMinUser << " < y [cm] < " << m_yMaxUser << "\n";
177 std::cout << " " << m_zMinUser << " < z [cm] < " << m_zMaxUser << "\n";
178 if (std::isinf(m_xMinUser) || std::isinf(m_xMaxUser)) {
179 std::cerr << m_className << "::SetArea:\n";
180 std::cerr << " Warning: infinite x-range\n";
181 }
182 if (std::isinf(m_yMinUser) || std::isinf(m_yMaxUser)) {
183 std::cerr << m_className << "::SetArea:\n";
184 std::cerr << " Warning: infinite x-range\n";
185 }
186 if (std::isinf(m_zMinUser) || std::isinf(m_zMaxUser)) {
187 std::cerr << m_className << "::SetArea:\n";
188 std::cerr << " Warning: infinite x-range\n";
189 }
190 m_hasUserArea = true;
191 return true;
192}
193
194bool Sensor::SetArea(const double xmin, const double ymin, const double zmin,
195 const double xmax, const double ymax, const double zmax) {
196
197 if (fabs(xmax - xmin) < Small || fabs(ymax - ymin) < Small ||
198 fabs(zmax - zmin) < Small) {
199 std::cerr << m_className << "::SetArea:\n";
200 std::cerr << " Invalid range.\n";
201 return false;
202 }
203
204 m_xMinUser = xmin;
205 m_yMinUser = ymin;
206 m_zMinUser = zmin;
207 m_xMaxUser = xmax;
208 m_yMaxUser = ymax;
209 m_zMaxUser = zmax;
210
211 if (xmin > xmax) {
212 m_xMinUser = xmax;
213 m_xMaxUser = xmin;
214 }
215 if (ymin > ymax) {
216 m_yMinUser = ymax;
217 m_yMaxUser = ymin;
218 }
219 if (zmin > zmax) {
220 m_zMinUser = zmax;
221 m_zMaxUser = zmin;
222 }
223 m_hasUserArea = true;
224 return true;
225}
226
227bool Sensor::GetArea(double& xmin, double& ymin, double& zmin, double& xmax,
228 double& ymax, double& zmax) {
229
230 if (m_hasUserArea) {
231 xmin = m_xMinUser;
232 ymin = m_yMinUser;
233 zmin = m_zMinUser;
234 xmax = m_xMaxUser;
235 ymax = m_yMaxUser;
236 zmax = m_zMaxUser;
237 return true;
238 }
239
240 // User area bounds are not (yet) defined.
241 // Get the bounding box of the sensor.
242 if (!SetArea()) return false;
243
244 xmin = m_xMinUser;
245 ymin = m_yMinUser;
246 zmin = m_zMinUser;
247 xmax = m_xMaxUser;
248 ymax = m_yMaxUser;
249 zmax = m_zMaxUser;
250
251 return true;
252}
253
254bool Sensor::IsInArea(const double x, const double y, const double z) {
255
256 if (!m_hasUserArea) {
257 if (!SetArea()) {
258 std::cerr << m_className << "::IsInArea:\n";
259 std::cerr << " User area cannot be established.\n";
260 return false;
261 }
262 m_hasUserArea = true;
263 }
264
265 if (x >= m_xMinUser && x <= m_xMaxUser && y >= m_yMinUser &&
266 y <= m_yMaxUser && z >= m_zMinUser && z <= m_zMaxUser) {
267 return true;
268 }
269
270 if (m_debug) {
271 std::cout << m_className << "::IsInArea:\n" << std::endl;
272 std::cout << " (" << x << ", " << y << ", " << z << ") "
273 << " is outside.\n";
274 }
275 return false;
276}
277
278bool Sensor::IsWireCrossed(const double x0, const double y0, const double z0,
279 const double x1, const double y1, const double z1,
280 double& xc, double& yc, double& zc) {
281
282 for (int i = m_nComponents; i--;) {
283 if (m_components[i].comp->IsWireCrossed(x0, y0, z0,
284 x1, y1, z1, xc, yc, zc)) {
285 return true;
286 }
287 }
288 return false;
289}
290
291bool Sensor::IsInTrapRadius(double x0, double y0, double z0, double& xw,
292 double& yw, double& rw) {
293
294 for (int i = m_nComponents; i--;) {
295 if (m_components[i].comp->IsInTrapRadius(x0, y0, z0, xw, yw, rw)) {
296 return true;
297 }
298 }
299 return false;
300}
301
303
304 if (!comp) {
305 std::cerr << m_className << "::AddComponent:\n";
306 std::cerr << " Component pointer is null.\n";
307 return;
308 }
309
310 component newComponent;
311 newComponent.comp = comp;
312 m_components.push_back(newComponent);
313 ++m_nComponents;
314 if (m_nComponents == 1) m_lastComponent = 0;
315}
316
317void Sensor::AddElectrode(ComponentBase* comp, std::string label) {
318
319 if (!comp) {
320 std::cerr << m_className << "::AddElectrode:\n";
321 std::cerr << " Component pointer is null.\n";
322 return;
323 }
324
325 for (int i = m_nElectrodes; i--;) {
326 if (m_electrodes[i].label == label) {
327 std::cout << m_className << "::AddElectrode:\n";
328 std::cout << " Warning: An electrode with label \"" << label
329 << "\" exists already.\n";
330 std::cout << " Weighting fields will be summed up.\n";
331 break;
332 }
333 }
334
335 electrode newElectrode;
336 newElectrode.comp = comp;
337 newElectrode.label = label;
338 m_electrodes.push_back(newElectrode);
339 ++m_nElectrodes;
340 m_electrodes[m_nElectrodes - 1].signal.resize(m_nTimeBins);
341 m_electrodes[m_nElectrodes - 1].electronsignal.resize(m_nTimeBins);
342 m_electrodes[m_nElectrodes - 1].ionsignal.resize(m_nTimeBins);
343 std::cout << m_className << "::AddElectrode:\n";
344 std::cout << " Added readout electrode \"" << label << "\".\n";
345 std::cout << " All signals are reset.\n";
346 ClearSignal();
347}
348
350
351 m_components.clear();
352 m_nComponents = 0;
353 m_lastComponent = -1;
354 m_electrodes.clear();
355 m_nElectrodes = 0;
356 m_nTimeBins = 200;
357 m_tStart = 0.;
358 m_tStep = 10.;
359 m_nEvents = 0;
360 m_hasUserArea = false;
361}
362
363bool Sensor::GetVoltageRange(double& vmin, double& vmax) {
364
365 // We don't know the range yet.
366 bool set = false;
367 // Loop over the components.
368 double umin, umax;
369 for (int i = 0; i < m_nComponents; ++i) {
370 if (!m_components[i].comp->GetVoltageRange(umin, umax)) continue;
371 if (set) {
372 if (umin < vmin) vmin = umin;
373 if (umax > vmax) vmax = umax;
374 } else {
375 vmin = umin;
376 vmax = umax;
377 set = true;
378 }
379 }
380
381 // Warn if we still don't know the range.
382 if (!set) {
383 std::cerr << m_className << "::GetVoltageRange:\n";
384 std::cerr << " Sensor voltage range not known.\n";
385 vmin = vmax = 0.;
386 return false;
387 }
388
389 if (m_debug) {
390 std::cout << m_className << "::GetVoltageRange:\n";
391 std::cout << " Voltage range " << vmin << " < V < " << vmax << ".\n";
392 }
393 return true;
394}
395
397
398 for (int i = m_nElectrodes; i--;) {
399 m_electrodes[i].charge = 0.;
400 for (int j = m_nTimeBins; j--;) {
401 m_electrodes[i].signal[j] = 0.;
402 m_electrodes[i].electronsignal[j] = 0.;
403 m_electrodes[i].ionsignal[j] = 0.;
404 }
405 }
406 m_nEvents = 0;
407}
408
409void Sensor::AddSignal(const double& q, const double& t, const double& dt,
410 const double& x, const double& y, const double& z,
411 const double& vx, const double& vy, const double& vz) {
412
413 // Get the time bin.
414 if (t < m_tStart || dt <= 0.) {
415 if (m_debug) {
416 std::cerr << m_className << "::AddSignal:\n";
417 if (t < m_tStart) std::cerr << " Time " << t << " out of range.\n";
418 if (dt <= 0.) std::cerr << " Time step < 0.\n";
419 }
420 return;
421 }
422 const int bin = int((t - m_tStart) / m_tStep);
423 // Check if the starting time is outside the range
424 if (bin < 0 || bin >= m_nTimeBins) {
425 if (m_debug) {
426 std::cerr << m_className << "::AddSignal:\n";
427 std::cerr << " Bin " << bin << " out of range.\n";
428 }
429 return;
430 }
431 if (m_nEvents <= 0) m_nEvents = 1;
432
433 double wx = 0., wy = 0., wz = 0.;
434 if (m_debug) {
435 std::cout << m_className << "::AddSignal:\n";
436 std::cout << " Time: " << t << "\n";
437 std::cout << " Step: " << dt << "\n";
438 std::cout << " Charge: " << q << "\n";
439 std::cout << " Velocity: (" << vx << ", " << vy << ", " << vz << ")\n";
440 }
441 for (int i = m_nElectrodes; i--;) {
442 // Calculate the weighting field for this electrode
443 m_electrodes[i].comp->WeightingField(x, y, z, wx, wy, wz,
444 m_electrodes[i].label);
445 // Calculate the induced current
446 const double cur = -q * (wx * vx + wy * vy + wz * vz);
447 if (m_debug) {
448 std::cout << " Electrode " << m_electrodes[i].label << ":\n";
449 std::cout << " Weighting field: (" << wx << ", " << wy << ", " << wz
450 << ")\n";
451 std::cout << " Induced charge: " << cur* dt << "\n";
452 }
453 double delta = m_tStart + (bin + 1) * m_tStep - t;
454 // Check if the provided timestep extends over more than one time bin
455 if (dt > delta) {
456 m_electrodes[i].signal[bin] += cur * delta;
457 if (q < 0) {
458 m_electrodes[i].electronsignal[bin] += cur * delta;
459 } else {
460 m_electrodes[i].ionsignal[bin] += cur * delta;
461 }
462 delta = dt - delta;
463 int j = 1;
464 while (delta > m_tStep && bin + j < m_nTimeBins) {
465 m_electrodes[i].signal[bin + j] += cur * m_tStep;
466 if (q < 0) {
467 m_electrodes[i].electronsignal[bin + j] += cur * m_tStep;
468 } else {
469 m_electrodes[i].ionsignal[bin + j] += cur * m_tStep;
470 }
471 delta -= m_tStep;
472 ++j;
473 }
474 if (bin + j < m_nTimeBins) {
475 m_electrodes[i].signal[bin + j] += cur * delta;
476 if (q < 0) {
477 m_electrodes[i].electronsignal[bin + j] += cur * delta;
478 } else {
479 m_electrodes[i].ionsignal[bin + j] += cur * delta;
480 }
481 }
482 } else {
483 m_electrodes[i].signal[bin] += cur * dt;
484 if (q < 0) {
485 m_electrodes[i].electronsignal[bin] += cur * dt;
486 } else {
487 m_electrodes[i].ionsignal[bin] += cur * dt;
488 }
489 }
490 }
491}
492
493void Sensor::AddInducedCharge(const double q, const double x0, const double y0,
494 const double z0, const double x1, const double y1,
495 const double z1) {
496
497 if (m_debug) std::cout << m_className << "::AddInducedCharge:\n";
498 double w0 = 0., w1 = 0.;
499 for (int i = m_nElectrodes; i--;) {
500 // Calculate the weighting potential at the starting point.
501 w0 = m_electrodes[i]
502 .comp->WeightingPotential(x0, y0, z0, m_electrodes[i].label);
503 // Calculate the weighting potential at the end point.
504 w1 = m_electrodes[i]
505 .comp->WeightingPotential(x1, y1, z1, m_electrodes[i].label);
506 m_electrodes[i].charge += q * (w1 - w0);
507 if (m_debug) {
508 std::cout << " Electrode " << m_electrodes[i].label << ":\n";
509 std::cout << " Weighting potential at (" << x0 << ", " << y0 << ", "
510 << z0 << "): " << w0 << "\n";
511 std::cout << " Weighting potential at (" << x1 << ", " << y1 << ", "
512 << z1 << "): " << w1 << "\n";
513 std::cout << " Induced charge: " << m_electrodes[i].charge << "\n";
514 }
515 }
516}
517
518void Sensor::SetTimeWindow(const double tstart, const double tstep,
519 const int nsteps) {
520
521 m_tStart = tstart;
522 if (tstep <= 0.) {
523 std::cerr << m_className << "::SetTimeWindow:\n";
524 std::cerr << " Starting time out of range.\n";
525 } else {
526 m_tStep = tstep;
527 }
528
529 if (nsteps <= 0) {
530 std::cerr << m_className << "::SetTimeWindow:\n";
531 std::cerr << " Number of time bins out of range.\n";
532 } else {
533 m_nTimeBins = nsteps;
534 }
535
536 if (m_debug) {
537 std::cout << m_className << "::SetTimeWindow:\n";
538 std::cout << " " << m_tStart << " < t [ns] < "
539 << m_tStart + m_nTimeBins* m_tStep << "\n";
540 std::cout << " Step size: " << m_tStep << " ns\n";
541 }
542
543 std::cout << m_className << "::SetTimeWindow:\n";
544 std::cout << " Resetting all signals.\n";
545 for (int i = m_nElectrodes; i--;) {
546 m_electrodes[i].signal.clear();
547 m_electrodes[i].signal.resize(m_nTimeBins);
548 m_electrodes[i].electronsignal.clear();
549 m_electrodes[i].electronsignal.resize(m_nTimeBins);
550 m_electrodes[i].ionsignal.clear();
551 m_electrodes[i].ionsignal.resize(m_nTimeBins);
552 }
553 m_nEvents = 0;
554}
555
556double Sensor::GetElectronSignal(const std::string label, const int bin) {
557
558 if (m_nEvents <= 0) return 0.;
559 if (bin < 0 || bin >= m_nTimeBins) return 0.;
560 double sig = 0.;
561 for (int i = m_nElectrodes; i--;) {
562 if (m_electrodes[i].label == label)
563 sig += m_electrodes[i].electronsignal[bin];
564 }
565 if (m_debug) {
566 std::cout << m_className << "::GetElectronSignal:\n";
567 std::cout << " Electrode: " << label << "\n";
568 std::cout << " Bin: " << bin << "\n";
569 std::cout << " ElectronSignal: " << sig / m_tStep << "\n";
570 }
571 return m_signalConversion * sig / (m_nEvents * m_tStep);
572}
573
574double Sensor::GetIonSignal(const std::string label, const int bin) {
575
576 if (m_nEvents <= 0) return 0.;
577 if (bin < 0 || bin >= m_nTimeBins) return 0.;
578 double sig = 0.;
579 for (int i = m_nElectrodes; i--;) {
580 if (m_electrodes[i].label == label) sig += m_electrodes[i].ionsignal[bin];
581 }
582 if (m_debug) {
583 std::cout << m_className << "::GetIonSignal:\n";
584 std::cout << " Electrode: " << label << "\n";
585 std::cout << " Bin: " << bin << "\n";
586 std::cout << " IonSignal: " << sig / m_tStep << "\n";
587 }
588 return m_signalConversion * sig / (m_nEvents * m_tStep);
589}
590
591double Sensor::GetSignal(const std::string label, const int bin) {
592
593 if (m_nEvents <= 0) return 0.;
594 if (bin < 0 || bin >= m_nTimeBins) return 0.;
595 double sig = 0.;
596 for (int i = m_nElectrodes; i--;) {
597 if (m_electrodes[i].label == label) sig += m_electrodes[i].signal[bin];
598 }
599 if (m_debug) {
600 std::cout << m_className << "::GetSignal:\n";
601 std::cout << " Electrode: " << label << "\n";
602 std::cout << " Bin: " << bin << "\n";
603 std::cout << " Signal: " << sig / m_tStep << "\n";
604 }
605 return m_signalConversion * sig / (m_nEvents * m_tStep);
606}
607
608double Sensor::GetInducedCharge(const std::string label) {
609
610 if (m_nEvents <= 0) return 0.;
611 double charge = 0.;
612 for (int i = m_nElectrodes; i--;) {
613 if (m_electrodes[i].label == label) charge += m_electrodes[i].charge;
614 }
615 if (m_debug) {
616 std::cout << m_className << "::GetInducedCharge:\n";
617 std::cout << " Electrode: " << label << "\n";
618 std::cout << " Charge: " << charge / m_tStep << "\n";
619 }
620
621 return charge / m_nEvents;
622}
623
624void Sensor::SetTransferFunction(double (*f)(double t)) {
625
626 if (!f) {
627 std::cerr << m_className << "::SetTransferFunction:\n";
628 std::cerr << " Function pointer is null.\n";
629 return;
630 }
631 m_fTransfer = f;
632 m_hasTransferFunction = true;
633 m_transferFunctionTimes.clear();
634 m_transferFunctionValues.clear();
635}
636
637void Sensor::SetTransferFunction(std::vector<double> times,
638 std::vector<double> values) {
639
640 if (times.empty() || values.empty()) {
641 std::cerr << m_className << "::SetTransferFunction:\n";
642 std::cerr << " Time and value vectors must not be empty.\n";
643 return;
644 } else if (times.size() != values.size()) {
645 std::cerr << m_className << "::SetTransferFunction:\n";
646 std::cerr << " Time and value vectors must have same size.\n";
647 return;
648 }
649 m_transferFunctionTimes = times;
650 m_transferFunctionValues = values;
651 m_fTransfer = 0;
652 m_hasTransferFunction = true;
653}
654
655double Sensor::InterpolateTransferFunctionTable(double t) {
656
657 if (m_transferFunctionTimes.empty() || m_transferFunctionValues.empty()) {
658 return 0.;
659 }
660 // Don't extrapolate beyond the range defined in the table.
661 if (t < m_transferFunctionTimes.front() ||
662 t > m_transferFunctionTimes.back()) {
663 return 0.;
664 }
665 // Find the proper interval in the table.
666 int iLow = 0;
667 int iUp = m_transferFunctionTimes.size() - 1;
668 int iM;
669 while (iUp - iLow > 1) {
670 iM = (iUp + iLow) >> 1;
671 if (t >= m_transferFunctionTimes[iM]) {
672 iLow = iM;
673 } else {
674 iUp = iM;
675 }
676 }
677 // Linear interpolation.
678 return m_transferFunctionValues[iLow] +
679 (t - m_transferFunctionTimes[iLow]) *
680 (m_transferFunctionValues[iUp] - m_transferFunctionValues[iLow]) /
681 (m_transferFunctionTimes[iUp] - m_transferFunctionTimes[iLow]);
682}
683
684double Sensor::GetTransferFunction(const double t) {
685
686 if (!m_hasTransferFunction) return 0.;
687 if (m_fTransfer) return m_fTransfer(t);
688 return InterpolateTransferFunctionTable(t);
689}
690
692
693 if (!m_hasTransferFunction) {
694 std::cerr << m_className << "::ConvoluteSignal:\n";
695 std::cerr << " No transfer function available.\n";
696 return false;
697 }
698 if (m_nEvents <= 0) {
699 std::cerr << m_className << "::ConvoluteSignal:\n";
700 std::cerr << " No signals present.\n";
701 return false;
702 }
703
704 // Set the range where the transfer function is valid.
705 double cnvMin = 0.;
706 double cnvMax = 1.e10;
707
708 std::vector<double> cnvTab;
709 cnvTab.resize(2 * m_nTimeBins);
710 int iOffset = m_nTimeBins;
711 // Evaluate the transfer function.
712 for (int i = 0; i < m_nTimeBins; ++i) {
713 // Negative time part.
714 double t = -i * m_tStep;
715 if (t < cnvMin || t > cnvMax) {
716 cnvTab[iOffset - i] = 0.;
717 } else if (m_fTransfer) {
718 cnvTab[iOffset - i] = m_fTransfer(t);
719 } else {
720 cnvTab[iOffset - i] = InterpolateTransferFunctionTable(t);
721 }
722 if (i < 1) continue;
723 // Positive time part.
724 t = i * m_tStep;
725 if (t < cnvMin || t > cnvMax) {
726 cnvTab[iOffset + i] = 0.;
727 } else if (m_fTransfer) {
728 cnvTab[iOffset + i] = m_fTransfer(t);
729 } else {
730 cnvTab[iOffset + i] = InterpolateTransferFunctionTable(t);
731 }
732 }
733
734 std::vector<double> tmpSignal;
735 tmpSignal.resize(m_nTimeBins);
736 // Loop over all electrodes.
737 for (int i = 0; i < m_nElectrodes; ++i) {
738 for (int j = 0; j < m_nTimeBins; ++j) {
739 tmpSignal[j] = 0.;
740 for (int k = 0; k < m_nTimeBins; ++k) {
741 tmpSignal[j] +=
742 m_tStep * cnvTab[iOffset + j - k] * m_electrodes[i].signal[k];
743 }
744 }
745 for (int j = 0; j < m_nTimeBins; ++j) {
746 m_electrodes[i].signal[j] = tmpSignal[j];
747 }
748 }
749 return true;
750}
751
753
754 if (m_nEvents <= 0) {
755 std::cerr << m_className << "::IntegrateSignal:\n";
756 std::cerr << " No signals present.\n";
757 return false;
758 }
759
760 for (int i = 0; i < m_nElectrodes; ++i) {
761 for (int j = 0; j < m_nTimeBins; ++j) {
762 m_electrodes[i].signal[j] *= m_tStep;
763 m_electrodes[i].electronsignal[j] *= m_tStep;
764 m_electrodes[i].ionsignal[j] *= m_tStep;
765 if (j > 0) {
766 m_electrodes[i].signal[j] += m_electrodes[i].signal[j - 1];
767 m_electrodes[i].electronsignal[j] +=
768 m_electrodes[i].electronsignal[j - 1];
769 m_electrodes[i].ionsignal[j] += m_electrodes[i].ionsignal[j - 1];
770 }
771 }
772 }
773 return true;
774}
775
776void Sensor::SetNoiseFunction(double (*f)(double t)) {
777
778 if (f == 0) {
779 std::cerr << m_className << "::SetNoiseFunction:\n";
780 std::cerr << " Function pointer is null.\n";
781 return;
782 }
783 m_fNoise = f;
784 m_hasNoiseFunction = true;
785}
786
788
789 if (!m_hasNoiseFunction) {
790 std::cerr << m_className << "::AddNoise:\n";
791 std::cerr << " Noise function is not defined.\n";
792 return;
793 }
794 if (m_nEvents <= 0) m_nEvents = 1;
795
796 for (int i = m_nElectrodes; i--;) {
797 for (int j = m_nTimeBins; j--;) {
798 const double t = m_tStart + (j + 0.5) * m_tStep;
799 m_electrodes[i].signal[j] += m_fNoise(t);
800 // Adding noise to both channels might be wrong,
801 // maybe an extended option
802 // where to add noise would be an idea?
803 m_electrodes[i].electronsignal[j] += m_fNoise(t);
804 m_electrodes[i].ionsignal[j] += m_fNoise(t);
805 }
806 }
807}
808
810 const std::string label, int& n) {
811
812 // Reset the list of threshold crossings.
813 m_thresholdCrossings.clear();
814 m_nThresholdCrossings = n = 0;
815 m_thresholdLevel = thr;
816
817 // Set the interpolation order.
818 int iOrder = 1;
819
820 if (m_nEvents <= 0) {
821 std::cerr << m_className << "::ComputeThresholdCrossings:\n";
822 std::cerr << " No signals present.\n";
823 return false;
824 }
825
826 // Compute the total signal.
827 std::vector<double> signal;
828 signal.resize(m_nTimeBins);
829 for (int i = m_nTimeBins; i--;) signal[i] = 0.;
830 // Loop over the electrodes.
831 bool foundLabel = false;
832 for (int j = m_nElectrodes; j--;) {
833 if (m_electrodes[j].label == label) {
834 foundLabel = true;
835 for (int i = m_nTimeBins; i--;) {
836 signal[i] += m_electrodes[j].signal[i];
837 }
838 }
839 }
840 if (!foundLabel) {
841 std::cerr << m_className << "::ComputeThresholdCrossings:\n";
842 std::cerr << " Electrode " << label << " not found.\n";
843 return false;
844 }
845 for (int i = m_nTimeBins; i--;) {
846 signal[i] *= m_signalConversion / (m_nEvents * m_tStep);
847 }
848
849 // Establish the range.
850 double vMin = signal[0];
851 double vMax = signal[0];
852 for (int i = m_nTimeBins; i--;) {
853 if (signal[i] < vMin) vMin = signal[i];
854 if (signal[i] > vMax) vMax = signal[i];
855 }
856 if (thr < vMin && thr > vMax) {
857 if (m_debug) {
858 std::cout << m_className << "::ComputeThresholdCrossings:\n";
859 std::cout << " Threshold outside the range [" << vMin << ", " << vMax
860 << "]\n";
861 }
862 return true;
863 }
864
865 // Check for rising edges.
866 bool rise = true;
867 bool fall = false;
868
869 while (rise || fall) {
870 if (m_debug) {
871 if (rise) {
872 std::cout << m_className << "::ComputeThresholdCrossings:\n";
873 std::cout << " Hunting for rising edges.\n";
874 } else if (fall) {
875 std::cout << m_className << "::ComputeThresholdCrossings:\n";
876 std::cout << " Hunting for falling edges.\n";
877 }
878 }
879 // Initialise the vectors.
880 std::vector<double> times;
881 std::vector<double> values;
882 times.clear();
883 values.clear();
884 times.push_back(m_tStart);
885 values.push_back(signal[0]);
886 int nValues = 1;
887 // Scan the signal.
888 for (int i = 1; i < m_nTimeBins; ++i) {
889 // Compute the vector element.
890 const double tNew = m_tStart + i * m_tStep;
891 const double vNew = signal[i];
892 // If still increasing or decreasing, add to the vector.
893 if ((rise && vNew > values.back()) || (fall && vNew < values.back())) {
894 times.push_back(tNew);
895 values.push_back(vNew);
896 ++nValues;
897 // Otherwise see whether we crossed the threshold level.
898 } else if ((values[0] - thr) * (thr - values.back()) >= 0. &&
899 nValues > 1 && ((rise && values.back() > values[0]) ||
900 (fall && values.back() < values[0]))) {
901 // Compute the crossing time.
902 double tcr = Numerics::Divdif(times, values, nValues, thr, iOrder);
903 thresholdCrossing newCrossing;
904 newCrossing.time = tcr;
905 newCrossing.rise = rise;
906 m_thresholdCrossings.push_back(newCrossing);
907 ++m_nThresholdCrossings;
908 times.clear();
909 values.clear();
910 times.push_back(tNew);
911 values.push_back(vNew);
912 nValues = 1;
913 } else {
914 // No crossing, simply reset the vector.
915 times.clear();
916 values.clear();
917 times.push_back(tNew);
918 values.push_back(vNew);
919 nValues = 1;
920 }
921 }
922 // Check the final vector.
923 if ((values[0] - thr) * (thr - values.back()) >= 0. && nValues > 1 &&
924 ((rise && values.back() > values[0]) ||
925 (fall && values.back() < values[0]))) {
926 double tcr = Numerics::Divdif(times, values, nValues, thr, iOrder);
927 thresholdCrossing newCrossing;
928 newCrossing.time = tcr;
929 newCrossing.rise = rise;
930 m_thresholdCrossings.push_back(newCrossing);
931 ++m_nThresholdCrossings;
932 }
933 if (rise) {
934 rise = false;
935 fall = true;
936 } else if (fall) {
937 rise = fall = false;
938 }
939 }
940 n = m_nThresholdCrossings;
941
942 if (m_debug) {
943 std::cout << m_className << "::ComputeThresholdCrossings:\n";
944 std::cout << " Found " << m_nThresholdCrossings << " crossings.\n";
945 if (m_nThresholdCrossings > 0) {
946 std::cout << " Time [ns] Direction\n";
947 }
948 for (int i = 0; i < m_nThresholdCrossings; ++i) {
949 std::cout << " " << m_thresholdCrossings[i].time << " ";
950 if (m_thresholdCrossings[i].rise) {
951 std::cout << "rising\n";
952 } else {
953 std::cout << "falling\n";
954 }
955 }
956 }
957 // Seems to have worked.
958 return true;
959}
960
961bool Sensor::GetThresholdCrossing(const int i, double& time, double& level,
962 bool& rise) {
963
964 level = m_thresholdLevel;
965
966 if (i < 0 || i >= m_nThresholdCrossings) {
967 std::cerr << m_className << "::GetThresholdCrossing:\n";
968 std::cerr << " Index (" << i << ") out of range.\n";
969 time = m_tStart + m_nTimeBins * m_tStep;
970 return false;
971 }
972
973 time = m_thresholdCrossings[i].time;
974 rise = m_thresholdCrossings[i].rise;
975 return true;
976}
977
978bool Sensor::GetBoundingBox(double& xmin, double& ymin, double& zmin,
979 double& xmax, double& ymax, double& zmax) {
980
981 // We don't know the range yet
982 bool set = false;
983 // Loop over the fields
984 double x0, y0, z0, x1, y1, z1;
985 for (int i = m_nComponents; i--;) {
986 if (!m_components[i].comp->GetBoundingBox(x0, y0, z0, x1, y1, z1)) continue;
987 if (set) {
988 if (x0 < xmin) xmin = x0;
989 if (y0 < ymin) ymin = y0;
990 if (z0 < zmin) zmin = z0;
991 if (x1 > xmax) xmax = x1;
992 if (y1 > ymax) ymax = y1;
993 if (z1 > zmax) zmax = z1;
994 } else {
995 xmin = x0;
996 ymin = y0;
997 zmin = z0;
998 xmax = x1;
999 ymax = y1;
1000 zmax = z1;
1001 set = true;
1002 }
1003 }
1004
1005 // Warn if we still don't know the range
1006 if (!set) {
1007 std::cerr << m_className << "::GetBoundingBox:\n";
1008 std::cerr << " Sensor bounding box not known.\n";
1009 xmin = 0.;
1010 ymin = 0.;
1011 zmin = 0.;
1012 xmax = 0.;
1013 ymax = 0.;
1014 zmax = 0.;
1015 return false;
1016 }
1017
1018 if (m_debug) {
1019 std::cout << m_className << "::GetBoundingBox:\n";
1020 std::cout << " " << xmin << " < x [cm] < " << xmax << "\n";
1021 std::cout << " " << ymin << " < y [cm] < " << ymax << "\n";
1022 std::cout << " " << zmin << " < z [cm] < " << zmax << "\n";
1023 }
1024 return true;
1025}
1026}
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616
void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
Definition: Sensor.cc:95
bool GetVoltageRange(double &vmin, double &vmax)
Definition: Sensor.cc:363
bool ComputeThresholdCrossings(const double thr, const std::string label, int &n)
Definition: Sensor.cc:809
void AddElectrode(ComponentBase *comp, std::string label)
Definition: Sensor.cc:317
bool IntegrateSignal()
Definition: Sensor.cc:752
bool GetThresholdCrossing(const int i, double &time, double &level, bool &rise)
Definition: Sensor.cc:961
double GetTransferFunction(const double t)
Definition: Sensor.cc:684
void AddSignal(const double &q, const double &t, const double &dt, const double &x, const double &y, const double &z, const double &vx, const double &vy, const double &vz)
Definition: Sensor.cc:409
double WeightingPotential(const double x, const double y, const double z, const std::string label)
Definition: Sensor.cc:128
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string label)
Definition: Sensor.cc:110
double GetInducedCharge(const std::string label)
Definition: Sensor.cc:608
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&medium, int &status)
Definition: Sensor.cc:44
bool SetArea()
Definition: Sensor.cc:165
void SetNoiseFunction(double(*f)(double t))
Definition: Sensor.cc:776
bool IsInArea(const double x, const double y, const double z)
Definition: Sensor.cc:254
bool IsInTrapRadius(double x0, double y0, double z0, double &xw, double &yw, double &rw)
Definition: Sensor.cc:291
bool GetMedium(const double x, const double y, const double z, Medium *&medium)
Definition: Sensor.cc:141
void Clear()
Definition: Sensor.cc:349
double GetElectronSignal(const std::string label, const int bin)
Definition: Sensor.cc:556
void SetTransferFunction(double(*f)(double t))
Definition: Sensor.cc:624
bool IsWireCrossed(const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
Definition: Sensor.cc:278
bool ConvoluteSignal()
Definition: Sensor.cc:691
void AddComponent(ComponentBase *comp)
Definition: Sensor.cc:302
void SetTimeWindow(const double tstart, const double tstep, const int nsteps)
Definition: Sensor.cc:518
void AddNoise()
Definition: Sensor.cc:787
void ClearSignal()
Definition: Sensor.cc:396
void AddInducedCharge(const double q, const double x0, const double y0, const double z0, const double x1, const double y1, const double z1)
Definition: Sensor.cc:493
double GetIonSignal(const std::string label, const int bin)
Definition: Sensor.cc:574
double GetSignal(const std::string label, const int bin)
Definition: Sensor.cc:591
bool GetArea(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Definition: Sensor.cc:227
double Divdif(const std::vector< double > &f, const std::vector< double > &a, int nn, double x, int mm)
Definition: Numerics.cc:650