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.
12 // #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
13 #include <cmath>
14 
15 #include <iostream>
16 #include <fstream>
17 
19 : theMCParams(0),
20  theRecoParams(0),
21  theShapes()
22 {
23 /*
24 
25 Reco MC
26 --------------------------------------------------------------------------------------
27 000 not used (reserved)
28 101 101 hpdShape_ HPD (original version)
29 102 102 =101 HPD BV 30 volts in HBP iphi54
30 103 123 hpdShape_v2, hpdShapeMC_v2 HPD (2011. oct version)
31 104 124 hpdBV30Shape_v2, hpdBV30ShapeMC_v2 HPD bv30 in HBP iph54
32 105 125 hpdShape_v2, hpdShapeMC_v2 HPD (2011.11.12 version)
33 201 201 siPMShape_ SiPMs Zecotec shape (HO)
34 202 202 =201, SiPMs Hamamatsu shape (HO)
35 301 301 hfShape_ regular HF PMT shape
36 401 401 regular ZDC shape
37 --------------------------------------------------------------------------------------
38 
39 */
40 
41 
42  float ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3;
43 
44  // HPD Shape Version 1 (used before CMSSW5, until Oct 2011)
45  ts1=8. ; ts2=10. ; ts3=29.3; thpd=4.0; tpre=9.0; wd1=2.0; wd2=0.7; wd3=1.0;
46  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdShape_);
47  theShapes[101] = &hpdShape_;
48  theShapes[102] = theShapes[101];
49 
50  // HPD Shape Version 2 for CMSSW 5. Nov 2011 (RECO and MC separately)
51  ts1=8. ; ts2=10. ; ts3=25.0; thpd=4.0; tpre=9.0; wd1=2.0; wd2=0.7; wd3=1.0;
52  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdShape_v2);
53  theShapes[103] = &hpdShape_v2;
54 
55  ts1=8. ; ts2=10. ; ts3=29.3; thpd=4.0; tpre=7.0; wd1=2.0; wd2=0.7; wd3=1.0;
56  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdShapeMC_v2);
57  theShapes[123] = &hpdShapeMC_v2;
58 
59  // HPD Shape Version 3 for CMSSW 5. Nov 2011 (RECO and MC separately)
60  ts1=8. ; ts2=19. ; ts3=29.3; thpd=4.0; tpre=9.0; wd1=2.0; wd2=0.7; wd3=0.32;
61  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdShape_v3);
62  theShapes[105] = &hpdShape_v3;
63 
64  ts1=8. ; ts2=10. ; ts3=22.3; thpd=4.0; tpre=7.0; wd1=2.0; wd2=0.7; wd3=1.0;
65  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdShapeMC_v3);
66  theShapes[125] = &hpdShapeMC_v3;
67 
68  // HPD with Bias Voltage 30 volts, wider pulse. (HBPlus iphi54)
69 
70  ts1=8. ; ts2=12. ; ts3=31.7; thpd=9.0; tpre=9.0; wd1=2.0; wd2=0.7; wd3=1.0;
71  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdBV30Shape_v2);
72  theShapes[104] = &hpdBV30Shape_v2;
73 
74  ts1=8. ; ts2=12. ; ts3=31.7; thpd=9.0; tpre=9.0; wd1=2.0; wd2=0.7; wd3=1.0;
75  computeHPDShape(ts1,ts2,ts3,thpd,tpre,wd1,wd2,wd3, hpdBV30ShapeMC_v2);
77 
78  // HF and SiPM
79 
82 
83  theShapes[201] = &siPMShape_;
84  theShapes[202] = theShapes[201];
85  theShapes[301] = &hfShape_;
86  //theShapes[401] = new CaloCachedShapeIntegrator(&theZDCShape);
87 
88  /*
89  // backward-compatibility with old scheme
90  theShapes[0] = theShapes[101];
91  //FIXME "special" HB
92  theShapes[1] = theShapes[101];
93  theShapes[2] = theShapes[201];
94  theShapes[3] = theShapes[301];
95  //theShapes[4] = theShapes[401];
96  */
97 }
98 
99 
101  delete theMCParams;
102  delete theRecoParams;
103 }
104 
105 
107 {
109  es.get<HcalMCParamsRcd>().get(p);
110  theMCParams = new HcalMCParams(*p.product());
111 
113  es.get<HcalRecoParamsRcd>().get(q);
114  theRecoParams = new HcalRecoParams(*q.product());
115 
116 // std::cout<<" skdump in HcalPulseShapes::beginRun dupm MCParams "<<std::endl;
117 // std::ofstream skfile("skdumpMCParamsNewFormat.txt");
118 // HcalDbASCIIIO::dumpObject(skfile, (*theMCParams) );
119 }
120 
121 
123 {
124  delete theMCParams;
125  theMCParams = 0;
126 
127  delete theRecoParams;
128  theRecoParams = 0;
129 }
130 
131 
132 //void HcalPulseShapes::computeHPDShape()
133 void HcalPulseShapes::computeHPDShape(float ts1, float ts2, float ts3, float thpd, float tpre,
134  float wd1, float wd2, float wd3, Shape &tmphpdShape_)
135 {
136 
137  /*
138  std::cout << "o HcalPulseShapes::computeHPDShape "
139  << " ts1, ts2, ts3, thpd, tpre, w1, w2, w3 ="
140  << ts1 << ", " << ts2 << ", " << ts3 << ", "
141  << thpd << ", " << tpre << ", " << wd1 << ", " << wd2
142  << ", " << wd3 << std::endl;
143  */
144 
145 // pulse shape time constants in ns
146 /*
147  const float ts1 = 8.; // scintillation time constants : 1,2,3
148  const float ts2 = 10.;
149  const float ts3 = 29.3;
150  const float thpd = 4.; // HPD current collection drift time
151  const float tpre = 9.; // preamp time constant (refit on TB04 data)
152 
153  const float wd1 = 2.; // relative weights of decay exponents
154  const float wd2 = 0.7;
155  const float wd3 = 1.;
156 */
157  // pulse shape componnts over a range of time 0 ns to 255 ns in 1 ns steps
158  unsigned int nbin = 256;
159  tmphpdShape_.setNBin(nbin);
160  std::vector<float> ntmp(nbin,0.0); // zeroing output pulse shape
161  std::vector<float> nth(nbin,0.0); // zeroing HPD drift shape
162  std::vector<float> ntp(nbin,0.0); // zeroing Binkley preamp shape
163  std::vector<float> ntd(nbin,0.0); // zeroing Scintillator decay shape
164 
165  unsigned int i,j,k;
166  float norm;
167 
168  // HPD starts at I and rises to 2I in thpd of time
169  norm=0.0;
170  for(j=0;j<thpd && j<nbin;j++){
171  nth[j] = 1.0 + ((float)j)/thpd;
172  norm += nth[j];
173  }
174  // normalize integrated current to 1.0
175  for(j=0;j<thpd && j<nbin;j++){
176  nth[j] /= norm;
177  }
178 
179  // Binkley shape over 6 time constants
180  norm=0.0;
181  for(j=0;j<6*tpre && j<nbin;j++){
182  ntp[j] = ((float)j)*exp(-((float)(j*j))/(tpre*tpre));
183  norm += ntp[j];
184  }
185  // normalize pulse area to 1.0
186  for(j=0;j<6*tpre && j<nbin;j++){
187  ntp[j] /= norm;
188  }
189 
190 // ignore stochastic variation of photoelectron emission
191 // <...>
192 
193 // effective tile plus wave-length shifter decay time over 4 time constants
194  unsigned int tmax = 6 * (int)ts3;
195 
196  norm=0.0;
197  for(j=0;j<tmax && j<nbin;j++){
198  ntd[j] = wd1 * exp(-((float)j)/ts1) +
199  wd2 * exp(-((float)j)/ts2) +
200  wd3 * exp(-((float)j)/ts3) ;
201  norm += ntd[j];
202  }
203  // normalize pulse area to 1.0
204  for(j=0;j<tmax && j<nbin;j++){
205  ntd[j] /= norm;
206  }
207 
208  unsigned int t1,t2,t3,t4;
209  for(i=0;i<tmax && i<nbin;i++){
210  t1 = i;
211  // t2 = t1 + top*rand;
212  // ignoring jitter from optical path length
213  t2 = t1;
214  for(j=0;j<thpd && j<nbin;j++){
215  t3 = t2 + j;
216  for(k=0;k<4*tpre && k<nbin;k++){ // here "4" is set deliberately,
217  t4 = t3 + k; // as in test fortran toy MC ...
218  if(t4<nbin){
219  unsigned int ntb=t4;
220  ntmp[ntb] += ntd[i]*nth[j]*ntp[k];
221  }
222  }
223  }
224  }
225 
226  // normalize for 1 GeV pulse height
227  norm = 0.;
228  for(i=0;i<nbin;i++){
229  norm += ntmp[i];
230  }
231 
232  //cout << " Convoluted SHAPE ============== " << endl;
233  for(i=0; i<nbin; i++){
234  ntmp[i] /= norm;
235  // cout << " shape " << i << " = " << ntmp[i] << endl;
236  }
237 
238  for(i=0; i<nbin; i++){
239  tmphpdShape_.setShapeBin(i,ntmp[i]);
240  }
241 }
242 
244  // first create pulse shape over a range of time 0 ns to 255 ns in 1 ns steps
245  unsigned int nbin = 256;
246  hfShape_.setNBin(nbin);
247  std::vector<float> ntmp(nbin,0.0); //
248 
249  const float k0=0.7956; // shape parameters
250  const float p2=1.355;
251  const float p4=2.327;
252  const float p1=4.3; // position parameter
253 
254  float norm = 0.0;
255 
256  for(unsigned int j = 0; j < 25 && j < nbin; ++j){
257 
258  float r0 = j-p1;
259  float sigma0 = (r0<0) ? p2 : p2*p4;
260  r0 /= sigma0;
261  if(r0 < k0) ntmp[j] = exp(-0.5*r0*r0);
262  else ntmp[j] = exp(0.5*k0*k0-k0*r0);
263  norm += ntmp[j];
264  }
265  // normalize pulse area to 1.0
266  for(unsigned int j = 0; j < 25 && j < nbin; ++j){
267  ntmp[j] /= norm;
268  hfShape_.setShapeBin(j,ntmp[j]);
269  }
270 }
271 
272 
274 {
275  unsigned int nbin = 512;
276  siPMShape_.setNBin(nbin);
277  std::vector<float> nt(nbin,0.0); //
278 
279  double norm = 0.;
280  for (unsigned int j = 0; j < nbin; ++j) {
281  if (j <= 31.)
282  nt[j] = 2.15*j;
283  else if ((j > 31) && (j <= 96))
284  nt[j] = 102.1 - 1.12*j;
285  else
286  nt[j] = 0.0076*j - 6.4;
287  norm += (nt[j]>0) ? nt[j] : 0.;
288  }
289 
290  for (unsigned int j = 0; j < nbin; ++j) {
291  nt[j] /= norm;
292  siPMShape_.setShapeBin(j,nt[j]);
293  }
294 }
295 
296 
298 HcalPulseShapes::getShape(int shapeType) const
299 {
300 
301  // std::cout << "- HcalPulseShapes::Shape for type "<< shapeType
302  // << std::endl;
303 
304  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
305  if(shapeMapItr == theShapes.end()) {
306  throw cms::Exception("HcalPulseShapes") << "unknown shapeType";
307  return hpdShape_; // should not return this, but...
308  } else {
309  return *(shapeMapItr->second);
310  }
311 }
312 
313 
315 HcalPulseShapes::shape(const HcalDetId & detId) const
316 {
317  if(!theMCParams) {
318  return defaultShape(detId);
319  }
320  int shapeType = theMCParams->getValues(detId)->signalShape();
321  /*
322  int sub = detId.subdet();
323  int depth = detId.depth();
324  int inteta = detId.ieta();
325  int intphi = detId.iphi();
326 
327  std::cout << " HcalPulseShapes::shape cell:"
328  << " sub, ieta, iphi, depth = "
329  << sub << " " << inteta << " " << intphi
330  << " " << depth << " => ShapeId "<< shapeType
331  << std::endl;
332  */
333  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
334  if(shapeMapItr == theShapes.end()) {
335  return defaultShape(detId);
336  } else {
337  return *(shapeMapItr->second);
338  }
339 }
340 
343 {
344  if(!theRecoParams) {
345  return defaultShape(detId);
346  }
347  int shapeType = theRecoParams->getValues(detId.rawId())->pulseShapeID();
348  /*
349  int sub = detId.subdet();
350  int depth = detId.depth();
351  int inteta = detId.ieta();
352  int intphi = detId.iphi();
353 
354  std::cout << ">> HcalPulseShapes::shapeForReco cell:"
355  << " sub, ieta, iphi, depth = "
356  << sub << " " << inteta << " " << intphi
357  << " " << depth << " => ShapeId "<< shapeType
358  << std::endl;
359  */
360 
361  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
362  if(shapeMapItr == theShapes.end()) {
363  return defaultShape(detId);
364  } else {
365  return *(shapeMapItr->second);
366  }
367 }
368 
369 
372 {
373  edm::LogWarning("HcalPulseShapes") << "Cannot find HCAL MC Params ";
374  HcalSubdetector subdet = detId.subdet();
375  switch(subdet) {
376  case HcalBarrel:
377  return hbShape();
378  case HcalEndcap:
379  return heShape();
380  case HcalForward:
381  return hfShape();
382  case HcalOuter:
383  //FIXME doesn't look for SiPMs
384  return hoShape(false);
385  default:
386  throw cms::Exception("HcalPulseShapes") << "unknown detId";
387  break;
388  }
389 }
390 
const Shape & getShape(int shapeType) const
int i
Definition: DBlmapReader.cc:9
void setNBin(int n)
const Shape & heShape() const
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:32
void setShapeBin(int i, float f)
const Shape & hoShape(bool sipm=false) const
void beginRun(edm::EventSetup const &es)
const Shape & shapeForReco(const HcalDetId &detId) const
const Shape & shape(const HcalDetId &detId) const
automatically figures out which shape to return
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
double p4[4]
Definition: TauolaWrapper.h:92
HcalSubdetector
Definition: HcalAssistant.h:32
int j
Definition: DBlmapReader.cc:9
double p2[4]
Definition: TauolaWrapper.h:90
int nt
Definition: AMPTWrapper.h:32
int k[5][pyjets_maxn]
static const double tmax[3]
const Shape & defaultShape(const HcalDetId &detId) const
in case of conditions problems
const Shape & hfShape() const
void computeHPDShape(float, float, float, float, float, float, float, float, Shape &)
const HcalRecoParams * theRecoParams
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
const HcalMCParams * theMCParams
unsigned int signalShape() const
Definition: HcalMCParam.h:38
const Shape & hbShape() const
double p1[4]
Definition: TauolaWrapper.h:89
const Item * getValues(DetId fId) const