ATLAS Offline Software
TileRawChannelBuilderFitFilter.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 // Tile includes
6 #include "TileEvent/TileRawChannel.h"
7 #include "TileRecUtils/TileRawChannelBuilderFitFilter.h"
8 #include "TileEvent/TileRawChannelContainer.h"
9 #include "TileEvent/TileDigitsContainer.h"
10 #include "TileEvent/TileDigits.h"
11 #include "CaloIdentifier/TileID.h"
12 #include "TileIdentifier/TileHWID.h"
13 #include "TileIdentifier/TileTrigType.h"
14 #include "TileConditions/TileInfo.h"
15 #include "TileConditions/TileCablingService.h"
16 
17 // Atlas includes
18 #include "AthAllocators/DataPool.h"
19 #include "AthenaKernel/errorcheck.h"
20 
21 // Gaudi includes
22 #include "Gaudi/Property.h"
23 
24 
25 // lang include
26 #include <algorithm>
27 #include <cmath>
28 
29 
30 static const InterfaceID IID_ITileRawChannelBuilderFitFilter
31  ("TileRawChannelBuilderFitFilter", 1, 0);
32 
33 const InterfaceID& TileRawChannelBuilderFitFilter::interfaceID( ) {
34  return IID_ITileRawChannelBuilderFitFilter;
35  //return TileRawChannelBuilder::interfaceID();
36 }
37 
41 TileRawChannelBuilderFitFilter::TileRawChannelBuilderFitFilter(const std::string& type
42  , const std::string& name, const IInterface *parent)
43  : TileRawChannelBuilder(type, name, parent)
44  , m_t0Fit(0.0)
45  , m_iPeak0(0)
46  , m_minTime(0.0)
47  , m_maxTime(0.0)
48  , m_minTau(0)
49  , m_maxTau(0.0)
50  , m_pulseShapes(0)
51 {
52  //declare interfaces
53  declareInterface< TileRawChannelBuilder >( this );
54  declareInterface< TileRawChannelBuilderFitFilter >(this);
55 
56  m_rawChannelContainerKey = "TileRawChannelFit";
57 
58  //declare properties
59  declareProperty("FrameLength",m_frameLength = 9);
60  declareProperty("MaxIterate",m_maxIterate = 9);
61  declareProperty("NoiseLowGain",m_noiseLow = 0.6);
62  declareProperty("NoiseHighGain",m_noiseHigh = 1.5);
63  declareProperty("RMSChannelNoise",m_channelNoiseRMS = 3);
64  declareProperty("ExtraSamplesLeft",m_extraSamplesLeft=0); // increase window on left side
65  declareProperty("ExtraSamplesRight",m_extraSamplesRight=0); // increase window on right side
66  declareProperty("SaturatedSample",m_saturatedSample = -1.0);
67  declareProperty("SaturatedSampleError",m_saturatedSampleError = 6.0);
68  declareProperty("ZeroSampleError",m_zeroSampleError = 100.0);
69  declareProperty("NoiseThresholdRMS",m_noiseThresholdRMS = 3.0);
70  declareProperty("MaxTimeFromPeak",m_maxTimeFromPeak = 250.0);
71 
72  declareProperty("DisableNegativeAmp",m_disableNegativeAmp = false);
73 
74  declareProperty("SpecialDemoShape",m_specialDemoShape = -1); // if >=0 - pulse shape for Demo is stored in non-default structures
75 }
76 
80 TileRawChannelBuilderFitFilter::~TileRawChannelBuilderFitFilter(){
81 }
82 
86 StatusCode TileRawChannelBuilderFitFilter::initialize() {
87 
88  ATH_MSG_INFO( "TileRawChannelBuilderFitFilter::initialize()" );
89 
90  m_rChType = TileFragHash::FitFilter;
91 
92  // init in superclass
93  ATH_CHECK( TileRawChannelBuilder::initialize() );
94  ATH_CHECK( m_tileToolNoiseSample.retrieve() );
95 
96  if (m_bestPhase) {
97  //=== get TileToolTiming
98  // TileToolTiming can be disabled in the TileRawChannelBuilder
99  if (!m_tileToolTiming.isEnabled()) {
100  m_tileToolTiming.enable();
101  }
102  ATH_CHECK( m_tileToolTiming.retrieve() );
103  }
104 
105  // Get number of samples from TileInfo - ignore jobOptions settings
106  m_frameLength = m_tileInfo->NdigitSamples();
107  ATH_MSG_DEBUG( "NSamples = " << m_frameLength );
108 
109  // Get pulse shapes from TileInfo
110  m_pulseShapes = m_tileInfo->getPulseShapes();
111 
112  // Determine peak sample position
113  // peak sample position defines t=0
114  // m_iPeak0 = (int) (m_frameLength) / 2 + (m_frameLength % 2) - 1;
115  m_iPeak0 = (int) m_tileInfo->ItrigSample();
116 
117  // Min and max time are now calculated based on m_framelength - i.e. on the
118  // number of 25 ns samples read out. Don't forget t=0 corresponds to
119  // m_ipeak0-th sample
120  m_minTime = DTIME * (0 - m_iPeak0 - m_extraSamplesLeft);
121  m_maxTime = DTIME * (m_frameLength - 1 - m_iPeak0 + m_extraSamplesRight);
122  // maximal jump during one iteration
123  m_maxTau = (m_maxTime - m_minTime);
124  m_minTau = -m_maxTau;
125 
126  ATH_MSG_DEBUG( " ipeak0=" << m_iPeak0
127  << " min_time=" << m_minTime
128  << " max_time=" << m_maxTime
129  << " min_tau=" << m_minTau
130  << " max_tau=" << m_maxTau );
131 
132  if ( !m_demoFragIDs.empty() ) {
133  switch (m_specialDemoShape) {
134  case 1:
135  ATH_MSG_DEBUG( "Demonstrator channels use pulse shape from physics structures"); break;
136  case 2:
137  ATH_MSG_DEBUG( "Demonstrator channels use pulse shape from laser structures"); break;
138  case 3:
139  ATH_MSG_DEBUG( "Demonstrator channels use pulse shape from laser(for 100pF) and physics(for 5.2pF) structures"); break;
140  case 8:
141  ATH_MSG_DEBUG( "Demonstrator channels use pulse shape from cis structures"); break;
142  default:
143  ATH_MSG_DEBUG( "Demonstrator channels use the same pulse shape as legacy"); break;
144  }
145  }
146 
147  // Speedup for physics processing (max_iter=1):
148  // read initial pulse shapes into arrays
149  m_fnParameters[0] = 0.0;
150  m_fnParameters[1] = 0.0;
151  m_fnParameters[2] = 1.0;
152  m_t0Fit = 0.0;
153 
154  m_gPhysLo.reserve(m_frameLength);
155  m_dgPhysLo.reserve(m_frameLength);
156  m_gPhysHi.reserve(m_frameLength);
157  m_dgPhysHi.reserve(m_frameLength);
158 
159  double rsamp;
160  for (int isamp = 0; isamp < m_frameLength; ++isamp) {
161  rsamp = (isamp) * 1.0;
162  m_gPhysLo.push_back(scaledPulse(rsamp, &(m_pulseShapes->m_tlphys), &(m_pulseShapes->m_ylphys)));
163  m_dgPhysLo.push_back(pulse(rsamp, &(m_pulseShapes->m_tdlphys), &(m_pulseShapes->m_ydlphys)));
164  m_gPhysHi.push_back(scaledPulse(rsamp, &(m_pulseShapes->m_thphys), &(m_pulseShapes->m_yhphys)));
165  m_dgPhysHi.push_back(pulse(rsamp, &(m_pulseShapes->m_tdhphys), &(m_pulseShapes->m_ydhphys)));
166  }
167 
168  if (m_channelNoiseRMS < 0 || m_channelNoiseRMS > 3) {
169  m_channelNoiseRMS = 3;
170  }
171 
172  if (msgLvl(MSG::DEBUG)) {
173  switch (m_channelNoiseRMS) {
174  case 1:
175  msg(MSG::DEBUG) << " noise for all channels from noiseOrig array (OBSOLETE!!!) " << endmsg;
176  break;
177  case 2:
178  msg(MSG::DEBUG) << " noise for all channels from noiseNk array (OBSOLETE!!!) " << endmsg;
179  break;
180  case 3:
181  msg(MSG::DEBUG) << " noise for all channels from Conditions DB ";
182  if (m_cabling->getTestBeam()) {
183  const EventContext &ctx = Gaudi::Hive::currentContext();
184  msg(MSG::DEBUG) << " rmsLow(LBA01/0) = " << m_tileToolNoiseSample->getHfn(20, 0, TileID::LOWGAIN, TileRawChannelUnit::ADCcounts, ctx)
185  << " rmsHi(LBA01/0) = " << m_tileToolNoiseSample->getHfn(20, 0, TileID::HIGHGAIN, TileRawChannelUnit::ADCcounts, ctx)
186  << endmsg;
187  } else {
188  msg(MSG::DEBUG) << endmsg;
189  }
190  break;
191  default:
192  msg(MSG::DEBUG) << " common noise for all channels (OBSOLETE): " << endmsg;
193  msg(MSG::DEBUG) << " rmsLow = " << m_noiseLow
194  << " rmsHi = " << m_noiseHigh << endmsg;
195  break;
196  }
197  }
198 
199  // TileInfo
200  if (m_saturatedSample < 0) m_saturatedSample = m_f_ADCmax;
201 
202  return StatusCode::SUCCESS;
203 }
204 
205 StatusCode TileRawChannelBuilderFitFilter::finalize() {
206 
207  ATH_MSG_DEBUG( "Finalizing" );
208  return StatusCode::SUCCESS;
209 }
210 
211 TileRawChannel* TileRawChannelBuilderFitFilter::rawChannel(const TileDigits* digits, const EventContext& ctx) {
212 
213 
214  ++m_chCounter;
215 
216  const HWIdentifier adcId = digits->adc_HWID();
217 
218  ATH_MSG_VERBOSE( "Running FitFilter for TileRawChannel with HWID "
219  << m_tileHWID->to_string(adcId) );
220 
221  // tmp variables for filtering
222  double amplitude = 0.0;
223  double chi2 = 0.0;
224  double time = 0.0;
225  double pedestal = 0.0;
226 
227  // use fit filter
228  pulseFit(digits, amplitude, time, pedestal, chi2, ctx);
229 
230  unsigned int drawerIdx(0), channel(0), adc(0);
231  m_tileIdTransforms->getIndices(adcId, drawerIdx, channel, adc);
232 
233  // fit filter calib
234  // note that when called from TileROD_Decoder, m_calibrateEnergy is set
235  // from TileROD_Decoder...
236  if (m_calibrateEnergy) {
237  amplitude = m_tileToolEmscale->doCalibCis(drawerIdx, channel, adc, amplitude);
238  }
239 
240  ATH_MSG_VERBOSE( "Creating RawChannel"
241  << " a=" << amplitude
242  << " t=" << time
243  << " q=" << chi2
244  << " ped=" << pedestal );
245 
246 
247  // return new TileRawChannel
248  // TileRawChannel *rawCh = new TileRawChannel(adcId,amplitude,time,chi2,pedestal);
249 
250  DataPool<TileRawChannel> tileRchPool(m_dataPoollSize);
251  TileRawChannel *rawCh = tileRchPool.nextElementPtr();
252  rawCh->assign (adcId,
253  amplitude,
254  time,
255  chi2,
256  pedestal);
257 
258  if (m_correctTime && ((chi2 > 0) || m_bestPhase)) {
259  time -= m_tileToolTiming->getSignalPhase(drawerIdx, channel, adc);
260  rawCh->insertTime(time);
261  ATH_MSG_VERBOSE( "Correcting time, new time=" << rawCh->time() );
262  }
263 
264  int gain = m_tileHWID->adc(adcId);
265  if (TileID::HIGHGAIN == gain) {
266  ++m_nChH;
267  m_RChSumH += amplitude;
268  } else {
269  ++m_nChL;
270  m_RChSumL += amplitude;
271  }
272 
273  /*
274  ATH_MSG_ALWAYS("TileRawChannel Id: " << m_tileHWID->to_string(adcId)
275  << " => amp/ped/time: " << rawCh->amplitude() << "/" << rawCh->pedestal() << "/" << rawCh->time());
276  */
277  return rawCh;
278 }
279 
285 void TileRawChannelBuilderFitFilter::pulseFit(const TileDigits *digit
286  , double &amplitude, double &time, double &pedestal, double &chi2, const EventContext &ctx) {
287 
288  amplitude = 0.0;
289  time = 0.0;
290  pedestal = 0.0;
291  chi2 = MAX_CHI2;
292 
293  const HWIdentifier adcId = digit->adc_HWID();
294  int ros = m_tileHWID->ros(adcId);
295  int channel = m_tileHWID->channel(adcId);
296  int igain = m_tileHWID->adc(adcId);
297 
298  int drawer = m_tileHWID->drawer(adcId);
299  unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer);
300 
301  bool demo = (m_specialDemoShape > 0) && std::binary_search(m_demoFragIDs.begin(), m_demoFragIDs.end(), (ros << 8) | drawer);
302 
303  // Estimate channel noise
304  double rms = 0.0;
305  int noise_channel = (ros < 3) ? channel : channel + 48;
306 
307  if (igain == 0) {
308  switch (m_channelNoiseRMS) {
309  case 3:
310  rms = m_tileToolNoiseSample->getHfn(drawerIdx, channel, igain, TileRawChannelUnit::ADCcounts, ctx);
311  if (rms > 0.0) break;
312  /* FALLTHROUGH */
313  case 2:
314  rms = m_pulseShapes->m_noiseNkLo[noise_channel];
315  if (rms > 0.0) break;
316  /* FALLTHROUGH */
317  case 1:
318  rms = m_pulseShapes->m_noiseOrigLo[noise_channel];
319  if (rms > 0.0) break;
320  /* FALLTHROUGH */
321  default:
322  rms = m_noiseLow;
323  }
324  } else if (igain == 1) {
325  switch (m_channelNoiseRMS) {
326  case 3:
327  rms = m_tileToolNoiseSample->getHfn(drawerIdx, channel, igain, TileRawChannelUnit::ADCcounts, ctx);
328  if (rms > 0.0) break;
329  /* FALLTHROUGH */
330  case 2:
331  rms = m_pulseShapes->m_noiseNkHi[noise_channel];
332  if (rms > 0.0) break;
333  /* FALLTHROUGH */
334  case 1:
335  rms = m_pulseShapes->m_noiseOrigHi[noise_channel];
336  if (rms > 0.0) break;
337  /* FALLTHROUGH */
338  default:
339  rms = m_noiseHigh;
340  }
341  } else {
342  // neither low nor hi-gain, (e.g. adder data)
343  chi2 = -1.0;
344  return;
345  }
346 
347  ATH_MSG_VERBOSE( "PulseFit:"
348  << " ROS=" << ros
349  << " Channel=" << channel
350  << " gain=" << igain
351  << " RMS=" << rms
352  << " RMSChNoise=" << m_channelNoiseRMS );
353 
354  // First sample to fit, number of samples to fit
355  int ifit1 = 0;
356  int nfit = std::min(m_frameLength, digit->NtimeSamples());
357 
358  int ipeakMax = m_iPeak0;
359  int ipeakMin = m_iPeak0;
360 
361  ATH_MSG_VERBOSE( "ipeak0=" << m_iPeak0
362  << " ifit1=" << ifit1
363  << " ifit2=" << ifit1 + nfit - 1
364  << " nfit=" << nfit
365  << " idolas=" << m_idolas
366  << " idocis=" << m_idocis
367  << " CISchan=" << m_cischan
368  << " capdaq=" << m_capdaq
369  << " demoCh=" << ((demo)?"true":"false") );
370 
371  std::vector<float> samples = digit->samples();
372  samples.erase(samples.begin(),samples.begin()+m_firstSample);
373  samples.resize(m_frameLength);
374  double maxsamp = 0.0;
375  double minsamp = m_saturatedSample;
376 
377  double xvec[MAX_SAMPLES], yvec0[MAX_SAMPLES];
378  double yvec[MAX_SAMPLES], eyvec[MAX_SAMPLES];
379 
380  bool no_signal = true;
381  pedestal = samples[0];
382  const double delta = 1.e-6;
383 
384  for (int isamp = 0; isamp < nfit; ++isamp) {
385  int j = isamp + ifit1;
386  xvec[isamp] = j;
387  yvec0[isamp] = samples[j];
388  if (no_signal) no_signal = (fabs(samples[j] - pedestal) < delta);
389 
390  // Check for saturated or zero samples and de-weight accordingly
391  if ((yvec0[isamp] < m_saturatedSample) && (yvec0[isamp] > 0)) {
392  eyvec[isamp] = rms;
393  } else {
394  if (yvec0[isamp] >= m_saturatedSample) {
395  eyvec[isamp] = m_saturatedSampleError * rms;
396  ATH_MSG_VERBOSE( "Saturated ADC value yvec0[" << isamp << "]=" << yvec0[isamp]
397  << " (MAX=" << m_saturatedSample << " ) RMS=" << eyvec[isamp] );
398 
399  } else { // must be yvec0[isamp]==0
400  eyvec[isamp] = m_zeroSampleError * rms;
401  ATH_MSG_VERBOSE( "Zero ADC value yvec0[" << isamp << "]=" << yvec0[isamp]
402  << " RMS=" << eyvec[isamp] );
403  }
404  }
405 
406  if (yvec0[isamp] > maxsamp) {
407  // initial time guess based on
408  // sample with maximum value
409  maxsamp = yvec0[isamp];
410  ipeakMax = j;
411  }
412  if (yvec0[isamp] < minsamp) {
413  minsamp = yvec0[isamp];
414  ipeakMin = j;
415  }
416  ATH_MSG_VERBOSE( "isamp=" << isamp
417  << ", xvec=" << xvec[isamp]
418  << ", yvec0=" << yvec0[isamp]
419  << ", eyvec=" << eyvec[isamp] );
420  }
421 
422  if (no_signal) {
423  ATH_MSG_VERBOSE( "No signal detected" );
424  return;
425  }
426 
427  // Make an initial guess about pedestal
428  double pedg = (yvec0[0] + yvec0[nfit - 1]) / 2.0;
429 
430  // Position of max sample compared to nominal peak position
431  int delta_peak = 0;
432  // Time offset in pulse functions
433  m_t0Fit = 0.0;
434  // Flag for fixed time in fit
435  bool fixedTime = (m_maxIterate < 0);
436 
437  if (!fixedTime) {
438  if (m_disableNegativeAmp) { // try to reproduce Opt Filter behaviour
439  // take first sample as pedestal - like it is done in Opt Filter method
440  pedg = yvec0[0];
441  if (maxsamp - pedg > m_noiseThresholdRMS * rms) {
442  delta_peak = ipeakMax - m_iPeak0; // Adjust initial phase guess,
443  m_t0Fit = (delta_peak) * DTIME; // positive amplitude
444  } else {
445  fixedTime = true; // no signal above noise
446  m_t0Fit = 0.0; // fit with fixed time
447  }
448  } else {
449  if (maxsamp - pedg > m_noiseThresholdRMS * rms) {
450  delta_peak = ipeakMax - m_iPeak0; // Adjust initial phase guess,
451  m_t0Fit = (delta_peak) * DTIME; // positive amplitude
452  } else if (pedg - minsamp > m_noiseThresholdRMS * rms) {
453  delta_peak = ipeakMin - m_iPeak0; // Adjust initial phase guess,
454  m_t0Fit = (delta_peak) * DTIME; // negative amplitude
455  } else {
456  fixedTime = true; // no signal above noise
457  m_t0Fit = 0.0; // fit with fixed time
458  }
459  }
460  }
461 
462  double expectedTime = 0.;
463  if (fixedTime && m_bestPhase) {
464  expectedTime = m_tileToolTiming->getSignalPhase(drawerIdx, channel, igain);
465  delta_peak = std::round(expectedTime / DTIME); // Adjust initial phase guess
466  m_t0Fit = expectedTime;
467  }
468 
469  ATH_MSG_VERBOSE ( " initial value of"
470  << " t0fit=" << m_t0Fit
471  << " ipeakMax=" << ipeakMax
472  << " ipeakMin=" << ipeakMin
473  << " fixedTime=" << ((fixedTime) ? "true" : "false") );
474 
475  const std::vector<double>* tpulse = &m_dummy;
476  const std::vector<double>* ypulse = &m_dummy;
477  const std::vector<double>* tdpulse = &m_dummy;
478  const std::vector<double>* dpulse = &m_dummy;
479  const std::vector<double>* tleak = &m_dummy;
480  const std::vector<double>* yleak = &m_dummy;
481  const std::vector<double>* tdleak = &m_dummy;
482  const std::vector<double>* dleak = &m_dummy;
483 
484 
485  bool docis = m_idocis;
486  bool dolas = m_idolas;
487  if (demo) { // special treatment for Demo drawers - select different pulse shape
488  dolas = ((m_specialDemoShape == 2) || (m_specialDemoShape == 3 && m_capdaq > 10));
489  docis = (m_specialDemoShape == 8);
490  }
491 
492  if (docis && ((m_cischan == -1) || (channel == m_cischan) || demo)) { // CIS pulse
493  if (igain == 0) { // low gain
494  if (m_capdaq > 10) { // 100 pF capacitor
495  tpulse = &(m_pulseShapes->m_tlcis);
496  ypulse = &(m_pulseShapes->m_ylcis);
497  tdpulse = &(m_pulseShapes->m_tdlcis);
498  dpulse = &(m_pulseShapes->m_ydlcis);
499  tleak = &(m_pulseShapes->m_tleaklo);
500  yleak = &(m_pulseShapes->m_leaklo);
501  tdleak = &(m_pulseShapes->m_tdleaklo);
502  dleak = &(m_pulseShapes->m_dleaklo);
503  } else { // 5.2 pF capacitor
504  tpulse = &(m_pulseShapes->m_tslcis);
505  ypulse = &(m_pulseShapes->m_yslcis);
506  tdpulse = &(m_pulseShapes->m_tdslcis);
507  dpulse = &(m_pulseShapes->m_ydslcis);
508  tleak = &(m_pulseShapes->m_tsleaklo);
509  yleak = &(m_pulseShapes->m_sleaklo);
510  tdleak = &(m_pulseShapes->m_tdsleaklo);
511  dleak = &(m_pulseShapes->m_dsleaklo);
512  }
513  } else { // igain==1 => high-gain
514  if (m_capdaq > 10) { // 100 pF capacitor
515  tpulse = &(m_pulseShapes->m_thcis);
516  ypulse = &(m_pulseShapes->m_yhcis);
517  tdpulse = &(m_pulseShapes->m_tdhcis);
518  dpulse = &(m_pulseShapes->m_ydhcis);
519  tleak = &(m_pulseShapes->m_tleakhi);
520  yleak = &(m_pulseShapes->m_leakhi);
521  tdleak = &(m_pulseShapes->m_tdleakhi);
522  dleak = &(m_pulseShapes->m_dleakhi);
523  } else { // 5.2 pF capacitor
524  tpulse = &(m_pulseShapes->m_tshcis);
525  ypulse = &(m_pulseShapes->m_yshcis);
526  tdpulse = &(m_pulseShapes->m_tdshcis);
527  dpulse = &(m_pulseShapes->m_ydshcis);
528  tleak = &(m_pulseShapes->m_tsleakhi);
529  yleak = &(m_pulseShapes->m_sleakhi);
530  tdleak = &(m_pulseShapes->m_tdsleakhi);
531  dleak = &(m_pulseShapes->m_dsleakhi);
532  }
533  }
534  } else {
535  if (dolas) { // laser pulse
536  if (igain == 0) { // low gain
537  tpulse = &(m_pulseShapes->m_tllas);
538  ypulse = &(m_pulseShapes->m_yllas);
539  tdpulse = &(m_pulseShapes->m_tdllas);
540  dpulse = &(m_pulseShapes->m_ydllas);
541  } else { // igain==1 => high-gain
542  tpulse = &(m_pulseShapes->m_thlas);
543  ypulse = &(m_pulseShapes->m_yhlas);
544  tdpulse = &(m_pulseShapes->m_tdhlas);
545  dpulse = &(m_pulseShapes->m_ydhlas);
546  }
547  } else { // physics pulse
548  if (igain == 0) { // low gain
549  tpulse = &(m_pulseShapes->m_tlphys);
550  ypulse = &(m_pulseShapes->m_ylphys);
551  tdpulse = &(m_pulseShapes->m_tdlphys);
552  dpulse = &(m_pulseShapes->m_ydlphys);
553  } else { // igain==1 => high-gain
554  tpulse = &(m_pulseShapes->m_thphys);
555  ypulse = &(m_pulseShapes->m_yhphys);
556  tdpulse = &(m_pulseShapes->m_tdhphys);
557  dpulse = &(m_pulseShapes->m_ydhphys);
558  }
559  }
560  }
561 
562  if (demo) {
563  docis = false; // never fit demo channels as CIS with signal and leakage pulses
564  dolas = (m_idocis || m_idolas); // use laser option, i.e. just single pulse without leakage pulse
565  }
566 
567  // Variables used for iterative fitting
568  double gval, gpval, sy, syg, sygp, sg, sgp, sgg, sgpgp, sggp, serr, err2;
569  double dgg0, dgg, dggp, dgpgp, dyg, dygp, dg, dc, xd;
570  double sllp, sylp, slplp, dleakage, leakage;
571  double fixtau = 0.0, fixped = 0.0, fixampl = 0.0, fixchi2 = MAX_CHI2;
572  double leaktau = 0.0, leakped = 0.0, leakampl = 0.0, leakchi2 = MAX_CHI2;
573  double cistau = 0.0, cisped = 0.0, cisampl = 0.0, cisatau = 0.0, cischi2 = MAX_CHI2;
574  double tau, ped, ampl, atau = 0.0, tempChi2 = MAX_CHI2, oldchi2 = MAX_CHI2 / 2;
575 
576  // number of iterations
577  int niter = 0;
578  do {
579  ++niter;
580  ATH_MSG_VERBOSE ( "niter=" << niter << " maxIterate=" << m_maxIterate );
581 
582  if (tempChi2 < oldchi2) oldchi2 = tempChi2; // best chi2 up to now
583 
584  // parameters for pulse shape functions
585  // 0. phase
586  m_fnParameters[0] = 0.0;
587  // 1. pedestal
588  m_fnParameters[1] = 0.0;
589  // 2. amplitude
590  m_fnParameters[2] = 1.0;
591 
592  // CIS events linear fit
593  if (docis && ((m_cischan == -1) || (channel == m_cischan))) {
594  ATH_MSG_VERBOSE ( "Fit time with leakage" );
595  // CIS Part (A): fit for time using leakage pulse
596  sllp = 0.0;
597  sylp = 0.0;
598  slplp = 0.0;
599  for (int isamp = 0; isamp < nfit; ++isamp) {
600  ATH_MSG_VERBOSE ( "Lo gain leakage xvec[" << isamp << "]=" << xvec[isamp] );
601 
602  leakage = pulse(xvec[isamp], tleak, yleak);
603  dleakage = derivative(xvec[isamp], tdleak, dleak);
604 
605  // Samples with pedestal subtracted
606  yvec[isamp] = yvec0[isamp] - pedg;
607 
608  ATH_MSG_VERBOSE( " yvec[" << isamp << "]=" << yvec[isamp]
609  << " yvec0[" << isamp << "]=" << yvec0[isamp]
610  << " pedg=" << pedg);
611 
612  sllp += leakage * dleakage;
613  sylp += yvec[isamp] * dleakage;
614  slplp += dleakage * dleakage;
615  }
616  // Also allow for fixed-time fit to CIS events
617  if (fabs(slplp) > EPS_DG && !fixedTime) {
618  leaktau = (sllp - sylp) / slplp;
619  // Also have to check the range for leaktau
620  if (leaktau > m_maxTau)
621  leaktau = m_maxTau;
622  else if (leaktau < m_minTau) leaktau = m_minTau;
623  } else {
624  leaktau = 0.0;
625  }
626 
627  ATH_MSG_VERBOSE( " sllp=" << sllp
628  << " sylp=" << sylp
629  << " slplp=" << slplp
630  << " leaktau=" << leaktau);
631 
632  // CIS Part (B): using phase determined in part (A),
633  // subtract leakage pedestal and fit for amplitude, pedestal
634  m_fnParameters[0] = leaktau;
635  sy = 0.0;
636  sg = 0.0;
637  syg = 0.0;
638  sgg = 0.0;
639  serr = 0.0;
640  for (int isamp = 0; isamp < nfit; ++isamp) {
641  leakage = pulse(xvec[isamp], tleak, yleak);
642  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
643  yvec[isamp] = yvec0[isamp] - leakage;
644 
645  ATH_MSG_VERBOSE ( " yvec[" << isamp << "]=" << yvec[isamp]
646  << " yvec0[" << isamp << "]=" << yvec0[isamp]
647  << " leakage=" << leakage );
648 
649  err2 = eyvec[isamp] * eyvec[isamp];
650  sy += yvec[isamp] / err2;
651  sg += gval / err2;
652  syg += yvec[isamp] * gval / err2;
653  sgg += gval * gval / err2;
654  serr += 1.0 / err2;
655  }
656 
657  dgg0 = sg * sg - serr * sgg;
658  if (fabs(dgg0) > EPS_DG) {
659  leakampl = (sy * sg - serr * syg) / dgg0;
660  leakped = (syg * sg - sy * sgg) / dgg0;
661  } else {
662  leakampl = 0.0;
663  leakped = sy / serr;
664  }
665 
666  // Determine Chi2 for corresponding function for CIS leakage + pulse
667  ATH_MSG_VERBOSE ( " Determine Chi2 for CIS leakage + pulse" );
668 
669  leakchi2 = 0.0;
670  m_fnParameters[0] = leaktau;
671  m_fnParameters[1] = leakped;
672  m_fnParameters[2] = leakampl;
673 
674  for (int isamp = 0; isamp < nfit; ++isamp) {
675  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
676  leakage = pulse(xvec[isamp], tleak, yleak);
677  xd = yvec0[isamp] - (gval + leakage);
678  leakchi2 = leakchi2 + (xd * xd) / (eyvec[isamp] * eyvec[isamp]);
679 
680  ATH_MSG_VERBOSE ( " isamp=" << isamp
681  << " yvec0[" << isamp << "]=" << yvec0[isamp]
682  << " gval=" << gval
683  << ", leakage=" << leakage
684  << ", xd=" << xd );
685  }
686 
687  leakchi2 = leakchi2/(nfit-3.0);
688 
689  ATH_MSG_VERBOSE ( " leaktau=" << leaktau
690  << " leakped=" << leakped
691  << " leakampl=" << leakampl
692  << " leakchi2=" << leakchi2 );
693 
694  // CIS Part C: Least-squares fit with 3 parameters for pulse+leakage
695  if (!fixedTime) {
696  m_fnParameters[0] = 0.0;
697  m_fnParameters[1] = 0.0;
698  m_fnParameters[2] = 0.0;
699 
700  for (int isamp = 0; isamp < nfit; ++isamp) {
701  leakage = pulse(xvec[isamp], tleak, yleak);
702 
703  // Subtract leakage from samples
704  yvec[isamp] = yvec0[isamp] - leakage;
705 
706  ATH_MSG_VERBOSE ( " isamp=" << isamp
707  << " yvec0[" << isamp << "]=" << yvec0[isamp]
708  << " leakage=" << leakage
709  << " yvec[" << isamp << "]=" << yvec[isamp] );
710  }
711 
712  m_fnParameters[0] = 0.0;
713  m_fnParameters[1] = 0.0;
714  m_fnParameters[2] = 1.0;
715  sy = 0.0;
716  sg = 0.0;
717  sgp = 0.0;
718  syg = 0.0;
719  sygp = 0.0;
720  sgg = 0.0;
721  sggp = 0.0;
722  sgpgp = 0.0;
723  serr = 0.0;
724 
725  for (int isamp = 0; isamp < nfit; ++isamp) {
726  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
727  gpval = derivative(xvec[isamp], tdpulse, dpulse);
728  err2 = eyvec[isamp] * eyvec[isamp];
729  sy += yvec[isamp] / err2;
730  sg += gval / err2;
731  sgp += gpval / err2;
732  syg += yvec[isamp] * gval / err2;
733  sygp += yvec[isamp] * gpval / err2;
734  sgg += gval * gval / err2;
735  sggp += gval * gpval / err2;
736  sgpgp += gpval * gpval / err2;
737  serr += 1.0 / err2;
738  }
739 
740  dgg = sgg - sg * sg / serr;
741  dggp = sggp - sg * sgp / serr;
742  dgpgp = sgpgp - sgp * sgp / serr;
743  dyg = syg - sy * sg / serr;
744  dygp = sygp - sy * sgp / serr;
745  dg = dgg * dgpgp - dggp * dggp;
746 
747  if (fabs(dg) > EPS_DG) {
748  cisampl = (dyg * dgpgp - dygp * dggp) / dg;
749  cisatau = (dyg * dggp - dygp * dgg) / dg;
750  cisped = (sy
751  - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) / dg)
752  / serr;
753 
754  if (fabs(cisampl) > EPS_DG) {
755  cistau = cisatau / cisampl;
756  if (cistau > m_maxTau)
757  cistau = m_maxTau;
758  else if (cistau < m_minTau) cistau = m_minTau;
759  } else {
760  cistau = 0.0;
761  }
762  } else {
763  cisampl = 0.0;
764  cisatau = 0.0;
765  cistau = 0.0;
766  cisped = sy / serr;
767  }
768 
769  if (msgLvl(MSG::VERBOSE)) {
770  msg(MSG::VERBOSE) << " sy=" << sy
771  << " sg=" << sg
772  << " sgp=" << sgp
773  << " syg=" << syg
774  << " sygp=" << sygp
775  << " sgg=" << sgg
776  << " sggp=" << sggp
777  << " sgpgp=" << sgpgp << endmsg;
778 
779  msg(MSG::VERBOSE) << " dgg=" << dgg
780  << " dggp=" << dggp
781  << " sgpgp=" << sgpgp
782  << " dyg=" << dyg
783  << " dygp=" << dygp
784  << " dg=" << dg << endmsg;
785 
786  msg(MSG::VERBOSE) << " cistau=" << cistau
787  << " cisped=" << cisped
788  << " cisampl=" << cisampl << endmsg;
789  }
790 
791  // Determine Chi2 for pulse shape + leakage fit CIS Part C
792  cischi2 = 0.0;
793  m_fnParameters[0] = cistau;
794  m_fnParameters[1] = cisped;
795  m_fnParameters[2] = cisampl;
796 
797  for (int isamp = 0; isamp < nfit; ++isamp) {
798  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
799  leakage = pulse(xvec[isamp], tleak, yleak);
800  // Subtract leakage from samples
801  yvec[isamp] = yvec0[isamp] - leakage;
802  xd = yvec[isamp] - gval;
803  cischi2 = cischi2 + (xd * xd) / (eyvec[isamp] * eyvec[isamp]);
804 
805  ATH_MSG_VERBOSE ( " yvec0[" << isamp << "]=" << yvec0[isamp]
806  << " yvec[" << isamp << "]=" << yvec[isamp]
807  << " leakage=" << leakage
808  << " gval=" << gval
809  << " xd=" << xd );
810  }
811 
812  cischi2 = cischi2 / (nfit - 3.0);
813 
814  ATH_MSG_VERBOSE ( " cischi2=" << cischi2 );
815  }
816 
817  // Determine which set of parameters to use from CIS fit methods based on minimum chi2
818  if ((cischi2 < leakchi2) && !fixedTime) {
819  tau = cistau;
820  ped = cisped;
821  ampl = cisampl;
822  tempChi2 = cischi2;
823  } else {
824  tau = leaktau;
825  ped = leakped;
826  ampl = leakampl;
827  tempChi2 = leakchi2;
828  }
829  // End of fit for CIS events
830  } else {
831  // Physics and laser events
832 
833  if (niter == 1) {
834  /* For first iteration, also calculate 2-Parameter Fit for pedestal and amplitude
835  */
836  double t0fit_old = m_t0Fit;
837  m_t0Fit = expectedTime;
838 
839  sy = 0.0;
840  sg = 0.0;
841  sgg = 0.0;
842  syg = 0.0;
843  serr = 0.0;
844 
845  for (int isamp = 0; isamp < nfit; ++isamp) {
846  if (!dolas) {
847  // Use initial values for speeding up the physics events
848  int jsamp = (int) xvec[isamp] - delta_peak;
849  if (jsamp < 0)
850  jsamp = 0;
851  else if (jsamp >= nfit) jsamp = nfit - 1;
852 
853  if (igain == 0) {
854  gval = m_gPhysLo[jsamp];
855  } else {
856  gval = m_gPhysHi[jsamp];
857  }
858  } else {
859  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
860  }
861  err2 = eyvec[isamp] * eyvec[isamp];
862  sy += yvec0[isamp] / err2;
863  sg += gval / err2;
864  syg += yvec0[isamp] * gval / err2;
865  sgg += gval * gval / err2;
866  serr += 1.0 / err2;
867  }
868 
869  fixtau = 0.0;
870  dgg0 = sg * sg - serr * sgg;
871  if (fabs(dgg0) > EPS_DG) {
872  fixampl = (sy * sg - serr * syg) / dgg0;
873  fixped = (syg * sg - sy * sgg) / dgg0;
874  ATH_MSG_VERBOSE ( " 2-par fit: fixampl = " << fixampl
875  << " fixped = " << fixped );
876  } else {
877  fixampl = 0.0;
878  fixped = sy / serr;
879  ATH_MSG_VERBOSE ( " 2-par fit: small dgg0 = " << dgg0
880  << ", fixampl = " << fixampl
881  << " fixped = " << fixped );
882  }
883 
884  m_fnParameters[0] = fixtau; /* 2-Par fit Calculate chi2 for physics events */
885  m_fnParameters[1] = fixped;
886  m_fnParameters[2] = fixampl;
887  fixchi2 = 0.0;
888  for (int isamp = 0; isamp < nfit; ++isamp) {
889  dc = yvec0[isamp] - scaledPulse(xvec[isamp], tpulse, ypulse);
890  fixchi2 = fixchi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
891  ATH_MSG_VERBOSE ( " isamp= " << isamp
892  << " yvec0[" << isamp << "]= " << yvec0[isamp]
893  << " eyvec[" << isamp << "]= " << eyvec[isamp]
894  << " fixchi2= " << fixchi2 );
895  }
896  fixchi2 = fixchi2 / (nfit - 2.0);
897  ATH_MSG_VERBOSE ( " fixchi2/(nfit-2.0)=" << fixchi2
898  << " nfit=" << nfit );
899 
900  m_t0Fit = t0fit_old;
901 
902  // restore initial parameters for pulse shape functions - to be used in 3-par fit
903  m_fnParameters[0] = 0.0;
904  m_fnParameters[1] = 0.0;
905  m_fnParameters[2] = 1.0;
906  } /* end of 2-par fit in first iteration */
907 
908  if (fixedTime) {
909  m_t0Fit = expectedTime;
910  tau = fixtau;
911  ped = fixped;
912  ampl = fixampl;
913  tempChi2 = oldchi2 = -fabs(fixchi2);
914  } else {
915 
916  sy = 0.0;
917  sg = 0.0;
918  sgp = 0.0;
919  syg = 0.0;
920  sygp = 0.0;
921  sgg = 0.0;
922  sggp = 0.0;
923  sgpgp = 0.0;
924  serr = 0.0;
925 
926  for (int isamp = 0; isamp < nfit; ++isamp) {
927  if ((niter == 1) && (!dolas)) {
928  // Use initial function values stored in array for niter=1 physics
929  // XXX: double->int
930  int jsamp = (int) xvec[isamp] - delta_peak;
931  if (jsamp < 0)
932  jsamp = 0;
933  else if (jsamp >= nfit) jsamp = nfit - 1;
934 
935  if (igain == 0) { // igain ==0 => low-gain
936  gval = m_gPhysLo[jsamp];
937  gpval = m_dgPhysLo[jsamp];
938  } else { // must be igain==1 => high-gain
939  gval = m_gPhysHi[jsamp];
940  gpval = m_dgPhysHi[jsamp];
941  }
942  } else {
943  // Use the respective function values
944  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
945  gpval = derivative(xvec[isamp], tdpulse, dpulse);
946  }
947 
948  err2 = eyvec[isamp] * eyvec[isamp];
949  sy += yvec0[isamp] / err2;
950  sg += gval / err2;
951  sgp += gpval / err2;
952  syg += yvec0[isamp] * gval / err2;
953  sygp += yvec0[isamp] * gpval / err2;
954  sgg += gval * gval / err2;
955  sggp += gval * gpval / err2;
956  sgpgp += gpval * gpval / err2;
957  serr += 1.0 / err2;
958 
959  ATH_MSG_VERBOSE ( " isamp=" << isamp
960  << " gval=" << gval
961  << " sg=" << sg
962  << " gpval=" << gpval
963  << " sgp=" << sgp );
964  }
965 
966  dgg = sgg - sg * sg / serr;
967  dggp = sggp - sg * sgp / serr;
968  dgpgp = sgpgp - sgp * sgp / serr;
969  dyg = syg - sy * sg / serr;
970  dygp = sygp - sy * sgp / serr;
971  dg = dgg * dgpgp - dggp * dggp;
972 
973  // Fit for time, pedestal, and amplitude
974  if (fabs(dg) > EPS_DG) {
975  // Amplitude : ampl
976  ampl = (dyg * dgpgp - dygp * dggp) / dg;
977  // and Amplitude * time: atau
978  atau = (dyg * dggp - dygp * dgg) / dg;
979  // Pedestal
980  ped = (sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) / dg)
981  / serr;
982 
983  if (fabs(ampl) > EPS_DG) {
984  // Time
985  tau = atau / ampl;
986  if (tau > m_maxTau)
987  tau = m_maxTau;
988  else if (tau < m_minTau) tau = m_minTau;
989  } else {
990  tau = 0.0;
991  }
992  } else {
993  ampl = 0.0;
994  atau = 0.0;
995  tau = 0.0;
996  ped = sy / serr;
997  }
998 
999  if (msgLvl(MSG::VERBOSE)) {
1000  msg(MSG::VERBOSE) << " ped=" << ped << endmsg;
1001  msg(MSG::VERBOSE) << " sy=" << sy
1002  << " sg=" << sg
1003  << " sgp=" << sgp << endmsg;
1004 
1005  msg(MSG::VERBOSE) << " syg=" << syg
1006  << " sygp=" << sygp
1007  << " sgg=" << sgg << endmsg;
1008 
1009  msg(MSG::VERBOSE) << " sggp=" << sggp
1010  << " sgpgp=" << sgpgp << endmsg;
1011 
1012  msg(MSG::VERBOSE) << " ampl = (dyg*dgpgp - dygp*dggp)= " << ampl << endmsg;
1013 
1014  msg(MSG::VERBOSE) << " dyg=" << dyg
1015  << " dgpgp=" << dgpgp
1016  << " dyg*dgpgp=" << (dyg * dgpgp) << endmsg;
1017 
1018  msg(MSG::VERBOSE) << " dygp=" << dygp
1019  << " dggp=" << dggp
1020  << " dygp*dggp=" << (dygp * dggp) << endmsg;
1021 
1022  msg(MSG::VERBOSE) << " dyg=" << dyg
1023  << " dggp=" << dggp
1024  << " dyg*dggp=" << (dyg * dggp) << endmsg;
1025 
1026  msg(MSG::VERBOSE) << " dygp=" << dygp
1027  << " dgg=" << dgg
1028  << " dygp*dgg=" << (dygp * dgg) << endmsg;
1029 
1030  msg(MSG::VERBOSE) << " dg=" << dg
1031  << " atau=" << atau
1032  << " tau=" << tau << endmsg;
1033  }
1034 
1035  m_fnParameters[0] = tau;
1036  m_fnParameters[1] = ped;
1037  m_fnParameters[2] = ampl;
1038 
1039  tempChi2 = 0;
1040  // Calculate chi2 for physics and laser events
1041  for (int isamp = 0; isamp < nfit; ++isamp) {
1042  dc = yvec0[isamp] - scaledPulse(xvec[isamp], tpulse, ypulse);
1043  tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1044  ATH_MSG_VERBOSE ( " isamp=" << isamp
1045  << " yvec0[" << isamp << "]=" << yvec0[isamp]
1046  << " eyvec[" << isamp << "]=" << eyvec[isamp]
1047  << " dc=" << dc
1048  << " chi2=" << tempChi2 );
1049  }
1050 
1051  tempChi2 = tempChi2 / (nfit - 3.0);
1052  ATH_MSG_VERBOSE ( " chi2/(nfit-3.0)=" << tempChi2
1053  << " nfit=" << nfit );
1054  } // end if fixedTime
1055  } // end of physics and laser specific part
1056 
1057  if (msgLvl(MSG::VERBOSE))
1058  msg(MSG::VERBOSE) << " t0fit: " << m_t0Fit << ((tau < 0.0) ? " - " : " + ") << fabs(tau);
1059 
1060  // avoid infinite loop at the boundary
1061  if (fabs(m_t0Fit - m_maxTime) < 0.001 && tau >= 0.0) { // trying to go outside max boundary second time
1062  m_t0Fit = fixtau + (m_minTime - fixtau) * niter / m_maxIterate; // jump to negative time
1063  if (msgLvl(MSG::VERBOSE))
1064  msg(MSG::VERBOSE) << " jumping to " << m_t0Fit << endmsg;
1065  tempChi2 = MAX_CHI2;
1066  } else if (fabs(m_t0Fit - m_minTime) < 0.001 && tau <= 0.0) { // trying to go outside min boundary second time
1067  m_t0Fit = fixtau + (m_maxTime - fixtau) * niter / m_maxIterate; // jump to positive time
1068  if (msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << " jumping to " << m_t0Fit << endmsg;
1069  tempChi2 = MAX_CHI2;
1070  } else {
1071 
1072  // Iteration with parameter for time
1073  m_t0Fit += tau;
1074  if (msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << " = " << m_t0Fit << endmsg;
1075 
1076  // Check if total time does not exceed the limits:
1077  if (m_t0Fit > m_maxTime) {
1078  m_t0Fit = m_maxTime;
1079  tempChi2 = MAX_CHI2;
1080  } else if (m_t0Fit < m_minTime) {
1081  m_t0Fit = m_minTime;
1082  tempChi2 = MAX_CHI2;
1083  }
1084  }
1085 
1086  // save values of the best iteration
1087  if (tempChi2 < chi2) {
1088  time = m_t0Fit;
1089  pedestal = ped;
1090  amplitude = ampl;
1091  chi2 = tempChi2;
1092  }
1093 
1094  if (!fixedTime && tempChi2 < MAX_CHI2 && fabs(tau) < EPS_DG) { // too small tau
1095  if (m_t0Fit > fixtau)
1096  m_t0Fit = fixtau + (m_minTime - fixtau) * niter / m_maxIterate; // jump to negative time
1097  else
1098  m_t0Fit = fixtau + (m_maxTime - fixtau) * niter / m_maxIterate; // jump to positive time
1099 
1100  ATH_MSG_VERBOSE( " too small tau - jump to " << m_t0Fit);
1101  }
1102 
1103  ATH_MSG_VERBOSE ( " iter=" << niter
1104  << " t0fit=" << m_t0Fit
1105  << " phase=" << tau
1106  << " ped=" << ped
1107  << " ampl=" << ampl
1108  << " chi2=" << tempChi2 );
1109 
1110  } while (fabs(tempChi2 - oldchi2) > DELTA_CHI2 && (niter < m_maxIterate));
1111 
1112 // NGO never use the 2par fit result if non-pedestal event was detected!
1113 // if ((fabs(fixchi2) <= fabs(p_chi2))
1114 // && !(docis && ((m_cischan == -1) || (channel == m_cischan)))) {
1115 // /* results from 2-par fit */
1116 // p_time = fixtau;
1117 // p_pedestal = fixped;
1118 // p_amplitude = fixampl;
1119 // p_chi2 = -fabs(fixchi2);
1120 // }
1121 
1122 // TD: fit converged, now one extra iteration, leaving out the samples that
1123 // are more then m_MaxTimeFromPeak ns from the peak. We want to avoid to
1124 // extrapolate the pulse shape beyond the region where it is known.
1125 
1126 // Do it only in case when at least last sample is beyond m_MaxTimeFromPeak:
1127  if ((nfit - 1 - m_iPeak0) * DTIME > time + m_maxTimeFromPeak) {
1128 
1129  ATH_MSG_VERBOSE ( "Result before last iteration:"
1130  << " Time=" << time
1131  << " Ped=" << pedestal
1132  << " Amplitude=" << amplitude
1133  << " Chi2=" << chi2 );
1134 
1135  m_t0Fit = time;
1136  int nfit_real;
1137 
1138  // parameters for pulse shape functions
1139  // 0. phase
1140  m_fnParameters[0] = 0.0;
1141  // 1. pedestal
1142  m_fnParameters[1] = 0.0;
1143  // 2. amplitude
1144  m_fnParameters[2] = 1.0;
1145 
1146  // CIS events linear fit
1147  if (docis && ((m_cischan == -1) || (channel == m_cischan))) {
1148  if (!fixedTime) {
1149  ATH_MSG_VERBOSE ( "Fit time with leakage" );
1150  // CIS Part (A): fit for time using leakage pulse
1151  sllp = 0.0;
1152  sylp = 0.0;
1153  slplp = 0.0;
1154  for (int isamp = 0; (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1155 
1156  ATH_MSG_VERBOSE ( "Lo gain leakage xvec[" << isamp << "]=" << xvec[isamp] );
1157 
1158  leakage = pulse(xvec[isamp], tleak, yleak);
1159  dleakage = derivative(xvec[isamp], tdleak, dleak);
1160 
1161  // Samples with pedestal subtracted
1162  yvec[isamp] = yvec0[isamp] - pedg;
1163 
1164  ATH_MSG_VERBOSE ( " yvec[" << isamp << "]=" << yvec[isamp]
1165  << " yvec0[" << isamp << "]=" << yvec0[isamp]
1166  << " pedg=" << pedg );
1167 
1168  sllp += leakage * dleakage;
1169  sylp += yvec[isamp] * dleakage;
1170  slplp += dleakage * dleakage;
1171  }
1172  // Also allow for fixed-time fit to CIS events
1173  if (fabs(slplp) > EPS_DG && !fixedTime) {
1174  leaktau = (sllp - sylp) / slplp;
1175  // Also have to check the range for leaktau
1176  if (leaktau > m_maxTau)
1177  leaktau = m_maxTau;
1178  else if (leaktau < m_minTau) leaktau = m_minTau;
1179  } else {
1180  leaktau = 0.0;
1181  }
1182 
1183  ATH_MSG_VERBOSE ( " sllp=" << sllp
1184  << " sylp=" << sylp
1185  << " slplp=" << slplp
1186  << " leaktau=" << leaktau );
1187 
1188  // CIS Part (B): using phase determined in part (A),
1189  // subtract leakage pedestal and fit for amplitude, pedestal
1190  m_fnParameters[0] = leaktau;
1191  sy = 0.0;
1192  sg = 0.0;
1193  syg = 0.0;
1194  sgg = 0.0;
1195  serr = 0.0;
1196 
1197  for (int isamp = 0;
1198  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1199 
1200  leakage = pulse(xvec[isamp], tleak, yleak);
1201  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
1202  yvec[isamp] = yvec0[isamp] - leakage;
1203 
1204  ATH_MSG_VERBOSE ( " yvec[" << isamp << "]=" << yvec[isamp]
1205  << " yvec0[" << isamp << "]=" << yvec0[isamp]
1206  << " leakage=" << leakage );
1207 
1208  err2 = eyvec[isamp] * eyvec[isamp];
1209  sy += yvec[isamp] / err2;
1210  sg += gval / err2;
1211  syg += yvec[isamp] * gval / err2;
1212  sgg += gval * gval / err2;
1213  serr += 1.0 / err2;
1214  }
1215 
1216  dgg0 = sg * sg - serr * sgg;
1217  if (fabs(dgg0) > EPS_DG) {
1218  leakampl = (sy * sg - serr * syg) / dgg0;
1219  leakped = (syg * sg - sy * sgg) / dgg0;
1220  } else {
1221  leakampl = 0.0;
1222  leakped = sy / serr;
1223  }
1224 
1225  // Determine Chi2 for corresponding function for CIS leakage + pulse
1226  ATH_MSG_VERBOSE ( " Determine Chi2 for CIS leakage + pulse" );
1227 
1228  leakchi2 = 0.0;
1229  nfit_real = 0;
1230  m_fnParameters[0] = leaktau;
1231  m_fnParameters[1] = leakped;
1232  m_fnParameters[2] = leakampl;
1233 
1234  for (int isamp = 0;
1235  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1236 
1237  ++nfit_real;
1238  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
1239  leakage = pulse(xvec[isamp], tleak, yleak);
1240  xd = yvec0[isamp] - (gval + leakage);
1241  leakchi2 = leakchi2 + (xd * xd) / (eyvec[isamp] * eyvec[isamp]);
1242 
1243  ATH_MSG_VERBOSE ( " isamp=" << isamp
1244  << " yvec0[" << isamp << "]=" << yvec0[isamp]
1245  << " gval=" << gval
1246  << ", leakage=" << leakage
1247  << ", xd=" << xd );
1248  }
1249 
1250  leakchi2 = leakchi2 / (nfit_real - 3.0);
1251 
1252  ATH_MSG_VERBOSE ( " leaktau=" << leaktau
1253  << " leakped=" << leakped
1254  << " leakampl=" << leakampl
1255  << " leakchi2=" << leakchi2 );
1256 
1257  // CIS Part C: Least-squares fit with 3 parameters for pulse+leakage
1258  m_fnParameters[0] = 0.0;
1259  m_fnParameters[1] = 0.0;
1260  m_fnParameters[2] = 0.0;
1261 
1262  for (int isamp = 0;
1263  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1264 
1265  leakage = pulse(xvec[isamp], tleak, yleak);
1266 
1267  // Subtract leakage from samples
1268  yvec[isamp] = yvec0[isamp] - leakage;
1269 
1270  ATH_MSG_VERBOSE ( " isamp=" << isamp
1271  << " yvec0[" << isamp << "]=" << yvec0[isamp]
1272  << " leakage=" << leakage
1273  << " yvec[" << isamp << "]=" << yvec[isamp] );
1274  }
1275 
1276  m_fnParameters[0] = 0.0;
1277  m_fnParameters[1] = 0.0;
1278  m_fnParameters[2] = 1.0;
1279  sy = 0.0;
1280  sg = 0.0;
1281  sgp = 0.0;
1282  syg = 0.0;
1283  sygp = 0.0;
1284  sgg = 0.0;
1285  sggp = 0.0;
1286  sgpgp = 0.0;
1287  serr = 0.0;
1288 
1289  for (int isamp = 0;
1290  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1291 
1292  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
1293  gpval = derivative(xvec[isamp], tdpulse, dpulse);
1294  err2 = eyvec[isamp] * eyvec[isamp];
1295  sy += yvec[isamp] / err2;
1296  sg += gval / err2;
1297  sgp += gpval / err2;
1298  syg += yvec[isamp] * gval / err2;
1299  sygp += yvec[isamp] * gpval / err2;
1300  sgg += gval * gval / err2;
1301  sggp += gval * gpval / err2;
1302  sgpgp += gpval * gpval / err2;
1303  serr += 1.0 / err2;
1304  }
1305 
1306  dgg = sgg - sg * sg / serr;
1307  dggp = sggp - sg * sgp / serr;
1308  dgpgp = sgpgp - sgp * sgp / serr;
1309  dyg = syg - sy * sg / serr;
1310  dygp = sygp - sy * sgp / serr;
1311  dg = dgg * dgpgp - dggp * dggp;
1312 
1313  if (fabs(dg) > EPS_DG) {
1314  cisampl = (dyg * dgpgp - dygp * dggp) / dg;
1315  cisatau = (dyg * dggp - dygp * dgg) / dg;
1316  cisped = (sy
1317  - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) / dg)
1318  / serr;
1319 
1320  if (fabs(cisampl) > EPS_DG) {
1321  cistau = cisatau / cisampl;
1322  if (cistau > m_maxTau)
1323  cistau = m_maxTau;
1324  else if (cistau < m_minTau) cistau = m_minTau;
1325  } else {
1326  cistau = 0.0;
1327  }
1328  } else {
1329  cisampl = 0.0;
1330  cisatau = 0.0;
1331  cistau = 0.0;
1332  cisped = sy / serr;
1333  }
1334 
1335  if (msgLvl(MSG::VERBOSE)) {
1336  msg(MSG::VERBOSE) << " sy=" << sy
1337  << " sg=" << sg
1338  << " sgp=" << sgp
1339  << " syg=" << syg
1340  << " sygp=" << sygp
1341  << " sgg=" << sgg
1342  << " sggp=" << sggp
1343  << " sgpgp=" << sgpgp << endmsg;
1344 
1345  msg(MSG::VERBOSE) << " dgg=" << dgg
1346  << " dggp=" << dggp
1347  << " sgpgp=" << sgpgp
1348  << " dyg=" << dyg
1349  << " dygp=" << dygp
1350  << " dg=" << dg << endmsg;
1351 
1352  msg(MSG::VERBOSE) << " cistau=" << cistau
1353  << " cisped=" << cisped
1354  << " cisampl=" << cisampl << endmsg;
1355  }
1356 
1357  // Determine Chi2 for pulse shape + leakage fit CIS
1358  cischi2 = 0.0;
1359  nfit_real = 0;
1360  m_fnParameters[0] = cistau;
1361  m_fnParameters[1] = cisped;
1362  m_fnParameters[2] = cisampl;
1363 
1364  for (int isamp = 0;
1365  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1366 
1367  ++nfit_real;
1368  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
1369  leakage = pulse(xvec[isamp], tleak, yleak);
1370  // Subtract leakage from samples
1371  yvec[isamp] = yvec0[isamp] - leakage;
1372  xd = yvec[isamp] - gval;
1373  cischi2 = cischi2 + (xd * xd) / (eyvec[isamp] * eyvec[isamp]);
1374 
1375  ATH_MSG_VERBOSE ( " yvec0[" << isamp << "]=" << yvec0[isamp]
1376  << " yvec[" << isamp << "]=" << yvec[isamp]
1377  << " leakage=" << leakage
1378  << " gval=" << gval
1379  << " xd=" << xd );
1380  }
1381 
1382  cischi2=cischi2/(nfit_real-3.0);
1383 
1384  ATH_MSG_VERBOSE ( " cischi2=" << cischi2 );
1385 
1386  // Determine which set of parameters to use from CIS fit methods based on minimum chi2
1387  if ((cischi2 < leakchi2) && !fixedTime) {
1388  tau = cistau;
1389  ped = cisped;
1390  ampl = cisampl;
1391  tempChi2 = cischi2;
1392  } else {
1393  tau = leaktau;
1394  ped = leakped;
1395  ampl = leakampl;
1396  tempChi2 = leakchi2;
1397  }
1398  // End of fit for CIS events
1399  }
1400  } else { // Physics and laser events
1401  if (!fixedTime) {
1402 
1403  // restore initial parameters for pulse shape functions - to be used in 3-par fit
1404  m_fnParameters[0] = 0.0;
1405  m_fnParameters[1] = 0.0;
1406  m_fnParameters[2] = 1.0;
1407 
1408  sy = 0.0;
1409  sg = 0.0;
1410  sgp = 0.0;
1411  syg = 0.0;
1412  sygp = 0.0;
1413  sgg = 0.0;
1414  sggp = 0.0;
1415  sgpgp = 0.0;
1416  serr = 0.0;
1417 
1418  for (int isamp = 0;
1419  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1420  // Use the respective function values
1421  gval = scaledPulse(xvec[isamp], tpulse, ypulse);
1422  gpval = derivative(xvec[isamp], tdpulse, dpulse);
1423 
1424  err2 = eyvec[isamp] * eyvec[isamp];
1425  sy += yvec0[isamp] / err2;
1426  sg += gval / err2;
1427  sgp += gpval / err2;
1428  syg += yvec0[isamp] * gval / err2;
1429  sygp += yvec0[isamp] * gpval / err2;
1430  sgg += gval * gval / err2;
1431  sggp += gval * gpval / err2;
1432  sgpgp += gpval * gpval / err2;
1433  serr += 1.0 / err2;
1434 
1435  ATH_MSG_VERBOSE ( " isamp=" << isamp
1436  << " gval=" << gval
1437  << " sg=" << sg
1438  << " gpval=" << gpval
1439  << " sgp=" << sgp );
1440  }
1441 
1442  dgg = sgg - sg * sg / serr;
1443  dggp = sggp - sg * sgp / serr;
1444  dgpgp = sgpgp - sgp * sgp / serr;
1445  dyg = syg - sy * sg / serr;
1446  dygp = sygp - sy * sgp / serr;
1447  dg = dgg * dgpgp - dggp * dggp;
1448 
1449  // Fit for time, pedestal, and amplitude
1450  if (fabs(dg) > EPS_DG) {
1451  // Amplitude : ampl
1452  ampl = (dyg * dgpgp - dygp * dggp) / dg;
1453  // and Amplitude * time: atau
1454  atau = (dyg * dggp - dygp * dgg) / dg;
1455  // Pedestal
1456  ped = (sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) / dg)
1457  / serr;
1458 
1459  if (fabs(ampl) > EPS_DG) {
1460  // Time
1461  tau = atau / ampl;
1462  if (tau > m_maxTau)
1463  tau = m_maxTau;
1464  else if (tau < m_minTau) tau = m_minTau;
1465  } else {
1466  tau = 0.0;
1467  }
1468  } else {
1469  ampl = 0.0;
1470  atau = 0.0;
1471  tau = 0.0;
1472  ped = sy / serr;
1473  }
1474 
1475  if (msgLvl(MSG::VERBOSE)) {
1476  msg(MSG::VERBOSE) << " ped=" << ped << endmsg;
1477  msg(MSG::VERBOSE) << " sy=" << sy
1478  << " sg=" << sg
1479  << " sgp=" << sgp << endmsg;
1480 
1481  msg(MSG::VERBOSE) << " syg=" << syg
1482  << " sygp=" << sygp
1483  << " sgg=" << sgg << endmsg;
1484 
1485  msg(MSG::VERBOSE) << " sggp=" << sggp
1486  << " sgpgp=" << sgpgp << endmsg;
1487 
1488  msg(MSG::VERBOSE) << " ampl = (dyg*dgpgp - dygp*dggp)= " << ampl << endmsg;
1489  msg(MSG::VERBOSE) << " dyg=" << dyg
1490  << " dgpgp=" << dgpgp
1491  << " dyg*dgpgp=" << (dyg * dgpgp) << endmsg;
1492 
1493  msg(MSG::VERBOSE) << " dygp=" << dygp
1494  << " dggp=" << dggp
1495  << " dygp*dggp=" << (dygp * dggp) << endmsg;
1496 
1497  msg(MSG::VERBOSE) << " dyg=" << dyg
1498  << " dggp=" << dggp
1499  << " dyg*dggp=" << (dyg * dggp) << endmsg;
1500 
1501  msg(MSG::VERBOSE) << " dygp=" << dygp
1502  << " dgg=" << dgg
1503  << " dygp*dgg=" << (dygp * dgg) << endmsg;
1504 
1505  msg(MSG::VERBOSE) << " dg=" << dg
1506  << " atau=" << atau
1507  << " tau=" << tau << endmsg;
1508  }
1509 
1510  m_fnParameters[0] = tau;
1511  m_fnParameters[1] = ped;
1512  m_fnParameters[2] = ampl;
1513 
1514  tempChi2 = 0;
1515  nfit_real = 0;
1516  // Calculate chi2 for physics and laser events
1517  for (int isamp = 0;
1518  (isamp < nfit) && (DTIME * (isamp - m_iPeak0) - m_t0Fit < m_maxTimeFromPeak); ++isamp) {
1519 
1520  ++nfit_real;
1521  dc = yvec0[isamp] - scaledPulse(xvec[isamp], tpulse, ypulse);
1522  tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1523  ATH_MSG_VERBOSE ( " isamp=" << isamp
1524  << " yvec0[" << isamp << "]=" << yvec0[isamp]
1525  << " eyvec[" << isamp << "]=" << eyvec[isamp]
1526  << " dc=" << dc
1527  << " chi2=" << tempChi2 );
1528  }
1529 
1530  tempChi2 = tempChi2 / (nfit_real - 3.0);
1531  ATH_MSG_VERBOSE ( " chi2/(nfit_real-3.0)=" << tempChi2
1532  << " nfit_real=" << nfit_real );
1533  } // end if fixedTime
1534  } // end of physics and laser specific part
1535 
1536  if (msgLvl(MSG::VERBOSE))
1537  msg(MSG::VERBOSE) << " t0fit: " << m_t0Fit << ((tau < 0.0) ? " - " : " + ") << fabs(tau);
1538 
1539  // Iteration with parameter for time
1540  m_t0Fit += tau;
1541  if ( msgLvl(MSG::VERBOSE) )
1542  msg(MSG::VERBOSE) << " = " << m_t0Fit << endmsg;
1543 
1544  // Check if total time does not exceed the limits:
1545  if (m_t0Fit > m_maxTime) {
1546  m_t0Fit = m_maxTime;
1547  tempChi2 = MAX_CHI2;
1548  } else if (m_t0Fit < m_minTime) {
1549  m_t0Fit = m_minTime;
1550  tempChi2 = MAX_CHI2;
1551  }
1552 
1553  if (tempChi2 < MAX_CHI2) {
1554  time = m_t0Fit;
1555  pedestal = ped;
1556  amplitude = ampl;
1557  chi2 = tempChi2;
1558  } // otherwise using the previous iteration
1559 
1560  } // end if to use extra iteration
1561 
1562  ATH_MSG_VERBOSE ( "Result: Time=" << time
1563  << " Ped=" << pedestal
1564  << " Amplitude=" << amplitude
1565  << " Chi2=" << chi2 );
1566 }
1567 
1568 
1572 double TileRawChannelBuilderFitFilter::pulse(double x, const std::vector<double> * xvec
1573  , const std::vector<double> * yvec, bool zeroOutside) const {
1574 
1575  int size1 = xvec->size() - 1;
1576  if (size1 < 1) return 0.0;
1577 
1578  const double delta = 1.e-6;
1579 
1580  double xpmin = xvec->at(0);
1581  double xpmax = xvec->at(size1);
1582 
1583  double xp = (x - m_iPeak0) * DTIME - m_t0Fit - m_fnParameters[0];
1584 
1585  double val = 0.0;
1586  double tdiv = (xpmax - xpmin) / size1;
1587 
1588  if (xp < xpmin + delta) {
1589  if (zeroOutside && xp < xpmin - delta)
1590  val = 0.0;
1591  else
1592  val = yvec->at(0);
1593 #ifdef EXTRAPOLATE_TO_ZERO
1594  if (xp < xpmin - delta && val != 0.0) {
1595  double newval = val + ((yvec->at(1) - val) / tdiv) * (xp - xpmin);
1596  if (val * newval < 0.0) {
1597  val = 0.0;
1598  } else if (fabs(newval) < fabs(val)) {
1599  val = newval;
1600  }
1601  }
1602 #endif
1603  } else if (xp > xpmax - delta) {
1604  if (zeroOutside && xp > xpmax + delta)
1605  val = 0.0;
1606  else
1607  val = yvec->at(size1);
1608 #ifdef EXTRAPOLATE_TO_ZERO
1609  if (xp > xpmax + delta && val != 0.0) {
1610  double newval = val + ((yvec->at(size1 - 1) - val) / tdiv) * (xp - xpmax);
1611  if (val * newval < 0.0) {
1612  val = 0.0;
1613  } else if (fabs(newval) < fabs(val)) {
1614  val = newval;
1615  }
1616  }
1617 #endif
1618  } else {
1619  int j = (int) ((xp - xpmin) / tdiv);
1620  val = yvec->at(j) + ((yvec->at(j + 1) - yvec->at(j)) / tdiv) * (xp - xvec->at(j));
1621  }
1622 
1623  return val;
1624 }
TilePulseShapesStruct::m_yllas
std::vector< double > m_yllas
Definition: TilePulseShapes.h:63
TileRawChannelBuilderFitFilter::m_gPhysHi
std::vector< double > m_gPhysHi
Definition: TileRawChannelBuilderFitFilter.h:103
TilePulseShapesStruct::m_ydhphys
std::vector< double > m_ydhphys
Definition: TilePulseShapes.h:97
TilePulseShapesStruct::m_tllas
std::vector< double > m_tllas
Low Gain Pulse Laser.
Definition: TilePulseShapes.h:62
TileRawChannelBuilderFitFilter.h
TileRawChannelBuilderFitFilter::m_dgPhysLo
std::vector< double > m_dgPhysLo
Definition: TileRawChannelBuilderFitFilter.h:102
TilePulseShapesStruct::m_tlphys
std::vector< double > m_tlphys
Low Gain Pulse Physics.
Definition: TilePulseShapes.h:54
fitman.sy
sy
Definition: fitman.py:524
TilePulseShapesStruct::m_tdleakhi
std::vector< double > m_tdleakhi
Definition: TilePulseShapes.h:86
TileRawChannelBuilderFitFilter::m_minTime
double m_minTime
Definition: TileRawChannelBuilderFitFilter.h:95
TileRawChannelBuilder::m_dataPoollSize
int m_dataPoollSize
Definition: TileRawChannelBuilder.h:204
TileRawChannel.h
TileRawChannelBuilderFitFilter::m_noiseThresholdRMS
double m_noiseThresholdRMS
Definition: TileRawChannelBuilderFitFilter.h:115
checkCoolLatestUpdate.dg
dg
Definition: checkCoolLatestUpdate.py:9
TilePulseShapesStruct::m_ylphys
std::vector< double > m_ylphys
Definition: TilePulseShapes.h:55
plotting.yearwise_efficiency.channel
channel
Definition: yearwise_efficiency.py:28
TileRawChannelBuilder::m_tileToolEmscale
ToolHandle< TileCondToolEmscale > m_tileToolEmscale
Definition: TileRawChannelBuilder.h:166
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
TileRawChannelBuilder::m_chCounter
unsigned int m_chCounter
Definition: TileRawChannelBuilder.h:196
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TileInfo::NdigitSamples
int NdigitSamples() const
Returns the number of sammples (digits) per event.
Definition: TileInfo.h:75
TileRawChannelBuilderFitFilter::m_maxTimeFromPeak
double m_maxTimeFromPeak
Definition: TileRawChannelBuilderFitFilter.h:116
TilePulseShapesStruct::m_tdllas
std::vector< double > m_tdllas
Low Gain Pulse Laser.
Definition: TilePulseShapes.h:100
TileRawChannelBuilder::m_calibrateEnergy
bool m_calibrateEnergy
Definition: TileRawChannelBuilder.h:142
TileRawChannelBuilder::m_cabling
const TileCablingService * m_cabling
TileCabling instance.
Definition: TileRawChannelBuilder.h:182
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
TilePulseShapesStruct::m_ydlphys
std::vector< double > m_ydlphys
Definition: TilePulseShapes.h:93
TilePulseShapesStruct::m_dleakhi
std::vector< double > m_dleakhi
Definition: TilePulseShapes.h:87
TileFragHash::FitFilter
@ FitFilter
Definition: TileFragHash.h:35
TileRawChannelBuilderFitFilter::m_channelNoiseRMS
int m_channelNoiseRMS
Definition: TileRawChannelBuilderFitFilter.h:108
TileRawChannelBuilder::m_tileHWID
const TileHWID * m_tileHWID
Definition: TileRawChannelBuilder.h:161
TilePulseShapesStruct::m_ylcis
std::vector< double > m_ylcis
Definition: TilePulseShapes.h:35
CaloCondBlobAlgs_fillNoiseFromASCII.gain
gain
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:110
TilePulseShapesStruct::m_yhphys
std::vector< double > m_yhphys
Definition: TilePulseShapes.h:59
MuonGM::round
float round(const float toRound, const unsigned int decimals)
Definition: Mdt.cxx:27
TileRawChannelBuilderFitFilter::m_bestPhase
Gaudi::Property< bool > m_bestPhase
Definition: TileRawChannelBuilderFitFilter.h:71
TileRawChannelBuilderFitFilter::pulseFit
void pulseFit(const TileDigits *digit, double &amplitude, double &time, double &pedestal, double &chi2, const EventContext &ctx)
Calculate energy, time and chi2 for one channel using fitted pulse shape.
Definition: TileRawChannelBuilderFitFilter.cxx:285
TileRawChannelBuilder::initialize
virtual StatusCode initialize()
Initializer.
Definition: TileRawChannelBuilder.cxx:98
TilePulseShapesStruct::m_tsleaklo
std::vector< double > m_tsleaklo
Definition: TilePulseShapes.h:40
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
TileInfo.h
TilePulseShapesStruct::m_tdlphys
std::vector< double > m_tdlphys
Low Gain Pulse Physics.
Definition: TilePulseShapes.h:92
Tile_Base_ID::HIGHGAIN
@ HIGHGAIN
Definition: Tile_Base_ID.h:57
TilePulseShapesStruct::m_tleakhi
std::vector< double > m_tleakhi
Definition: TilePulseShapes.h:48
TilePulseShapesStruct::m_noiseNkHi
std::vector< double > m_noiseNkHi
Definition: TilePulseShapes.h:115
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
checkRpcDigits.digit
digit
Definition: checkRpcDigits.py:186
TilePulseShapesStruct::m_tdhcis
std::vector< double > m_tdhcis
Hi Gain Pulse CIS.
Definition: TilePulseShapes.h:82
TileRawData::adc_HWID
HWIdentifier adc_HWID(void) const
Definition: TileRawData.h:53
HWIdentifier
Definition: HWIdentifier.h:13
TileRawChannelBuilderFitFilter::m_noiseHigh
double m_noiseHigh
Definition: TileRawChannelBuilderFitFilter.h:119
TileRawChannelBuilderFitFilter::m_dummy
std::vector< double > m_dummy
Definition: TileRawChannelBuilderFitFilter.h:86
x
#define x
TileRawChannelBuilder::m_demoFragIDs
Gaudi::Property< std::vector< int > > m_demoFragIDs
Definition: TileRawChannelBuilder.h:184
TileRawChannelBuilderFitFilter::interfaceID
static const InterfaceID & interfaceID()
AlgTool InterfaceID.
Definition: TileRawChannelBuilderFitFilter.cxx:33
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
TilePulseShapesStruct::m_thphys
std::vector< double > m_thphys
Hi Gain Pulse Physics.
Definition: TilePulseShapes.h:58
TilePulseShapesStruct::m_yhlas
std::vector< double > m_yhlas
Definition: TilePulseShapes.h:67
TileRawChannelBuilder::m_cischan
int m_cischan
Definition: TileRawChannelBuilder.h:192
TileHWID::channel
int channel(const HWIdentifier &id) const
extract channel field from HW identifier
Definition: TileHWID.h:189
TileID.h
TileRawChannelBuilderFitFilter::m_extraSamplesLeft
int m_extraSamplesLeft
Definition: TileRawChannelBuilderFitFilter.h:110
TilePulseShapesStruct::m_tslcis
std::vector< double > m_tslcis
Definition: TilePulseShapes.h:36
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
TileRawChannelBuilderFitFilter::m_pulseShapes
const TilePulseShapesStruct * m_pulseShapes
Definition: TileRawChannelBuilderFitFilter.h:122
TilePulseShapesStruct::m_ydlcis
std::vector< double > m_ydlcis
Definition: TilePulseShapes.h:73
TileRawChannelContainer.h
TilePulseShapesStruct::m_dsleaklo
std::vector< double > m_dsleaklo
Definition: TilePulseShapes.h:79
TilePulseShapesStruct::m_ydslcis
std::vector< double > m_ydslcis
Definition: TilePulseShapes.h:75
TilePulseShapesStruct::m_yhcis
std::vector< double > m_yhcis
Definition: TilePulseShapes.h:45
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
TilePulseShapesStruct::m_tshcis
std::vector< double > m_tshcis
Definition: TilePulseShapes.h:46
TilePulseShapesStruct::m_dleaklo
std::vector< double > m_dleaklo
Definition: TilePulseShapes.h:77
TileRawChannelBuilderFitFilter::m_t0Fit
double m_t0Fit
Definition: TileRawChannelBuilderFitFilter.h:89
TileInfo::ItrigSample
int ItrigSample() const
The sample at which the pulse should ideally peak.
Definition: TileInfo.h:77
TileRawChannelBuilderFitFilter::rawChannel
virtual TileRawChannel * rawChannel(const TileDigits *digits, const EventContext &ctx) override
Builder virtual method to be implemented by subclasses.
Definition: TileRawChannelBuilderFitFilter.cxx:211
TileRawChannelBuilderFitFilter::m_fnParameters
double m_fnParameters[3]
Definition: TileRawChannelBuilderFitFilter.h:92
EPS_DG
#define EPS_DG
Definition: TilePulseShapes.h:26
TileRawChannelBuilder::m_nChL
int m_nChL
Definition: TileRawChannelBuilder.h:198
TilePulseShapesStruct::m_yshcis
std::vector< double > m_yshcis
Definition: TilePulseShapes.h:47
TileHWID.h
TileRawChannelBuilderFitFilter::m_iPeak0
int m_iPeak0
Definition: TileRawChannelBuilderFitFilter.h:93
TilePulseShapesStruct::m_ydllas
std::vector< double > m_ydllas
Definition: TilePulseShapes.h:101
TileCablingService.h
TilePulseShapesStruct::m_noiseNkLo
std::vector< double > m_noiseNkLo
(2) Noise with resistors added to PMT channels (so-called noise-killers)
Definition: TilePulseShapes.h:114
TilePulseShapesStruct::m_leakhi
std::vector< double > m_leakhi
Definition: TilePulseShapes.h:49
TileDigitsContainer.h
TileRawChannelBuilderFitFilter::m_tileToolNoiseSample
ToolHandle< TileCondToolNoiseSample > m_tileToolNoiseSample
Definition: TileRawChannelBuilderFitFilter.h:127
TileRawChannelBuilder::m_capdaq
double m_capdaq
Definition: TileRawChannelBuilder.h:193
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
TilePulseShapesStruct::m_thcis
std::vector< double > m_thcis
Hi Gain Pulse CIS.
Definition: TilePulseShapes.h:44
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
01SubmitToGrid.samples
samples
Definition: 01SubmitToGrid.py:58
TileRawChannelBuilder::m_tileIdTransforms
ToolHandle< TileCondIdTransforms > m_tileIdTransforms
Definition: TileRawChannelBuilder.h:172
TileRawChannelBuilderFitFilter::m_maxTau
double m_maxTau
Definition: TileRawChannelBuilderFitFilter.h:98
TileRawChannelBuilderFitFilter::m_saturatedSample
double m_saturatedSample
Definition: TileRawChannelBuilderFitFilter.h:112
TileRawChannelBuilder::m_nChH
int m_nChH
Definition: TileRawChannelBuilder.h:199
TilePulseShapesStruct::m_leaklo
std::vector< double > m_leaklo
Definition: TilePulseShapes.h:39
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
TileRawChannelBuilderFitFilter::m_noiseLow
double m_noiseLow
Definition: TileRawChannelBuilderFitFilter.h:118
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TileRawChannelBuilder::m_RChSumL
double m_RChSumL
Definition: TileRawChannelBuilder.h:200
TilePulseShapesStruct::m_ydhlas
std::vector< double > m_ydhlas
Definition: TilePulseShapes.h:105
Tile_Base_ID::LOWGAIN
@ LOWGAIN
Definition: Tile_Base_ID.h:57
WriteCellNoiseToCool.igain
igain
Definition: WriteCellNoiseToCool.py:338
TilePulseShapesStruct::m_noiseOrigHi
std::vector< double > m_noiseOrigHi
Definition: TilePulseShapes.h:111
TileRawChannelBuilderFitFilter::m_saturatedSampleError
double m_saturatedSampleError
Definition: TileRawChannelBuilderFitFilter.h:113
TileRawChannelBuilderFitFilter::m_gPhysLo
std::vector< double > m_gPhysLo
Definition: TileRawChannelBuilderFitFilter.h:101
TileRawChannelBuilder
Definition: TileRawChannelBuilder.h:59
maskDeadModules.ros
ros
Definition: maskDeadModules.py:35
plotBeamSpotCompare.xd
xd
Definition: plotBeamSpotCompare.py:220
DataPool::nextElementPtr
pointer nextElementPtr()
obtain the next available element in pool by pointer pool is resized if its limit has been reached On...
TilePulseShapesStruct::m_dsleakhi
std::vector< double > m_dsleakhi
Definition: TilePulseShapes.h:89
TileRawChannelBuilderFitFilter::pulse
double pulse(double x, const std::vector< double > *xvec, const std::vector< double > *yvec, bool zeroOutside=false) const
pulse interpolation
Definition: TileRawChannelBuilderFitFilter.cxx:1572
TilePulseShapesStruct::m_tdhlas
std::vector< double > m_tdhlas
Hi Gain Pulse Laser.
Definition: TilePulseShapes.h:104
DataPool.h
TilePulseShapesStruct::m_tdsleaklo
std::vector< double > m_tdsleaklo
Definition: TilePulseShapes.h:78
min
#define min(a, b)
Definition: cfImp.cxx:40
TileRawChannelBuilder::m_correctTime
bool m_correctTime
Definition: TileRawChannelBuilder.h:145
TilePulseShapesStruct::m_tdhphys
std::vector< double > m_tdhphys
Hi Gain Pulse Physics.
Definition: TilePulseShapes.h:96
TileRawChannelBuilder::m_tileInfo
const TileInfo * m_tileInfo
Definition: TileRawChannelBuilder.h:216
TileRawChannelBuilderFitFilter::m_specialDemoShape
int m_specialDemoShape
Definition: TileRawChannelBuilderFitFilter.h:125
TileRawChannelBuilderFitFilter::m_extraSamplesRight
int m_extraSamplesRight
Definition: TileRawChannelBuilderFitFilter.h:111
TileRawChannelBuilderFitFilter::initialize
virtual StatusCode initialize() override
Initializer.
Definition: TileRawChannelBuilderFitFilter.cxx:86
DTIME
#define DTIME
Definition: TileRawChannelBuilderFitFilter.h:116
TileRawChannelBuilder::m_f_ADCmax
float m_f_ADCmax
Definition: TileRawChannelBuilder.h:218
TileInfo::getPulseShapes
const TilePulseShapesStruct * getPulseShapes() const
Return pointer to TilePulseShapes.
Definition: TileInfo.h:253
TilePulseShapesStruct::m_ydshcis
std::vector< double > m_ydshcis
Definition: TilePulseShapes.h:85
TileCablingService::getTestBeam
bool getTestBeam() const
Definition: TileCablingService.h:274
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
TileTrigType.h
TileDigits
Definition: TileDigits.h:30
errorcheck.h
Helpers for checking error return status codes and reporting errors.
TileRawChannel::insertTime
void insertTime(float time)
Definition: TileRawChannel.cxx:76
TileRawChannelBuilderFitFilter::m_dgPhysHi
std::vector< double > m_dgPhysHi
Definition: TileRawChannelBuilderFitFilter.h:104
DELTA_CHI2
#define DELTA_CHI2
Definition: TilePulseShapes.h:28
TilePulseShapesStruct::m_tleaklo
std::vector< double > m_tleaklo
Definition: TilePulseShapes.h:38
TilePulseShapesStruct::m_tdshcis
std::vector< double > m_tdshcis
Definition: TilePulseShapes.h:84
Example_ReadSampleNoise.ped
ped
Definition: Example_ReadSampleNoise.py:45
MAX_SAMPLES
#define MAX_SAMPLES
Definition: TilePulseShapes.h:14
TilePulseShapesStruct::m_tlcis
std::vector< double > m_tlcis
Low Gain Pulse CIS.
Definition: TilePulseShapes.h:34
TileRawChannelBuilderFitFilter::m_maxIterate
int m_maxIterate
Definition: TileRawChannelBuilderFitFilter.h:109
TilePulseShapesStruct::m_sleaklo
std::vector< double > m_sleaklo
Definition: TilePulseShapes.h:41
TilePulseShapesStruct::m_tsleakhi
std::vector< double > m_tsleakhi
Definition: TilePulseShapes.h:50
TileRawChannelBuilderFitFilter::m_zeroSampleError
double m_zeroSampleError
Definition: TileRawChannelBuilderFitFilter.h:114
TilePulseShapesStruct::m_yslcis
std::vector< double > m_yslcis
Definition: TilePulseShapes.h:37
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
ReadFloatFromCool.adc
adc
Definition: ReadFloatFromCool.py:48
TilePulseShapesStruct::m_tdslcis
std::vector< double > m_tdslcis
Definition: TilePulseShapes.h:74
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
TilePulseShapesStruct::m_sleakhi
std::vector< double > m_sleakhi
Definition: TilePulseShapes.h:51
TileRawChannelBuilder::m_RChSumH
double m_RChSumH
Definition: TileRawChannelBuilder.h:201
AthCommonMsg< AlgTool >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
beamspotnt.rms
rms
Definition: bin/beamspotnt.py:1266
DataPool
a typed memory pool that saves time spent allocation small object. This is typically used by containe...
Definition: DataPool.h:47
TileRawChannelBuilder::m_firstSample
int m_firstSample
Definition: TileRawChannelBuilder.h:140
TileRawChannelBuilderFitFilter::m_minTau
double m_minTau
Definition: TileRawChannelBuilderFitFilter.h:97
TileRawChannelBuilderFitFilter::derivative
double derivative(double x, const std::vector< double > *xvec, const std::vector< double > *yvec) const
Definition: TileRawChannelBuilderFitFilter.h:82
TileCalibUtils::getDrawerIdx
static unsigned int getDrawerIdx(unsigned int ros, unsigned int drawer)
Returns a drawer hash.
Definition: TileCalibUtils.cxx:60
TileRawChannelBuilder::m_idocis
bool m_idocis
Definition: TileRawChannelBuilder.h:191
TilePulseShapesStruct::m_tdsleakhi
std::vector< double > m_tdsleakhi
Definition: TilePulseShapes.h:88
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
python.output.AtlRunQueryRoot.yvec
yvec
Definition: AtlRunQueryRoot.py:989
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
TileRawChannelBuilderFitFilter::~TileRawChannelBuilderFitFilter
~TileRawChannelBuilderFitFilter()
Destructor.
Definition: TileRawChannelBuilderFitFilter.cxx:80
TileRawChannelBuilderFitFilter::m_disableNegativeAmp
bool m_disableNegativeAmp
Definition: TileRawChannelBuilderFitFilter.h:124
MAX_CHI2
#define MAX_CHI2
Definition: TilePulseShapes.h:29
TileDigits.h
TileRawChannelBuilderFitFilter::scaledPulse
double scaledPulse(double x, const std::vector< double > *xvec, const std::vector< double > *yvec) const
Definition: TileRawChannelBuilderFitFilter.h:78
TileRawChannelBuilderFitFilter::finalize
virtual StatusCode finalize() override
Definition: TileRawChannelBuilderFitFilter.cxx:205
TilePulseShapesStruct::m_thlas
std::vector< double > m_thlas
Hi Gain Pulse Laser.
Definition: TilePulseShapes.h:66
TilePulseShapesStruct::m_tdlcis
std::vector< double > m_tdlcis
Low Gain Pulse CIS.
Definition: TilePulseShapes.h:72
TileRawChannelBuilder::m_idolas
bool m_idolas
Definition: TileRawChannelBuilder.h:189
TileRawChannelBuilderFitFilter::TileRawChannelBuilderFitFilter
TileRawChannelBuilderFitFilter(const std::string &type, const std::string &name, const IInterface *parent)
Constructor.
Definition: TileRawChannelBuilderFitFilter.cxx:41
TilePulseShapesStruct::m_tdleaklo
std::vector< double > m_tdleaklo
Definition: TilePulseShapes.h:76
TileRawChannelUnit::ADCcounts
@ ADCcounts
Definition: TileRawChannelUnit.h:17
TileRawChannelBuilderFitFilter::m_maxTime
double m_maxTime
Definition: TileRawChannelBuilderFitFilter.h:96
TileRawChannelBuilderFitFilter::m_frameLength
int m_frameLength
Definition: TileRawChannelBuilderFitFilter.h:107
TilePulseShapesStruct::m_ydhcis
std::vector< double > m_ydhcis
Definition: TilePulseShapes.h:83
TilePulseShapesStruct::m_noiseOrigLo
std::vector< double > m_noiseOrigLo
(1) Original noise
Definition: TilePulseShapes.h:110
Generated on Sun Jun 30 2024 21:29:06 for ATLAS Offline Software by   doxygen 1.8.18