105{
107
108
110 SG::ReadHandle<LArDigitContainer> digitContainer (
m_DataLocation, ctx);
111 ATH_MSG_DEBUG(
"1) LArDigitContainer container size = " << digitContainer->size() );
112
113 ATH_MSG_DEBUG(
"2) LArDigitContainer container size = " << digitContainer->size() );
114 if( digitContainer->empty() ) {
116 return StatusCode::SUCCESS;
117 }
118
119 auto larRawChannelContainer = std::make_unique<LArRawChannelContainer>();
120
121
122 const ILArPedestal* larPedestal=nullptr;
125 larPedestal=nullptr;
126 ATH_MSG_DEBUG(
"No pedestal found in database. Use default values." );
127 }
128 }
129
130 SG::ReadCondHandle<LArOnOffIdMapping> cablingHdl{
m_cablingKey, ctx};
131 const LArOnOffIdMapping*
cabling{*cablingHdl};
132 if(!cabling) {
134 return StatusCode::FAILURE;
135 }
136
137 const LArADC2MeV* adc2mev = nullptr;
139 SG::ReadCondHandle<LArADC2MeV> adc2mevH (
m_adc2mevKey, ctx);
140 adc2mev = *adc2mevH;
141 }
142
143
144
145
146
147
148
149 int nMinEM(0),nMinHEC(0),nMinFCAL(0);
150 std::vector<float> fSumEM,fSumHEC,fSumFCAL;
151 std::vector<float> *pSum = nullptr;
152
153 fSumEM.resize(0);
154 fSumHEC.resize(0);
155 fSumFCAL.resize(0);
156
157 for(int iloop=0;iloop<2;iloop++) {
158
159
161
162 for (const LArDigit* digit : *digitContainer) {
163
164
165 const std::vector<short>& samples =
digit->samples();
166 unsigned int nSamples = samples.size();
167 const HWIdentifier chid =
digit->channelID();
169
170
171 float thePedestal=-1;
172 if (larPedestal) {
173 float DBpedestal =larPedestal->
pedestal(chid,gain);
175 thePedestal=DBpedestal;
176 }
177 if (thePedestal<0) {
179 ATH_MSG_DEBUG(
"No pedestal found for this cell. Use default value " << thePedestal );
180 }
181
182 std::vector<float> mySamples;
183 mySamples.resize(samples.size());
184
185
186
187
188 for (
unsigned int i=0;
i<samples.size();
i++ )
189 {
190 mySamples[
i] = ((
float)samples[
i]) - thePedestal;
191 }
192
193
194 float GainFactor;
196 GainFactor = 9.8*9.8;
198 GainFactor = 9.8;
200 GainFactor = 1.0;
201 } else {
202 GainFactor = 1.0;
203 ATH_MSG_ERROR(
"Channel " << chid <<
"unknown gain: " << gain );
204 }
205
206
209 ATH_MSG_DEBUG( std::hex <<
" FebID / chid / channel = " << FebID <<
" / " << chid <<
" / " << channel );
210
211
212
214 bool isEM = false;
215 bool isHEC = false;
216 int nMin = 0;
217 unsigned int nAverage = 1;
218 float myScale = 1;
219 try {
220 const Identifier
id =
cabling->cnvToIdentifier(chid);
223 nMin = nMinFCAL;
224 pSum = &fSumFCAL;
227 }
228 else if (
m_emId->is_lar_em(
id)) {
229 isEM = true;
230 nMin = nMinEM;
231 pSum = &fSumEM;
234 }
235 else if (
m_hecId->is_lar_hec(
id)) {
236 isHEC = true;
237 nMin = nMinHEC;
238 pSum = &fSumHEC;
241 }
242 }
243
244 catch (LArID_Exception & execpt) {
246
247
248 continue;
249 }
250
251 if(!pSum) continue;
252
253
254
255
256
257 if ( nAverage > nSamples ) {
259 << nAverage << ") is larger than total number of samples ("
260 << nSamples << ")! adjusting nAverage ... " );
262 }
263
264
265 if ( iloop == 0 ) {
266
267 if ( pSum->empty())
268
269 pSum->resize(nSamples,0);
270
271 for(
unsigned i=0;
i<
nSamples-nAverage+1;
i++ ) {
272 for( unsigned j=0;j<nAverage;j++ ) {
273 (*pSum)[
i] += mySamples[
i+j];
274 }
275 }
276 }
277 else {
278 float maxADCPeak = 0.;
279 unsigned int iPeakSamp = 0;
280 float averagedADC = 0;
282 if ( maxADCPeak < mySamples[i] )
283 {
284 maxADCPeak = mySamples[
i];
286 }
287 if ( (int)i >= nMin && i < nMin+nAverage )
288 averagedADC += mySamples[
i];
289 }
290 averagedADC /= myScale;
291
292 bool CubicFailed = false;
293 float ADCPeak=0.;
295
296
297
299 {
300 ADCPeak = maxADCPeak;
301 time = ((
float)iPeakSamp) * 25.0 * nanosecond;
302 }
303
304
305
306 else {
307
308
309 if(isFCAL &&
m_mode ==
"CUBIC" &&
311
313
314 unsigned int it0;
315
316 const float invT[3][3]
317 = { { 1, 0, 0 },
318 { -1.5, 2, -0.5 },
319 { 0.5, -1, 0.5 } };
320
321
322 if ( iPeakSamp <= 1 ) {
323 it0 = 1;
324 } else if ( iPeakSamp >= nSamples - 1 ) {
326 } else {
327 it0 = iPeakSamp - 1;
328 }
329
330
331
332 float A[3] = {0, 0, 0};
333 float dtmax = 0.0;
334
335 for (int ia = 0; ia < 3; ia++)
336 for (
int it = 0;
it < 3;
it++)
337 A[ia] += invT[ia][it] * mySamples[it0+it];
338
339
340 if ( not ( CubicFailed = ( A[2] == 0 ) ) ) {
341 dtmax = -1.0 *
A[1] / 2.0 /
A[2];
342 if ( ! ( CubicFailed = ( dtmax < 0 || dtmax > 2 ) ) ) {
343
345 for(int ia = 0; ia < 3; ia++)
346 ADCPeak += A[ia] *
pow(dtmax, ia);
347 }
348 }
349
350
352 float weightSum = 0.;
353 float timeSum = 0.;
354 for(
int it=0;
it<3;
it++) {
355 timeSum +=
float(mySamples[it0+it]) *
float(it);
356 weightSum +=
float(mySamples[it0+it]);
357 }
358 time = (
float(it0) + timeSum / weightSum) * 25.0 * nanosecond;
359
361 CubicFailed = false;
362 }
363
364
365 }
else if(not isFCAL ||
m_mode !=
"CUBIC" ) {
366
367
368 if (
m_mode ==
"PARABOLA" &&
372 try {
373 const Identifier
id =
cabling->cnvToIdentifier(chid);
375 if (
m_emId->is_em_barrel(
id)) {
377 }
378 }
379 catch (LArID_Exception & execpt){
381 << " Cannot get offline identifier from online ID = " << chid );
382 }
383
385 if(peak.size() >1){
386 ADCPeak = peak[0]-thePedestal;
387 if (peak.size()==2)
time = peak[1];
389 }else{
390 ATH_MSG_DEBUG(
"No pic is computed from Parabola. Use scaled average of selected samples" );
391 ADCPeak = averagedADC;
393 }
394 }else{
395 ATH_MSG_FATAL(
"No parabola tool available ! Choose another mode" );
396 }
397 }
398
399
400 else if (
m_mode ==
"CUBIC" &&
402
403 unsigned int it0;
404 const float invT[4][4]
405 = { { 1, 0, 0, 0},
406 { -1.83333, 3, -1.5, 0.333333},
407 { 1, -2.5, 2, -0.5},
408 {-0.166666, 0.5, -0.5, 0.166666} };
409
410
411 if ( iPeakSamp <= 1 ) {
413 } else if ( iPeakSamp >= nSamples - 2 ) {
415 } else {
416 it0 = ( mySamples[iPeakSamp-2] > mySamples[iPeakSamp+2] )
417 ? iPeakSamp - 2
418 : iPeakSamp - 1;
419 }
420
421
422 float A[4] = {0, 0, 0, 0};
423 float dtmax = 0.0;
424 float disc;
425
426 for (int ia = 0; ia < 4; ia++)
427 for (
int it = 0;
it < 4;
it++)
428 A[ia] += invT[ia][it] * mySamples[it0+it];
429
430
431 disc =
A[2]*
A[2] - 3*
A[1]*
A[3];
432 if ( ! ( CubicFailed = ( disc < 0 || A[3] == 0 ) ) ) {
433 dtmax = (-
A[2]-std::sqrt(disc))/(A[3]*3);
434 if ( ! ( CubicFailed = ( dtmax < 0 || dtmax > 3 ) ) ) {
435 time = (
float(it0) + dtmax) * 25.0 * nanosecond;
436 for(int ia = 0; ia < 4; ia++)
437 ADCPeak += A[ia] *
pow(dtmax, ia);
438 }
439 }
440 }
441 }
442 }
443
444
445 if(
m_mode ==
"FIXED" || CubicFailed ) {
446 ADCPeak = averagedADC;
448 }
449
451 << CubicFailed
452 << " Detector: "
453 << (isEM?"EM":(isHEC?"HEC":(isFCAL?"FCAL":"none")))
454 << " Mode: "
456 << " Signal: "
457 << ADCPeak
458 << " Peak: "
459 << maxADCPeak
460 << " PeakSample: "
461 << iPeakSamp
462 );
463
466 float ADCPeakPower=ADCPeak;
467
468 LArVectorProxy ramp = adc2mev->
ADC2MEV(chid,gain);
469
470 if (ramp.size()>1 && ramp[1]<500 && ramp[1]>0) {
472 for (
unsigned i=1;
i<ramp.size();
i++)
474
475 ADCPeakPower*=ADCPeak;
476 }
477 }
478 }
479 if (energy==-9999) {
480 ATH_MSG_DEBUG(
"No Ramp found for this cell. Use default values" );
481
482
483 float ADCtoMeV = 10.0;
484
485
486
487
488
489
490 try {
491 const Identifier
id =
cabling->cnvToIdentifier(chid);
492
494
495 if (
m_emId->is_em_barrel(
id)) {
499 if (layer==2 &&
eta<32)
500 ADCtoMeV=ADCtoMeV*(12./18.);
501 ATH_MSG_DEBUG(
" in EMB s="<< layer <<
", using ADCtoMeV = " << ADCtoMeV );
502 }
else if (
m_emId->is_em_endcap_inner(
id)) {
503
505 ATH_MSG_DEBUG(
" in EMEC inner s="<<
m_emId->sampling(
id)<<
", using ADCtoMeV = " << ADCtoMeV );
506 }
else if (
m_emId->is_em_endcap_outer(
id)) {
507
508
510 ATH_MSG_DEBUG(
" in EMEC outer s="<<
m_emId->sampling(
id)<<
", using ADCtoMeV = " << ADCtoMeV );
511 }
else if (
m_fcalId->is_lar_fcal(
id)) {
512
515 }
else if (
m_hecId->is_lar_hec(
id)) {
516
519 }
520 }
521
522 catch (LArID_Exception & execpt)
523 {
525 ADCtoMeV=1;
526 }
527 energy = ADCPeak * ADCtoMeV * GainFactor;
528 }
531
533
534 LArRawChannel larRawChannel(chid,(int)energy,(int)time,iquality,iprovenance, gain);
535 larRawChannelContainer->add(larRawChannel);
536 }
537 }
538 if ( iloop == 0 ) {
539
541 float tmpSum = 0;
542
543 if ( !fSumEM.empty() ) {
544 for(i=0;
i<fSumEM.size();
i++) {
545 if (i == 0 || fSumEM[i] > tmpSum ) {
548 }
549 }
550 ATH_MSG_DEBUG(
"Found best EM window starting at sample <" << nMinEM <<
">" );
551 }
552
553 for(i=0;
i<fSumHEC.size();
i++) {
554 if (i == 0 || fSumHEC[i] > tmpSum ) {
557 }
558 }
559 if ( !fSumHEC.empty() ) {
560 ATH_MSG_DEBUG(
"Found best HEC window starting at sample <" << nMinHEC <<
">" );
561 }
562
563 for(i=0;
i<fSumFCAL.size();
i++) {
564 if (i == 0 || fSumFCAL[i] > tmpSum ) {
566 tmpSum = fSumFCAL[
i];
567 }
568 }
569 if ( !fSumFCAL.empty() ) {
570 ATH_MSG_DEBUG(
"Found best FCAL window starting at sample <" << nMinFCAL <<
">" );
571 }
572
573 }
574 }
575
577
578
580
581 return StatusCode::SUCCESS;
582}
Scalar eta() const
pseudorapidity method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_WARNING(x)
constexpr int pow(int base, int exp) noexcept
const ServiceHandle< StoreGateSvc > & detStore() const
virtual float pedestal(const HWIdentifier &id, int gain) const =0
const LArVectorProxy ADC2MEV(const HWIdentifier &id, int gain) const
SG::ReadHandleKey< LArDigitContainer > m_DataLocation
SG::ReadCondHandleKey< LArADC2MeV > m_adc2mevKey
SG::ReadCondHandleKey< LArOnOffIdMapping > m_cablingKey
SG::WriteHandleKey< LArRawChannelContainer > m_ChannelContainerName
time(flags, cells_name, *args, **kw)
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())
retrieve(aClass, aKey=None)
bool isFCAL(const xAOD::CaloCluster *cluster)
return true if the cluster (or the majority of its energy) is in the FCAL0