PulseShaper.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Code to shape pulses include the effect of the pre-amplifier
3 /// and the energy dependent fall time.
4 /// \author aurisaam@ucmail.uc.edu
5 /// \date
6 ////////////////////////////////////////////////////////////////////////
7 #include <cmath>
8 #include <iterator>
9 #include <algorithm>
10 
11 #include "PulseShaper.h"
12 
13 namespace rsim{
14 
16  IPulseShaper(pset),
17  fI( pset.get< double >("ASICPreAmpTime")),
18  fR( pset.get< double >("ASICRiseTime")),
19  fFInt( pset.get< double >("ASICFallInt")),
20  fFSlope( pset.get< double >("ASICFallSlope")),
21  fMaxStep( pset.get< double > ("ASICMaxStep") ),
22  fTime(0.0)
23  {
24  fIInv = 1.0/fI;
25  fRInv = 1.0/fR;
26  fInput.resize(2);
27  fPre.resize(7);
28  fRise.resize(4);
29  fRiseDot.resize(3);
30  fShape.resize(2);
31  Reset();
32  }
33 
35  {
36  for (int i = 0; i < 2; ++i) fInput[i] = 0.0;
37  for (int i = 0; i < 7; ++i) fPre[i] = 0.0;
38  for (int i = 0; i < 4; ++i) fRise[i] = 0.0;
39  for (int i = 0; i < 3; ++i) fRiseDot[i] = 0.0;
40  for (int i = 0; i < 2; ++i) fShape[i] = 0.0;
41  }
42 
43  double PulseShaper::Step(double xNext, double dt)
44  {
45  //step forward each variable
46  fInput[0] = fInput[1];
47  fPre[0] = fPre[4];
48  fRise[0] = fRise[2];
49  fRiseDot[0] = fRiseDot[2];
50  fShape[0] = fShape[1];
51  //set the new time step and the next impulse
52  fDt = dt;
53  fTime += fDt;
54  fInput[1] = xNext;
55 
56  double w1(0.0), w2(0.0), w3(0.0), w4(0.0);
57  fPre[1] = fPre[0] - 0.25*fDt*fPre[0]/fI;
58  fPre[2] = fPre[1] - 0.25*fDt*fPre[0]/fI;
59  fPre[3] = fPre[2] - 0.25*fDt*fPre[0]/fI;
60  fPre[4] = fPre[3] - 0.25*fDt*fPre[0]/fI + fInput[1];
61  fPre[5] = fPre[4] - 0.25*fDt*fPre[0]/fI;
62  fPre[6] = fPre[5] - 0.25*fDt*fPre[0]/fI;
63 
64  w1 = 0.5*fDt*(fPre[0] - fRise[0])/fR;
65  w2 = 0.5*fDt*(fPre[1] - fRise[0] - 0.5*w1)/fR;
66  w3 = 0.5*fDt*(fPre[1] - fRise[0] - 0.5*w2)/fR;
67  w4 = 0.5*fDt*(fPre[2] - fRise[0] - w3)/fR;
68  fRise[1] = fRise[0] + (1.0/6.0)*(w1 + 2*w2 + 2*w3 + w4);
69 
70  w1 = 0.5*fDt*(fPre[2] - fRise[1])/fR;
71  w2 = 0.5*fDt*(fPre[3] - fRise[1] - 0.5*w1)/fR;
72  w3 = 0.5*fDt*(fPre[3] - fRise[1] - 0.5*w2)/fR;
73  w4 = 0.5*fDt*(fPre[4] - fRise[1] - w3)/fR;
74  fRise[2] = fRise[1] + (1.0/6.0)*(w1 + 2*w2 + 2*w3 + w4);
75 
76  w1 = 0.5*fDt*(fPre[4] - fRise[2])/fR;
77  w2 = 0.5*fDt*(fPre[5] - fRise[2] - 0.5*w1)/fR;
78  w3 = 0.5*fDt*(fPre[5] - fRise[2] - 0.5*w2)/fR;
79  w4 = 0.5*fDt*(fPre[6] - fRise[2] - w3)/fR;
80  fRise[3] = fRise[2] + (1.0/6.0)*(w1 + 2*w2 + 2*w3 + w4);
81 
82  fRiseDot[1] = (fRise[2] - fRise[0])/fDt;
83  fRiseDot[2] = (fRise[3] - fRise[1])/fDt;
84 
85  double FInv0 = (fPre[0] > 0) ? 1.0/(fFInt + fFSlope*fPre[0]) : 1.0/fFInt;
86  double FInv1 = (fPre[2] > 0) ? 1.0/(fFInt + fFSlope*fPre[2]) : 1.0/fFInt;
87  double FInv2 = (fPre[4] > 0) ? 1.0/(fFInt + fFSlope*fPre[4]) : 1.0/fFInt;
88  w1 = fDt*(fRiseDot[0] - fShape[0]*FInv0);
89  w2 = fDt*(fRiseDot[1] - (fShape[0] + 0.5*w1)*FInv1);
90  w3 = fDt*(fRiseDot[1] - (fShape[0] + 0.5*w2)*FInv1);
91  w4 = fDt*(fRiseDot[2] - (fShape[0] + w3)*FInv2);
92  fShape[1] = fShape[0] + (1.0/6.0)*(w1 + 2*w2 + 2*w3 + w4);
93 
94  return fShape[1];
95  }
96 
97  std::vector<double> PulseShaper::CreateTrace( std::list< std::pair<double,double> >& ADCPulsesIn, std::list< std::pair<double,double> >& SaggedPulses )
98  {
99  ClusterPulses( ADCPulsesIn, SaggedPulses );
100 
101  std::vector<double> trace(fNumClockticksInSpill, 0.0);
103 
104  if (fClusteredPulses.size() == 0) return trace;
105 
106  Reset();
107  int iDigitization(0);
108  double firstPulseTime = fClusteredPulses.front().first;
109  double firstDigitizationTime = fClocktick*fPhase;
110 
111  fTime = (firstPulseTime < firstDigitizationTime) ? firstPulseTime : firstDigitizationTime;
112  fTime -= fClocktick;
113 
114  while (fTime < fnsPerSpill) {
115  double timeToNextInput(9999999.0);
116  if (fClusteredPulses.size() > 0) timeToNextInput = fClusteredPulses.front().first - fTime;
117  double timeToNextDigitization = (iDigitization*fNumClockticksPerDigitization + fPhase)*fClocktick - fTime;
118 
119  double maxStep = fMaxStep;
120  if (std::fabs(fInput[1]) == 0)
121  {
122  double delta = std::fabs(NextRiseDot());
123  if (delta > 0) maxStep = 1.0/delta;
124  //else maxStep = fnsPerDigitization;
125  maxStep = std::min(maxStep, fMaxStep);
126  }
127 
128  //there is neither a clocktick nor an input before we hit the maximum step
129  if ( (maxStep < timeToNextInput) && (maxStep < timeToNextDigitization) ) {
130  Step(0.0, maxStep);
131  }
132 
133  //a clocktick and the input are at about the same time
134  else if ( std::fabs(timeToNextInput - timeToNextDigitization) < 1e-8 ) {
135  double input = fClusteredPulses.front().second;
136  double output = Step(input, timeToNextDigitization);
137  fClusteredPulses.pop_front();
138  if ( iDigitization*fNumClockticksPerDigitization + fPhase < fNumClockticksInSpill) {
139  trace[iDigitization*fNumClockticksPerDigitization + fPhase] += output;
140  }
141  ++iDigitization;
142  }
143 
144  //there is an input before the next clocktick
145  else if (timeToNextInput < timeToNextDigitization) {
146  double input = fClusteredPulses.front().second;
147  Step(input, timeToNextInput);
148  fClusteredPulses.pop_front();
149  }
150 
151  //the clocktick is before the input
152  else {
153  double output = Step(0.0, timeToNextDigitization);
154  if ( iDigitization*fNumClockticksPerDigitization + fPhase < fNumClockticksInSpill) {
155  trace[iDigitization*fNumClockticksPerDigitization + fPhase] += output;
156  }
157  ++iDigitization;
158  }
159  }
160  return trace;
161  }
162 }
163 
ofstream output
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double delta
Definition: runWimpSim.h:98
double NextRiseDot()
Definition: PulseShaper.h:34
std::vector< double > fRise
Definition: PulseShaper.h:53
Common configuration params for SimpleReadout, FPGAAlgorithms, NoiseMaker.
Definition: Cand.cxx:23
std::vector< double > fInput
Definition: PulseShaper.h:51
int fNumClockticksPerDigitization
Definition: IPulseShaper.h:38
std::vector< double > fPre
Definition: PulseShaper.h:52
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: trace.hpp:19
std::vector< double > CreateTrace(std::list< std::pair< double, double > > &ADCPulses, std::list< std::pair< double, double > > &SaggedPulses)
Definition: PulseShaper.cxx:97
std::list< std::pair< double, double > > fClusteredPulses
Definition: IPulseShaper.h:33
double Step(double xNext, double dt)
Definition: PulseShaper.cxx:43
PulseShaper(const fhicl::ParameterSet &pset)
Definition: PulseShaper.cxx:15
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
std::vector< double > fShape
Definition: PulseShaper.h:55
std::vector< double > fRiseDot
Definition: PulseShaper.h:54
void ClusterPulses(std::list< std::pair< double, double > > &ADCPulses, std::list< std::pair< double, double > > &SaggedPulses)