CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalPulseShapes.cc
Go to the documentation of this file.
2 #include <cmath>
3 
7 }
8 
9 
11 {
12 
13  // pulse shape time constants in ns
14  const float ts1 = 8.; // scintillation time constants : 1,2,3
15  const float ts2 = 10.;
16  const float ts3 = 29.3;
17  const float thpd = 4.; // HPD current collection drift time
18  const float tpre = 9.; // preamp time constant (refit on TB04 data)
19 
20  const float wd1 = 2.; // relative weights of decay exponents
21  const float wd2 = 0.7;
22  const float wd3 = 1.;
23 
24  // pulse shape componnts over a range of time 0 ns to 255 ns in 1 ns steps
25  int nbin = 256;
26  sh.setNBin(nbin);
27  std::vector<float> ntmp(nbin,0.0); // zeroing output pulse shape
28  std::vector<float> nth(nbin,0.0); // zeroing HPD drift shape
29  std::vector<float> ntp(nbin,0.0); // zeroing Binkley preamp shape
30  std::vector<float> ntd(nbin,0.0); // zeroing Scintillator decay shape
31 
32  int i,j,k;
33  float norm;
34 
35  // HPD starts at I and rises to 2I in thpd of time
36  norm=0.0;
37  for(j=0;j<thpd && j<nbin;j++){
38  nth[j] = 1.0 + ((float)j)/thpd;
39  norm += nth[j];
40  }
41  // normalize integrated current to 1.0
42  for(j=0;j<thpd && j<nbin;j++){
43  nth[j] /= norm;
44  }
45 
46  // Binkley shape over 6 time constants
47  norm=0.0;
48  for(j=0;j<6*tpre && j<nbin;j++){
49  ntp[j] = ((float)j)*exp(-((float)(j*j))/(tpre*tpre));
50  norm += ntp[j];
51  }
52  // normalize pulse area to 1.0
53  for(j=0;j<6*tpre && j<nbin;j++){
54  ntp[j] /= norm;
55  }
56 
57 // ignore stochastic variation of photoelectron emission
58 // <...>
59 
60 // effective tile plus wave-length shifter decay time over 4 time constants
61  int tmax = 6 * (int)ts3;
62 
63  norm=0.0;
64  for(j=0;j<tmax && j<nbin;j++){
65  ntd[j] = wd1 * exp(-((float)j)/ts1) +
66  wd2 * exp(-((float)j)/ts2) +
67  wd3 * exp(-((float)j)/ts3) ;
68  norm += ntd[j];
69  }
70  // normalize pulse area to 1.0
71  for(j=0;j<tmax && j<nbin;j++){
72  ntd[j] /= norm;
73  }
74 
75  int t1,t2,t3,t4;
76  for(i=0;i<tmax && i<nbin;i++){
77  t1 = i;
78  // t2 = t1 + top*rand;
79  // ignoring jitter from optical path length
80  t2 = t1;
81  for(j=0;j<thpd && j<nbin;j++){
82  t3 = t2 + j;
83  for(k=0;k<4*tpre && k<nbin;k++){ // here "4" is set deliberately,
84  t4 = t3 + k; // as in test fortran toy MC ...
85  if(t4<nbin){
86  int ntb=t4;
87  ntmp[ntb] += ntd[i]*nth[j]*ntp[k];
88  }
89  }
90  }
91  }
92 
93  // normalize for 1 GeV pulse height
94  norm = 0.;
95  for(i=0;i<nbin;i++){
96  norm += ntmp[i];
97  }
98 
99  //cout << " Convoluted SHAPE ============== " << endl;
100  for(i=0; i<nbin; i++){
101  ntmp[i] /= norm;
102  // cout << " shape " << i << " = " << ntmp[i] << endl;
103  }
104 
105  for(i=0; i<nbin; i++){
106  sh.setShapeBin(i,ntmp[i]);
107  }
108 }
109 
111 
112  // cout << endl << " ===== computeShapeHF !!! " << endl << endl;
113 
114  const float ts = 3.0; // time constant in t * exp(-(t/ts)**2)
115 
116  // first create pulse shape over a range of time 0 ns to 255 ns in 1 ns steps
117  int nbin = 256;
118  sh.setNBin(nbin);
119  std::vector<float> ntmp(nbin,0.0); //
120 
121  int j;
122  float norm;
123 
124  // HF SHAPE
125  norm = 0.0;
126  for( j = 0; j < 3 * ts && j < nbin; j++){
127  ntmp[j] = ((float)j)*exp(-((float)(j*j))/(ts*ts));
128  norm += ntmp[j];
129  }
130  // normalize pulse area to 1.0
131  for( j = 0; j < 3 * ts && j < nbin; j++){
132  ntmp[j] /= norm;
133 
134  // cout << " nt [" << j << "] = " << ntmp[j] << endl;
135  sh.setShapeBin(j,ntmp[j]);
136  }
137 }
138 
140  nbin_=0;
141  tpeak_=0;
142 }
143 
145  nbin_=n;
146  shape_=std::vector<float>(n,0.0f);
147 }
148 
150  if (i>=0 && i<nbin_) shape_[i]=f;
151 }
152 
154  // shape is in 1 ns steps
155  int i=(int)(t+0.5);
156  float rv=0;
157  if (i>=0 && i<nbin_) rv=shape_[i];
158  return rv;
159 }
160 
161 float HcalPulseShapes::Shape::at(double t) const {
162  // shape is in 1 ns steps
163  int i=(int)(t+0.5);
164  float rv=0;
165  if (i>=0 && i<nbin_) rv=shape_[i];
166  return rv;
167 }
168 
169 float HcalPulseShapes::Shape::integrate(double t1, double t2) const {
170  static const float int_delta_ns = 0.05f;
171  double intval = 0.0;
172 
173  for (double t = t1; t < t2; t+= int_delta_ns) {
174  float loedge = at(t);
175  float hiedge = at(t+int_delta_ns);
176  intval += (loedge+hiedge)*int_delta_ns/2.0;
177  }
178 
179  return (float)intval;
180 }
float operator()(double time) const
int i
Definition: DBlmapReader.cc:9
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
float integrate(double tmin, double tmax) const
void computeHFShape(Shape &s)
tuple norm
Definition: lumiNorm.py:78
int j
Definition: DBlmapReader.cc:9
double f[11][100]
int k[5][pyjets_maxn]
static const double tmax[3]
float at(double time) const
void setShapeBin(int i, float f)
void computeHPDShape(Shape &s)
list at
Definition: asciidump.py:428