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