ATLAS Offline Software
TileRawChannelBuilderOpt2Filter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 
7 //
8 // Base on the code of Ximo Poveda.
9 // Andrei.Artamonov@cern.ch, July 2008
10 //
11 // TileRawChannelBuilderOpt2Filter.cxx
12 //
13 // implementation of the Optimal Filtering based on Lagrange multipliers
14 // for energy/time reconstruction in TileCal
15 //
17 
18 // Tile includes
25 #include "TileEvent/TileDigits.h"
26 #include "CaloIdentifier/TileID.h"
29 
30 // Atlas includes
31 #include "AthAllocators/DataPool.h"
33 
34 // Gaudi includes
35 #include "Gaudi/Property.h"
36 
37 #include "CLHEP/Matrix/Matrix.h"
38 //using namespace std;
39 #include <algorithm>
40 
41 //interface stuff
42 static const InterfaceID IID_ITileRawChannelBuilderOpt2Filter("TileRawChannelBuilderOpt2Filter", 1, 0);
43 
44 
46  return IID_ITileRawChannelBuilderOpt2Filter;
47 }
48 
49 
50 #define TILE_Opt2FilterBUILDERVERBOSE false
51 
53  const std::string& name, const IInterface *parent)
55  , m_nSignal(0)
56  , m_nNegative(0)
57  , m_nCenter(0)
58  , m_nConst(0)
59  , m_nSamples(0)
60  , m_t0SamplePosition(0)
61 {
62  //declare interfaces
63  declareInterface< TileRawChannelBuilder >( this );
64  declareInterface< TileRawChannelBuilderOpt2Filter >(this);
65 
66  m_rawChannelContainerKey = "TileRawChannelOpt2";
67 
68  //declare properties
69  declareProperty("MaxIterations", m_maxIterations = 5);
70  declareProperty("PedestalMode", m_pedestalMode = 17);
71  declareProperty("TimeForConvergence", m_timeForConvergence = 0.5);
72  declareProperty("ConfTB", m_confTB = false);
73  declareProperty("OF2", m_of2 = true);
74  declareProperty("Minus1Iteration", m_minus1Iter = false);
75  declareProperty("AmplitudeCorrection", m_correctAmplitude = false);
76  declareProperty("TimeCorrection", m_correctTimeNI = false);
77  declareProperty("BestPhase", m_bestPhase = false);
78  declareProperty("EmulateDSP", m_emulateDsp = false);
79  declareProperty("NoiseThresholdHG", m_noiseThresholdHG = 5);
80  declareProperty("NoiseThresholdLG", m_noiseThresholdLG = 3);
81  declareProperty("MinTime", m_minTime = 0.0);
82  declareProperty("MaxTime", m_maxTime = -1.0);
83 }
84 
85 
87 }
88 
89 
91 
92  ATH_MSG_INFO( "initialize()" );
93 
94  m_rChType = TileFragHash::OptFilterOffline; // type for offline Opt Filter
95 
96  // init in superclass
98 
99  if (m_maxIterations != 1) m_correctTimeNI = false;
100 
101  // bits 12-15 - various options
102  // if (m_correctTimeNI) m_bsflags |= 0x1000;
103  if (m_correctAmplitude) m_bsflags |= 0x2000;
104  if (m_maxIterations > 1) m_bsflags |= 0x4000;
105  if (m_bestPhase) m_bsflags |= 0x8000;
106 
107  ATH_MSG_DEBUG( " MaxIterations=" << m_maxIterations
108  << " PedestalMode=" << m_pedestalMode
109  << " TimeForConvergence=" << m_timeForConvergence
110  << " ConfTB=" << m_confTB
111  << " OF2=" << m_of2
112  << " Minus1Iteration=" << m_minus1Iter
113  << " AmplitudeCorrection=" << m_correctAmplitude
114  << " TimeCorrection=" << m_correctTimeNI
115  << " Best Phase " << m_bestPhase );
116 
117  ATH_MSG_DEBUG( " NoiseThresholdHG=" << m_noiseThresholdHG
118  << " NoiseThresholdLG=" << m_noiseThresholdLG);
119 
122  if (m_maxTime < m_minTime) { // set time window if it was not set from jobOptions
123  m_maxTime = 25 * (m_nSamples - m_t0SamplePosition - 1);
125  }
126  ATH_MSG_DEBUG(" NSamples=" << m_nSamples
127  << " T0Sample=" << m_t0SamplePosition
128  << " minTime=" << m_minTime
129  << " maxTime=" << m_maxTime );
130 
131  if (m_pedestalMode % 10 > 2 && m_nSamples != m_pedestalMode % 10) {
132  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Changing PedestalMode from " << m_pedestalMode;
133  m_pedestalMode = (m_pedestalMode / 10) * 10 + m_nSamples;
134  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " to " << m_pedestalMode << endmsg;
135  }
136 
137 
138  if (m_nSamples != 7 && (m_pedestalMode == 71 || m_pedestalMode == 7621)) {
139  ATH_MSG_ERROR( "Incompatable pedestal mode [" << m_pedestalMode
140  << "] and number of samples [" << m_nSamples << "]" );
141  return StatusCode::FAILURE;
142  }
143 
144  m_nSignal = 0;
145  m_nNegative = 0;
146  m_nCenter = 0;
147  m_nConst = 0;
148 
149  //=== get TileCondToolOfc
150  ATH_CHECK( m_tileCondToolOfc.retrieve() );
151 
152  //=== get TileCondToolNoiseSample
153  ATH_CHECK( m_tileToolNoiseSample.retrieve(EnableTool{m_pedestalMode == -1}) );
154 
155  if (m_bestPhase) {
156  //=== get TileToolTiming
157  // TileToolTiming can be disabled in the TileRawChannelBuilder
158  if (!m_tileToolTiming.isEnabled()) {
159  m_tileToolTiming.enable();
160  }
161  ATH_CHECK( m_tileToolTiming.retrieve() );
162  }
163 
164  ATH_MSG_INFO( "initialization completed" );
165 
166  return StatusCode::SUCCESS;
167 }
168 
169 
171 
172  if (msgLvl(MSG::VERBOSE)) {
173  if (m_maxIterations == 1) { // Without iterations
174  msg(MSG::VERBOSE) << "Counters: Signal=" << m_nSignal
175  << " Constant=" << m_nConst
176  << " Total=" << m_nSignal + m_nConst << endmsg;
177  } else {
178  msg(MSG::VERBOSE) << "Counters: Signal=" << m_nSignal
179  << " Negat=" << m_nNegative
180  << " Center=" << m_nCenter
181  << " Constant=" << m_nConst
182  << " Total=" << m_nSignal + m_nNegative + m_nConst + m_nCenter << endmsg;
183  }
184  }
185 
186  ATH_MSG_DEBUG( "Finalizing" );
187 
188  return StatusCode::SUCCESS;
189 }
190 
191 
192 TileRawChannel * TileRawChannelBuilderOpt2Filter::rawChannel(const TileDigits* digits, const EventContext& ctx) {
193 
194  ++m_chCounter;
195 
196  double pedestal = 0.;
197  double energy = 0.;
198  double time = 0.;
199  double chi2 = 0.;
200  m_digits = digits->samples();
201  m_digits.erase(m_digits.begin(),m_digits.begin()+m_firstSample);
202  m_digits.resize(m_nSamples);
203  const HWIdentifier adcId = digits->adc_HWID();
204  int gain = m_tileHWID->adc(adcId);
205 
206  ATH_MSG_VERBOSE( "Building Raw Channel, with OptFilter, HWID:" << m_tileHWID->to_string(adcId)
207  << " gain=" << gain );
208 
209  int ros = m_tileHWID->ros(adcId);
210  int drawer = m_tileHWID->drawer(adcId);
211  int channel = m_tileHWID->channel(adcId);
212 
213  chi2 = filter(ros, drawer, channel, gain, pedestal, energy, time, ctx);
214 
215  unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer);
216 
217  if (m_calibrateEnergy) {
218  energy = m_tileToolEmscale->doCalibCis(drawerIdx, channel, gain, energy);
219  }
220 
221  if (msgLvl(MSG::VERBOSE)) {
222  msg(MSG::VERBOSE) << "Creating OptFilter RawChannel"
223  << " a=" << energy
224  << " t=" << time
225  << " ped=" << pedestal
226  << " q=" << chi2 << endmsg;
227 
228  msg(MSG::VERBOSE) << "digits:";
229 
230  for (unsigned int i = 0; i < m_digits.size(); ++i)
231  msg(MSG::VERBOSE) << " " << m_digits[i];
232 
233  msg(MSG::VERBOSE) << " " << endmsg;
234  }
235 
236  // return new TileRawChannel
237  // TileRawChannel *rawCh = new TileRawChannel(adcId,OptFilterEne,OptFilterTime,OptFilterChi2,OptFilterPed);
239  TileRawChannel *rawCh = tileRchPool.nextElementPtr();
240  rawCh->assign (adcId,
241  energy,
242  time,
243  chi2,
244  pedestal);
245 
246  if (m_correctTime
247  && (time != 0
248  && time < m_maxTime
249  && time > m_minTime)) {
250 
251  time -= m_tileToolTiming->getSignalPhase(drawerIdx, channel, gain);
252  rawCh->insertTime(time);
253 
254  ATH_MSG_VERBOSE( "Correcting time, new time=" << rawCh->time() );
255 
256  }
257 
258  if (TileID::HIGHGAIN == gain) {
259  ++m_nChH;
260  m_RChSumH += energy;
261  } else {
262  ++m_nChL;
263  m_RChSumL += energy;
264  }
265 
266  return rawCh;
267 }
268 
269 
271 
272  ATH_MSG_VERBOSE( " findMaxDigitPosition()" );
273 
274  int iMaxDigit = 0;
275  float maxDigit = 0.;
276  bool saturated = false;
277 
278  for (unsigned int i = 0; i < m_digits.size(); i++) {
279  if (m_digits[i] > m_ADCmaxMinusEps) saturated = true;
280  if (maxDigit < m_digits[i]) {
281  maxDigit = m_digits[i];
282  iMaxDigit = i;
283  }
284  }
285 
286  if (msgLvl(MSG::VERBOSE)) {
287  for (unsigned int i = 0; i < m_digits.size(); i++) {
288  msg(MSG::VERBOSE) << " " << m_digits[i];
289  }
290 
291  msg(MSG::VERBOSE) << "; Max: digit[" << iMaxDigit << "]=" << maxDigit << endmsg;
292 
293  if (saturated) msg(MSG::VERBOSE) << " Samples saturated" << endmsg;
294  }
295 
296  return iMaxDigit;
297 }
298 
299 
300 float TileRawChannelBuilderOpt2Filter::getPedestal(int ros, int drawer, int channel, int gain, const EventContext &ctx) {
301  float pedestal = 0.;
302 
303  switch (m_pedestalMode) {
304  case -1:
305  // use pedestal from conditions DB
307  break;
308  case 7:
309  pedestal = m_digits[6];
310  break;
311  case 9:
312  pedestal = m_digits[8];
313  break;
314  case 12:
315  pedestal = .5 * (m_digits[0] + m_digits[1]);
316  break;
317  case 17:
318  pedestal = .5 * (m_digits[0] + m_digits[6]);
319  break;
320  case 19:
321  pedestal = .5 * (m_digits[0] + m_digits[8]);
322  break;
323  case 71:
324  pedestal = std::min(m_digits[0], m_digits[6]);
325  break;
326  case 7621:
327  pedestal = 0.5 * std::min(m_digits[0] + m_digits[1], m_digits[5] + m_digits[6]);
328  break;
329  default:
330  pedestal = m_digits[0];
331  break;
332  }
333 
334  ATH_MSG_VERBOSE("getPedestal(): pedestal=" << pedestal);
335 
336  return pedestal;
337 }
338 
339 
341  , int &gain, double &pedestal, double &amplitude, double &time, const EventContext &ctx) {
342 
343  ATH_MSG_VERBOSE( "filter()" );
344 
345  amplitude = 0.;
346  time = 0.;
347  double chi2 = 0.;
348 
349  auto minMaxDigits = std::minmax_element(m_digits.begin(), m_digits.end());
350  float minDigit = *minMaxDigits.first;
351  float maxDigit = *minMaxDigits.second;
352 
353  if (maxDigit - minDigit < 0.01) { // constant value in all samples
354 
355  pedestal = minDigit;
356  chi2 = 0.;
357  ATH_MSG_VERBOSE( "CASE NO SIGNAL: maxdig-mindig = " << maxDigit << "-" << minDigit
358  << " = " << maxDigit - minDigit );
359 
360  m_nConst++;
361 
362  } else {
363 
364  pedestal = getPedestal(ros, drawer, channel, gain, ctx);
365  double phase = 0.;
366  int nIterations = 0;
367 
368  if (m_maxIterations == 1) { // Without iterations
369  unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer);
370  // AA 3.10.08 --- take best phase from COOL
371  if (m_bestPhase) {
372  // AS 19.11.09 - note minus sign here - time in DB is opposite to best phase
373  phase = -m_tileToolTiming->getSignalPhase(drawerIdx, channel, gain);
374  ATH_MSG_VERBOSE( "Best phase: " << phase
375  << " drawerIdx " << drawerIdx
376  << " channel " << channel );
377  }
378 
379  double ofcPhase(phase);
380  chi2 = compute(ros, drawer, channel, gain, pedestal, amplitude, time, ofcPhase, ctx);
381 
382  // If weights for tau=0 are used, deviations are seen in the amplitude =>
383  // function to correct the amplitude
385  && amplitude > m_ampMinThresh
386  && time > m_timeMinThresh
387  && time < m_timeMaxThresh) {
388 
389  amplitude *= correctAmp(time, m_of2);
390  ATH_MSG_VERBOSE( "Amplitude corrected by " << correctAmp(time, m_of2)
391  << " new amplitude is " << amplitude );
392  }
393 
394  if (m_correctTimeNI ) {
395  ATH_MSG_VERBOSE( "Time " << time
396  << " is corrected by " << correctTime(time, m_of2)
397  << ", new time is " << time + correctTime(time, m_of2));
398 
399  time += correctTime(time, m_of2);
400  }
401 
402  time += (phase - ofcPhase); // correct time if actual phase used in the calculation is different from required
403 
404  if (time > m_maxTime) time = m_maxTime;
405  if (time < m_minTime) time = m_minTime;
406 
407  m_nSignal++;
408 
409  } else { // With iterations => 3 cases defined for correct pedestal treatment
410 
411  // values used for pedestal-like events when iterations are performed
412  //const int sigma_hi = 5;
413  //const int sigma_lo = 3;
415  int digits_size_1 = m_digits.size() - 1;
416 
417  // Signal events: OF with iterations
418  if ((maxDigit - m_digits[0] > sigma)
419  || (m_digits[0] - m_digits[digits_size_1] > 4 * sigma)) {
420 
421  ATH_MSG_VERBOSE( "CASE Signal: maxdig-OptFilterDigits[0]="
422  << maxDigit - m_digits[0]
423  << " OptFilterDigits[0]-OptFilterDigits[ digits_size_1]="
424  << m_digits[0] - m_digits[digits_size_1] );
425 
426  nIterations = iterate(ros, drawer, channel, gain, pedestal, amplitude, time, chi2, ctx);
427 
428  ATH_MSG_VERBOSE( "number of iterations= " << nIterations );
429 
430  m_nSignal++;
431  } else if (m_digits[0] - minDigit > sigma) { //Pedestal events: OF with negative iterations
432 
433  if (msgLvl(MSG::VERBOSE)) {
434  msg(MSG::VERBOSE) << "CASE NEGATIVE: maxdig-OptFilterDigits[0]="
435  << maxDigit - m_digits[0]
436  << " OptFilterDigits[0]-OptFilterDigits[digits_size_1]="
437  << m_digits[0] - m_digits[digits_size_1] << endmsg;
438  msg(MSG::VERBOSE) << "OptFilterDigits[0]-mindig="
439  << m_digits[0] - minDigit << endmsg;
440  }
441 
442  for (int i = 0; i <= digits_size_1; i++) // Mirror around pedestal
443  m_digits[i] = pedestal - (m_digits[i] - pedestal);
444 
445  nIterations = iterate(ros, drawer, channel, gain, pedestal, amplitude, time, chi2, ctx);
446 
447  amplitude = -amplitude;
448 
449  ++m_nNegative;
450 
451  } else { // Gaussian center: no iterations and phase=0
452 
453  if (msgLvl(MSG::VERBOSE)) {
454  msg(MSG::VERBOSE) << "CASE CENTER: maxdig-OptFilterDigits[0]="
455  << maxDigit - m_digits[0]
456  << " OptFilterDigits[0]-OptFilterDigits[digits_size_1]="
457  << m_digits[0] - m_digits[digits_size_1] << endmsg;
458  msg(MSG::VERBOSE) << "OptFilterDigits[0]-mindig="
459  << m_digits[0] - minDigit
460  << endmsg;
461  }
462 
463  if (m_bestPhase) {
464  unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer);
465  // AS 19.11.09 - note minus sign here - time in DB is opposite to best phase
466  phase = -m_tileToolTiming->getSignalPhase(drawerIdx, channel, gain);
467  ATH_MSG_VERBOSE( "Best phase: " << phase
468  << " drawerIdx " << drawerIdx
469  << " channel " << channel );
470  }
471 
472  chi2 = compute(ros, drawer, channel, gain, pedestal, amplitude, time, phase, ctx);
473 
474  // If weights for tau=0 are used, deviations are seen in the amplitude =>
475  // function to correct the amplitude
477  && amplitude > m_ampMinThresh
478  && time > m_timeMinThresh
479  && time < m_timeMaxThresh) {
480 
481  amplitude *= correctAmp(time, m_of2);
482  ATH_MSG_VERBOSE( "Amplitude corrected by " << correctAmp(time, m_of2)
483  << " new amplitude is " << amplitude );
484  }
485 
486  if (m_bestPhase) {
487  time = -phase;
488  chi2 = -chi2;
489  } else {
490  time = 0.;
491  }
492 
493  m_nCenter++;
494 
495  }
496 
497  }
498 
499  }
500 
501  return chi2;
502 }
503 
504 
506  double &pedestal, double &amplitude, double &time, double &chi2, const EventContext &ctx) {
507 
508  ATH_MSG_VERBOSE( "iterate()" );
509 
510  int nIterations = 0;
511  double savePhase = 0.0;
512  double phase = 0.0;
513  time = -1000.;
514 
515  // Mythic -1 iteration or DSP emulation case
516  if (m_minus1Iter || (m_emulateDsp && (m_maxIterations > 1)))
518 
519  while ((time > m_timeForConvergence
520  || time < (-1.) * m_timeForConvergence
521  || m_emulateDsp)
522  && nIterations < m_maxIterations) {
523 
524  chi2 = compute(ros, drawer, channel, gain, pedestal, amplitude, time, phase, ctx);
525 
526  savePhase = phase;
527 
528  if (m_emulateDsp)
529  phase -= round(time); // rounding phase to integer like in DSP
530  else
531  phase -= time; // no rounding at all for OFC on the fly
532 
533  if (phase > m_maxTime) phase = m_maxTime;
534  if (phase < m_minTime) phase = m_minTime;
535 
536  ++nIterations;
537  ATH_MSG_VERBOSE( " OptFilter computed with phase=" << savePhase
538  << " Time=" << time
539  << " END ITER=" << nIterations
540  << " new phase=" << phase
541  << " chi2=" << chi2
542  << " Amp=" << amplitude );
543  }
544 
545  time -= savePhase;
546  if (time > m_maxTime) time = m_maxTime;
547  if (time < m_minTime) time = m_minTime;
548 
549  ATH_MSG_VERBOSE( "OptFilterEne=" << amplitude
550  << " Phase=" << savePhase
551  << " Absolute Time=" << time );
552 
553  return nIterations;
554 }
555 
556 
558  double &pedestal, double &amplitude, double &time, double& phase, const EventContext &ctx) {
559 
560  ATH_MSG_VERBOSE( "compute();"
561  << " ros=" << ros
562  << " drawer=" << drawer
563  << " channel=" << channel
564  << " gain=" << gain );
565 
566  int i = 0, digits_size = m_digits.size();
567  double chi2 = 0.;
568  double a[99];
569  double b[99];
570  double c[99];
571  double g[99];
572  double dg[99];
573 
574  amplitude = 0.;
575  time = 0.;
576  float ofcPhase = (float) phase;
577  unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer);
578  TileOfcWeightsStruct weights;
579  if (m_tileCondToolOfc->getOfcWeights(drawerIdx, channel, gain, ofcPhase, m_of2, weights, ctx).isFailure())
580  {
581  ATH_MSG_ERROR( "getOfcWeights fails" );
582  return 0;
583  }
584 
585  phase = ofcPhase;
586 
587  for (i = 0; i < digits_size; ++i) {
588  a[i] = weights.w_a[i];
589  b[i] = weights.w_b[i];
590  g[i] = weights.g[i];
591  dg[i] = weights.dg[i];
592  if (m_of2) c[i] = weights.w_c[i]; // [OptFilterPha+100];
593  }
594 
595  // for DSP emulation
596  short a_int[99], ascale = 0, calib = 0, calib_offset = 0;
597  short b_int[99], bscale = 0;
598  unsigned short scale = 0, round = 0, lut = 0;
599  int slope_scale = 0;
600  int dspEnergy = 0, dspTime = 0;
601  if (m_emulateDsp) {
602  dspEnergy = 0;
603  dspTime = 0;
604  slope_scale = (int) truncf(log(std::pow(2., 15) - 1.) / log(2.));
605  calib_offset = (short) roundf(std::pow(2., slope_scale));
606  ofc2int(digits_size, a, a_int, ascale);
607  ofc2int(digits_size, b, b_int, bscale);
608  calib = ascale + slope_scale;
609  }
610 
611  ATH_MSG_VERBOSE( "OptFilterPha=" << phase );
612 
613  if (m_of2) {
614  if (m_emulateDsp) {
615  pedestal = m_digits[0];
616  for (i = 0; i < digits_size; ++i) {
617  dspEnergy += a_int[i] * ((int) m_digits[i]);
618  dspTime += b_int[i] * ((int) m_digits[i]);
619  amplitude += a[i] * m_digits[i]; //aa temp
620  time += b[i] * m_digits[i]; //aa temp
621  }
622  } else {
623  pedestal = 0.;
624  for (i = 0; i < digits_size; ++i) {
625  amplitude += a[i] * m_digits[i];
626  time += b[i] * m_digits[i];
627  pedestal += c[i] * m_digits[i];
628  }
629  }
630  } else {
631  if (m_emulateDsp) pedestal = m_digits[0];
632  for (i = 0; i < digits_size; ++i) {
633  amplitude += a[i] * (m_digits[i] - pedestal);
634  time += b[i] * (m_digits[i] - pedestal);
635  }
636  }
637 
638  if (m_emulateDsp) {
640  round = 1 << (ascale - 2);
641  short e2Dsp = (unsigned short) ((dspEnergy + round) >> (ascale - 1));
642  size_t OptInd = abs(e2Dsp);
643  if (OptInd >= (sizeof(lookup) / sizeof(short))) OptInd = 0;
644  lut = lookup[OptInd];
645  scale = bscale - 4 + lookup[0] - 9;
646  round = 1 << (scale - 1);
647  // int told=dspTime;
648  dspTime = (((dspTime + 0x100) >> 9) * lut + round) >> scale;
649  // printf(" 1 OptFilterTime %f OptFilterTimeDSP %d e2 %d round %d lut %d bscale %d scale %d told %d\n",OptFilterTime/OptFilterEne,OptFilterTimeDSP,OptFilterE2DSP,round,lut,bscale,scale,told);
650  time = dspTime / 16.0;
651 
652  // printf(" 1 OptFilterEneDSP %d calib_offset %d calib %d \n",OptFilterEneDSP,calib_offset,calib);
653 
654  dspEnergy = (dspEnergy + m_i_ADCmaxPlus1) >> 11;
655  // printf(" 2 OptFilterEneDSP %d \n",OptFilterEneDSP);
656  dspEnergy = (dspEnergy * calib_offset + (1 << (calib - 15 - 1))) >> (calib - 15);
657  // printf(" 3 OptFilterEneDSP %d \n",OptFilterEneDSP);
658  double goffset = (gain == 0) ? 512 : 2048;
659  dspEnergy = (int) (dspEnergy + goffset);
660  amplitude = (dspEnergy - goffset) / 16.;
661  }
662 
663  bool goodEnergy = (fabs(amplitude) > 1.0e-04);
664  if (goodEnergy) {
665  if (msgLvl(MSG::VERBOSE)) {
666  msg(MSG::VERBOSE) << "OptFilterEne=" << amplitude << endmsg;
667  msg(MSG::VERBOSE) << "OptFilterTime*OptFilterEne=" << time << endmsg;
668  }
669  if (!m_emulateDsp) time /= amplitude;
670  } else {
671  if (msgLvl(MSG::VERBOSE)) {
672  msg(MSG::VERBOSE) << "OptFilterEne=" << amplitude
673  << " ... assuming 0.0" << endmsg;
674  msg(MSG::VERBOSE) << "OptFilterTime*OptFilterEne=" << time
675  << " ... assuming 0.0" << endmsg;
676  }
677  time = amplitude = 0.0;
678  }
679 
680  /* new QF in both cases now 02.2010 AA
681  if(!m_emulatedsp) {
682  for (i=0; i<digits_size; ++i){
683  if(OptFilterDigits[i]>0)
684  OptFilterChi2 += 50*fabs(OptFilterDigits[i]-(OptFilterEne*g[i]+OptFilterPed))/OptFilterDigits[i];
685  else
686  OptFilterChi2 += 50*fabs(OptFilterDigits[i]-(OptFilterEne*g[i]+OptFilterPed))/m_f_ADCmaxPlus1;
687  }
688  } else {
689  */
690  for (i = 0; i < digits_size; ++i) {
691  double dqf = m_digits[i] - amplitude * g[i] + amplitude * time * dg[i] - pedestal;
692  chi2 += dqf * dqf;
693  }
694  chi2 = sqrt(chi2);
695  // new QF }
696 
697  // std::cout << "emulate " << m_emulatedsp << " OptFilterEne " << OptFilterEne << " OptFilterDigits[3]" << OptFilterDigits[3] << " OptFilterTime="<<OptFilterTime<<" OptFilterPed="<<OptFilterPed<<" OptFilterChi2="<<OptFilterChi2<<" g 3 " << g[3] << " dg 1 3 5 " << dg[1] << " " << dg[3] << " " << dg[5] <<std::endl;
698  if (fabs(chi2) > 1.0e-04 || goodEnergy) {
699  if (msgLvl(MSG::VERBOSE)) {
700  msg(MSG::VERBOSE) << "OptFilterTime=" << time << endmsg;
701  msg(MSG::VERBOSE) << "OptFilterPed=" << pedestal << endmsg;
702  msg(MSG::VERBOSE) << "OptFilterChi2=" << chi2 << endmsg;
703  }
704  } else {
705  if (msgLvl(MSG::VERBOSE)) {
706  msg(MSG::VERBOSE) << "OptFilterTime=" << time << endmsg;
707  msg(MSG::VERBOSE) << "OptFilterPed=" << pedestal << endmsg;
708  msg(MSG::VERBOSE) << "OptFilterChi2=" << chi2
709  << " ... assuming 0.0" << endmsg;
710  }
711  chi2 = 0.0;
712  }
713  return chi2;
714 }
715 
716 
717 void TileRawChannelBuilderOpt2Filter::ofc2int(int nDigits, double* w_off, short* w_int, short &scale) {
718  // Number of bits of the integer word (signed -1 )
719  int numberBits = 16;
720  numberBits = numberBits - 1;
721 
722  // Get Absolute Maximum
723  double max = -10000.;
724  double sum = 0.;
725  for (int i = 0; i < nDigits; i++) {
726  sum += w_off[i];
727  if (fabs(w_off[i]) > max) max = fabs(w_off[i]);
728  // printf("idig = %d weight = %f \n",i, w_off[i]);
729  }
730  if (fabs(sum) > max) max = fabs(sum);
731 
732  // Get Scale at Maximum
733  scale = 0;
734  if (max != 0) scale = (int) truncf(log((std::pow(2., numberBits) - 1.) / max) / log(2.));
735 
736  // Convert to Integer the weights and the sum
737  for (int i = 0; i < nDigits; i++)
738  w_int[i] = (short) roundf(w_off[i] * std::pow(2., scale));
739  //aa w_sum_dsp = (short) roundf(sum*std::pow(2.,scale));
740 
741 
742 /*
743  if (msgLvl(MSG::VERBOSE)) {
744  printf("\nAbsolute Max value = %15.8f -> Scale = %3d\n",max,scale);
745 
746  for (int i=0; i<ndigits; i++)
747  {
748  if ( i == 0){
749  printf("\n Offline Off*Scale Dsp/scale Dsp Scale \n");
750  printf("----------------------------------------------------------------------------------\n");
751  printf(" %17.10f %17.10f %17.10f 0x%04X %3d \n",
752  w_off[i],w_off[i]*std::pow(2.,scale), w_int[i]/std::pow(2.,scale),(unsigned int) w_int[i] ,scale);
753  } else {
754  printf(" %17.10f %17.10f %17.10f 0x%04X \n",
755  w_off[i],w_off[i]*std::pow(2.,scale), w_int[i]/std::pow(2.,scale),(unsigned int)w_int[i]);
756  }
757  }
758  //aa printf(" %17.10f %17.10f %17.10f 0x%04X <- SUM\n",
759  //aa sum,sum*std::pow(2.,scale),w_sum_dsp/std::pow(2.,scale),(unsigned int)w_sum_dsp);
760  }
761 */
762  return;
763 }
764 
TileRawChannelBuilderOpt2Filter::m_minTime
double m_minTime
min allowed time = -25*(m_nSamples-1)/2
Definition: TileRawChannelBuilderOpt2Filter.h:106
ReadOfcFromCool.phase
phase
Definition: ReadOfcFromCool.py:127
TileRawChannelBuilderOpt2Filter::TileRawChannelBuilderOpt2Filter
TileRawChannelBuilderOpt2Filter(const std::string &type, const std::string &name, const IInterface *parent)
Constructor.
Definition: TileRawChannelBuilderOpt2Filter.cxx:52
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
xAOD::short
short
Definition: Vertex_v1.cxx:165
TileRawChannelBuilder::m_dataPoollSize
int m_dataPoollSize
Definition: TileRawChannelBuilder.h:204
TileRawChannel.h
max
#define max(a, b)
Definition: cfImp.cxx:41
checkCoolLatestUpdate.dg
dg
Definition: checkCoolLatestUpdate.py:9
TileRawChannelBuilderOpt2Filter::m_noiseThresholdLG
int m_noiseThresholdLG
Definition: TileRawChannelBuilderOpt2Filter.h:112
plotting.yearwise_efficiency.channel
channel
Definition: yearwise_efficiency.py:28
TileRawChannelBuilderOpt2Filter::m_correctTimeNI
bool m_correctTimeNI
If true, resulting time is corrected when using method without iteration.
Definition: TileRawChannelBuilderOpt2Filter.h:94
TileRawChannelBuilder::m_tileToolEmscale
ToolHandle< TileCondToolEmscale > m_tileToolEmscale
Definition: TileRawChannelBuilder.h:166
TileRawChannelBuilderOpt2::lookup
const unsigned short lookup[2401]
Definition: TileRawChannelBuilderOpt2FilterLookup.h:21
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
TileRawChannelBuilderOpt2Filter::compute
double compute(int ros, int drawer, int channel, int gain, double &pedestal, double &amplitude, double &time, double &phase, const EventContext &ctx)
Definition: TileRawChannelBuilderOpt2Filter.cxx:557
TileRawChannelBuilder::m_chCounter
unsigned int m_chCounter
Definition: TileRawChannelBuilder.h:196
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TileOfcWeightsStruct::g
float g[99]
Definition: ITileCondToolOfc.h:14
TileRawChannelBuilderOpt2FilterLookup.h
TileInfo::NdigitSamples
int NdigitSamples() const
Returns the number of sammples (digits) per event.
Definition: TileInfo.h:75
TileRawChannelBuilder::m_ampMinThresh
float m_ampMinThresh
correct amplitude if it's above amplitude threshold (in ADC counts)
Definition: TileRawChannelBuilder.h:152
TileRawChannelBuilder::m_calibrateEnergy
bool m_calibrateEnergy
Definition: TileRawChannelBuilder.h:142
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
TileOfcWeightsStruct::w_c
float w_c[99]
Definition: ITileCondToolOfc.h:14
TileRawChannelBuilderOpt2Filter::m_nSignal
int m_nSignal
internal counters
Definition: TileRawChannelBuilderOpt2Filter.h:98
TileRawChannelBuilder::m_tileHWID
const TileHWID * m_tileHWID
Definition: TileRawChannelBuilder.h:161
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
CaloCondBlobAlgs_fillNoiseFromASCII.gain
gain
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:110
MuonGM::round
float round(const float toRound, const unsigned int decimals)
Definition: Mdt.cxx:27
TileRawChannelBuilderOpt2Filter::m_pedestalMode
int m_pedestalMode
pedestal mode to use
Definition: TileRawChannelBuilderOpt2Filter.h:88
TileRawChannelBuilderOpt2Filter.h
TileRawChannelBuilder::initialize
virtual StatusCode initialize()
Initializer.
Definition: TileRawChannelBuilder.cxx:98
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
TileRawChannelBuilderOpt2Filter::initialize
virtual StatusCode initialize()
Initialize method.
Definition: TileRawChannelBuilderOpt2Filter.cxx:90
TileInfo.h
Tile_Base_ID::HIGHGAIN
@ HIGHGAIN
Definition: Tile_Base_ID.h:57
TileCalibUtils.h
TileRawChannelBuilder::m_ADCmaxMinusEps
float m_ADCmaxMinusEps
Definition: TileRawChannelBuilder.h:221
TileRawChannelBuilderOpt2Filter::m_tileCondToolOfc
ToolHandle< ITileCondToolOfc > m_tileCondToolOfc
Definition: TileRawChannelBuilderOpt2Filter.h:69
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
TileRawData::adc_HWID
HWIdentifier adc_HWID(void) const
Definition: TileRawData.h:53
HWIdentifier
Definition: HWIdentifier.h:13
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
TileRawChannelBuilderOpt2Filter::m_minus1Iter
bool m_minus1Iter
bool variable for whether to apply -1 iteration (initial phase guess)
Definition: TileRawChannelBuilderOpt2Filter.h:92
TileRawChannelBuilder::m_tileToolTiming
ToolHandle< TileCondToolTiming > m_tileToolTiming
Definition: TileRawChannelBuilder.h:169
TileRawChannel::time
float time(int ind=0) const
Definition: TileRawChannel.h:103
Example_ReadSampleNoise.drawer
drawer
Definition: Example_ReadSampleNoise.py:39
TileRawChannelBuilderOpt2Filter::interfaceID
static const InterfaceID & interfaceID()
AlgTool InterfaceID.
Definition: TileRawChannelBuilderOpt2Filter.cxx:45
TileHWID::channel
int channel(const HWIdentifier &id) const
extract channel field from HW identifier
Definition: TileHWID.h:189
TileID.h
TileRawChannelBuilderOpt2Filter::m_timeForConvergence
double m_timeForConvergence
minimum time difference to quit iteration procedure
Definition: TileRawChannelBuilderOpt2Filter.h:90
TileRawChannelBuilderOpt2Filter::m_tileToolNoiseSample
ToolHandle< TileCondToolNoiseSample > m_tileToolNoiseSample
Applies OF algorithm.
Definition: TileRawChannelBuilderOpt2Filter.h:72
TileRawChannelBuilderOpt2Filter::ofc2int
void ofc2int(int nDigits, double *w_off, short *w_int, short &scale)
Definition: TileRawChannelBuilderOpt2Filter.cxx:717
TileRawChannel::assign
void assign(const HWIdentifier &id, float amplitude, float time, float quality, float ped=0.0)
Definition: TileRawChannel.h:63
TileHWID::ros
int ros(const HWIdentifier &id) const
extract ros field from HW identifier
Definition: TileHWID.h:167
TileRawChannelBuilderOpt2Filter::m_of2
bool m_of2
bool variable for OF method: true=> OF2; false=> OF1
Definition: TileRawChannelBuilderOpt2Filter.h:91
TileFragHash::OptFilterOffline
@ OptFilterOffline
Definition: TileFragHash.h:34
TileRawChannelContainer.h
xAOD::saturated
setScaleOne setStatusOne saturated
Definition: gFexGlobalRoI_v1.cxx:51
TileRawChannelBuilder::m_rChType
TileFragHash::TYPE m_rChType
Definition: TileRawChannelBuilder.h:136
TileHWID::adc
int adc(const HWIdentifier &id) const
extract adc field from HW identifier
Definition: TileHWID.h:193
TileOfcWeightsStruct::dg
float dg[99]
Definition: ITileCondToolOfc.h:14
TileInfo::ItrigSample
int ItrigSample() const
The sample at which the pulse should ideally peak.
Definition: TileInfo.h:77
TileRawChannelBuilder::m_nChL
int m_nChL
Definition: TileRawChannelBuilder.h:198
TileRawChannelBuilderOpt2Filter::m_emulateDsp
bool m_emulateDsp
Definition: TileRawChannelBuilderOpt2Filter.h:97
TileHWID.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
ascale
void ascale(TH1F *h, double s_)
Definition: comparitor.cxx:410
lumiFormat.i
int i
Definition: lumiFormat.py:92
TileRawChannelBuilder::m_bsflags
unsigned int m_bsflags
Definition: TileRawChannelBuilder.h:138
TileDigitsContainer.h
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
TileRawChannelBuilder::m_rawChannelContainerKey
SG::WriteHandleKey< TileRawChannelContainer > m_rawChannelContainerKey
Definition: TileRawChannelBuilder.h:129
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TileRawChannelBuilder::m_nChH
int m_nChH
Definition: TileRawChannelBuilder.h:199
TileOfcWeightsStruct
Definition: ITileCondToolOfc.h:13
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
TileRawChannel
Definition: TileRawChannel.h:35
test_pyathena.parent
parent
Definition: test_pyathena.py:15
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TileRawChannelBuilder::m_RChSumL
double m_RChSumL
Definition: TileRawChannelBuilder.h:200
TileRawChannelBuilderOpt2Filter::filter
double filter(int ros, int drawer, int channel, int &gain, double &pedestal, double &amplitude, double &time, const EventContext &ctx)
Definition: TileRawChannelBuilderOpt2Filter.cxx:340
TileRawChannelBuilder
Definition: TileRawChannelBuilder.h:59
maskDeadModules.ros
ros
Definition: maskDeadModules.py:35
DataPool::nextElementPtr
pointer nextElementPtr()
obtain the next available element in pool by pointer pool is resized if its limit has been reached On...
TileDigits::samples
const std::vector< float > & samples() const
Definition: TileDigits.h:58
TileRawChannelBuilderOpt2Filter::iterate
int iterate(int ros, int drawer, int channel, int gain, double &pedestal, double &amplitude, double &time, double &chi2, const EventContext &ctx)
Computes A,time,ped using OF.
Definition: TileRawChannelBuilderOpt2Filter.cxx:505
DataPool.h
min
#define min(a, b)
Definition: cfImp.cxx:40
TileRawChannelBuilder::m_correctTime
bool m_correctTime
Definition: TileRawChannelBuilder.h:145
TileRawChannelBuilder::m_tileInfo
const TileInfo * m_tileInfo
Definition: TileRawChannelBuilder.h:216
TileRawChannelBuilderOpt2Filter::getPedestal
float getPedestal(int ros, int drawer, int channel, int gain, const EventContext &ctx)
Apply the number of iterations needed for reconstruction by calling the Filter method.
Definition: TileRawChannelBuilderOpt2Filter.cxx:300
TileRawChannelBuilderOpt2Filter::rawChannel
virtual TileRawChannel * rawChannel(const TileDigits *digits, const EventContext &ctx)
Builder virtual method to be implemented by subclasses.
Definition: TileRawChannelBuilderOpt2Filter.cxx:192
TileRawChannelBuilder::correctAmp
static double correctAmp(double phase, bool of2=true)
Amplitude correction factor according to the time when using weights for tau=0 without iterations.
Definition: TileRawChannelBuilder.cxx:646
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
TileDigits
Definition: TileDigits.h:30
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
errorcheck.h
Helpers for checking error return status codes and reporting errors.
TileRawChannel::insertTime
void insertTime(float time)
Definition: TileRawChannel.cxx:76
TileRawChannelBuilder::m_i_ADCmaxPlus1
int m_i_ADCmaxPlus1
Definition: TileRawChannelBuilder.h:219
TileRawChannelBuilder::m_timeMaxThresh
float m_timeMaxThresh
correct amplitude is time is below time max threshold
Definition: TileRawChannelBuilder.h:154
TileRawChannelBuilderOpt2Filter::m_nCenter
int m_nCenter
internal counters
Definition: TileRawChannelBuilderOpt2Filter.h:100
TileRawChannelBuilderOpt2Filter::m_maxTime
double m_maxTime
max allowed time = 25*(m_nSamples-1)/2
Definition: TileRawChannelBuilderOpt2Filter.h:105
PlotSFuncertainty.calib
calib
Definition: PlotSFuncertainty.py:110
TileRawChannelBuilderOpt2Filter::m_nConst
int m_nConst
internal counters
Definition: TileRawChannelBuilderOpt2Filter.h:101
TileRawChannelBuilder::correctTime
static double correctTime(double phase, bool of2=true)
Time correction factor.
Definition: TileRawChannelBuilder.cxx:705
TileRawChannelBuilderOpt2Filter::m_maxIterations
int m_maxIterations
maximum number of iteration to perform
Definition: TileRawChannelBuilderOpt2Filter.h:87
TileRawChannelBuilderOpt2Filter::m_confTB
bool m_confTB
use testbeam configuration
Definition: TileRawChannelBuilderOpt2Filter.h:89
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
a
TList * a
Definition: liststreamerinfos.cxx:10
TileHWID::drawer
int drawer(const HWIdentifier &id) const
extract drawer field from HW identifier
Definition: TileHWID.h:171
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
TileRawChannelBuilderOpt2Filter::m_digits
std::vector< float > m_digits
Definition: TileRawChannelBuilderOpt2Filter.h:109
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
TileRawChannelBuilder::m_RChSumH
double m_RChSumH
Definition: TileRawChannelBuilder.h:201
AthCommonMsg< AlgTool >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
DataPool
a typed memory pool that saves time spent allocation small object. This is typically used by containe...
Definition: DataPool.h:47
TileRawChannelBuilderOpt2Filter::m_correctAmplitude
bool m_correctAmplitude
If true, resulting amplitude is corrected when using weights for tau=0 without iteration.
Definition: TileRawChannelBuilderOpt2Filter.h:93
TileRawChannelBuilderOpt2Filter::m_t0SamplePosition
int m_t0SamplePosition
position of peak sample = (m_nSamples-1)/2
Definition: TileRawChannelBuilderOpt2Filter.h:104
TileRawChannelBuilderOpt2Filter::m_noiseThresholdHG
int m_noiseThresholdHG
Definition: TileRawChannelBuilderOpt2Filter.h:111
TileRawChannelBuilder::m_firstSample
int m_firstSample
Definition: TileRawChannelBuilder.h:140
TileRawChannelBuilderOpt2Filter::~TileRawChannelBuilderOpt2Filter
~TileRawChannelBuilderOpt2Filter()
Destructor.
Definition: TileRawChannelBuilderOpt2Filter.cxx:86
TileOfcWeightsStruct::w_a
float w_a[99]
Definition: ITileCondToolOfc.h:14
TileCalibUtils::getDrawerIdx
static unsigned int getDrawerIdx(unsigned int ros, unsigned int drawer)
Returns a drawer hash.
Definition: TileCalibUtils.cxx:60
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
TileHWID::to_string
std::string to_string(const HWIdentifier &id, int level=0) const
extract all fields from HW identifier HWIdentifier get_all_fields ( const HWIdentifier & id,...
Definition: TileHWID.cxx:49
TileDigits.h
TileRawChannelBuilder::m_timeMinThresh
float m_timeMinThresh
correct amplitude is time is above time min threshold
Definition: TileRawChannelBuilder.h:153
python.compressB64.c
def c
Definition: compressB64.py:93
TileRawChannelBuilderOpt2Filter::m_nSamples
int m_nSamples
number of samples in the data
Definition: TileRawChannelBuilderOpt2Filter.h:103
TileRawChannelBuilderOpt2Filter::m_bestPhase
bool m_bestPhase
Definition: TileRawChannelBuilderOpt2Filter.h:96
readCCLHist.float
float
Definition: readCCLHist.py:83
TileRawChannelBuilderOpt2Filter::findMaxDigitPosition
int findMaxDigitPosition()
Finds maximum digit position in the pulse.
Definition: TileRawChannelBuilderOpt2Filter.cxx:270
TileRawChannelUnit::ADCcounts
@ ADCcounts
Definition: TileRawChannelUnit.h:17
TileOfcWeightsStruct::w_b
float w_b[99]
Definition: ITileCondToolOfc.h:14
TileRawChannelBuilderOpt2Filter::m_nNegative
int m_nNegative
internal counters
Definition: TileRawChannelBuilderOpt2Filter.h:99
TileRawChannelBuilderOpt2Filter::finalize
virtual StatusCode finalize()
Finalize method.
Definition: TileRawChannelBuilderOpt2Filter.cxx:170