this method is purely virtual because every derived class needs to implement it.
129{
130 CaloLCCoeffHelper hp;
132
133 float wrong_dm_energy = 0.0;
134 float dm_total = 0.0;
135
136
137 double weightMom (1);
140 }
141
142
143
144
145 double eWeighted = theCluster->e();
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
164#ifdef MAKE_MOMENTS
165 set_zero_moments(theCluster);
166#endif
167 return StatusCode::SUCCESS;
168 }
169
170 double pi0Prob = 0;
173 ATH_MSG_ERROR(
"Cannot retrieve EM_PROBABILITY cluster moment, "
174 << " cluster energy " << theCluster->e() << " remains the same." );
175 return StatusCode::FAILURE;
176 }
177 if ( pi0Prob < 0 ) {
178 pi0Prob = 0;
179 } else if ( pi0Prob > 1 ) {
180 pi0Prob = 1;
181 }
183 pi0Prob = 1;
184 }
185
186 double center_lambda = 0;
188 ATH_MSG_ERROR(
"Cannot retrieve CENTER_LAMBDA cluster moment, "
189 << " cluster energy " << theCluster->e() << " remains the same." );
190 return StatusCode::FAILURE;
191 }
192
193 const CaloLocalHadCoeff*
data(
nullptr);
194 SG::ReadCondHandle<CaloLocalHadCoeff> rch(
m_key,ctx);
198 return StatusCode::FAILURE;
199 }
200
201 ATH_MSG_DEBUG(
"Cluster is selected for local DM calibration."
202 << " Old cluster energy:" << theCluster->e()
204
205
206 std::vector<Area> areas;
207 std::vector<Cell>
cells;
208 float smp_energy[CaloSampling::Unknown];
209 float cls_unweighted_energy = 0;
211 smp_energy, cls_unweighted_energy,
data);
212 if ( !
sc.isSuccess() ) {
213#ifdef MAKE_MOMENTS
214 set_zero_moments(theCluster);
215#endif
218 }
219
220 float cls_initial_energy = theCluster->e();
221 float cls_side = (theCluster->eta()<0?-1.0:1.0);
222 float cls_eta = fabs(theCluster->eta());
223 float cls_phi = theCluster->phi();
224 float log10ener = 0.0;
225 if(cls_unweighted_energy > 0.0) {
226 log10ener = log10(cls_unweighted_energy);
227 }
228 float log10lambda;
229 if(center_lambda > 0.0) {
230 log10lambda = log10(center_lambda);
231 if(log10lambda >=4.0) log10lambda = 3.9999;
232 }else{
233 log10lambda = 0.0;
234 }
235#ifdef DEBUG_DMTHINGS
236 std::cout << " cls_initial_energy: " << cls_initial_energy << " cls_eta: " << cls_eta << std::endl;
237 std::cout << " log10ener: " << log10ener << " log10lambda: " << log10lambda << std::endl;
238 for(int i_smp=0; i_smp<CaloSampling::Unknown; i_smp++){
239 std::cout << " i_smp: " << i_smp << " smp_energy: " << smp_energy[i_smp] << std::endl;
240 }
241#endif
242
243 std::vector<float > vars;
244 vars.resize(6, 0.0);
251
252 size_t n_dm =
data->getSizeAreaSet();
253
254
255 for(int i_mix=0; i_mix<2; i_mix++){
256 float mixWeight = (i_mix==0?(1-pi0Prob):pi0Prob);
257 if(mixWeight == 0) continue;
258
260
261
262 for(size_t i_dm=0; i_dm < n_dm; i_dm++){
263 if(areas[i_dm].eprep <= 0.0) continue;
264
265 const CaloLocalHadCoeff::LocalHadArea *
area =
data->getArea(i_dm);
267 if(log10ener > emax) log10ener = emax - 0.0001;
269
271 if( !pars ) continue;
272
273 float edm = 0.0;
274
279
281 }
282
283
289 }
290
291
293 const int nSmp=
pars->size()-1;
294
295 double ecalonew = 0.0;
296 double ecaloold = 0.0;
297 for(int i_smp=0; i_smp<nSmp; i_smp++){
298 ecaloold += smp_energy[i_smp];
299 ecalonew += smp_energy[i_smp] * (*pars)[i_smp];
300 }
301 ecalonew += (*pars)[nSmp];
302 edm = ecalonew - ecaloold;
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318 } else{
319 std::cout <<
"CaloLCDeadMaterialTool -> Error! Unknown correction type " <<
area->getType()
320 <<
" with number of parameters " <<
area->getNpars() << std::endl;
321 }
322 wrong_dm_energy += areas[i_dm].edm_wrong;
323
324#ifdef DEBUG_DMTHINGS
325 std::cout << " result> edm: " << edm << " edm_wrong:" << edm_wrong << std::endl;
326#endif
327 edm -= areas[i_dm].edm_wrong;
328 if(edm > 0.0) {
329 areas[i_dm].erec = areas[i_dm].erec + edm*mixWeight;
330 dm_total += edm;
331 }
332 }
333 }
334
335
336
337
338
339
340#ifdef DEBUG_DMTHINGS
341 std::cout << "wc> calculation of weights" << std::endl;
342#endif
343
344 if(dm_total > 0.0) {
346
347 theCluster->setE(cls_initial_energy + dm_total - wrong_dm_energy);
348
350
351
352 float sub_ener_old = 0.0;
353 for(
size_t i_c = 0; i_c <
cells.size(); i_c++)
354 {
357 sub_ener_old +=
cell.weight*
cell.energy;
358 }
359 if(sub_ener_old > 0.0) {
360 float corr = (sub_ener_old + dm_total)/sub_ener_old;
361 for(
size_t i_c = 0; i_c <
cells.size(); i_c++)
362 {
364 cells[i_c].weight *= corr;
365 }
366 }
367
369
370
371
372 for (
size_t cell_index = 0; cell_index <
cells.size(); ++cell_index)
373 {
376 float we =
cell.weight *
cell.energy;
378 areas[
cell.dm].sub_ener_old += we;
382 }
383
384 for(size_t i_dm=0; i_dm < areas.size(); i_dm++) {
385 Area& dma = areas[i_dm];
386 float corrfac = 1;
387 if (i_dm >= n_dm || dma.eprep <= 0)
388 corrfac = 1;
389 else if (dma.sub_ener_old > 0)
390 corrfac = (dma.sub_ener_old + dma.erec) / dma.sub_ener_old;
391
393 areas[
sDM_FCAL].sub_ener_old += (corrfac-1)*dma.sub_ener_old;
395 areas[
sDM_LEAKAGE].sub_ener_old += (corrfac-1)*dma.sub_ener_old;
397 areas[
sDM_UNCLASS].sub_ener_old += (corrfac-1)*dma.sub_ener_old;
398
399 dma.corr = corrfac;
400 }
401
403 for (
size_t cell_index = 0; cell_index <
cells.size(); ++cell_index)
404 {
409 cell.weight *= gblfac;
410 }
411 } else {
413 }
414
415
417 int cell_index = 0;
420 for (;itrCell!=itrCellEnd; ++itrCell) {
421
422 const float old_weight = itrCell.weight();
423 float new_weight =
cells[cell_index].weight;
425 theCluster->reweightCell(itrCell, new_weight);
426 }else{
427 new_weight = old_weight;
428 }
429#ifdef DEBUG_DMTHINGS
430 std::cout <<
" cells> " << cell_index <<
" ener: " <<
cells[cell_index].energy <<
" old_w: " << old_weight <<
" new_w:" << new_weight << std::endl;
431 if( new_weight > 100) std::cout << "DeadMaterialCorrectionTool2 -> Panic! Too large weight" << std::endl;
432#endif
433 cell_index++;
434 }
436 }
437
438 }
439
441 << " (contains wrongly estimated DM energy: " << wrong_dm_energy << ")"
442 << " calculated DM energy (to be added):" << dm_total
443 << " new cluster energy:" << theCluster->e());
444
445#ifdef MAKE_MOMENTS
446
448 double xmom;
449 bool result = theCluster->retrieveMoment(xtype, xmom,
false);
450
451 if(result) {
453 bool useLink = true;
454
455 x = (
double)smp_energy[CaloSampling::PreSamplerB];
456 theCluster->insertMoment(CaloCluster::ENG_RECO_EMB0,
x);
457
458 x = (
double)smp_energy[CaloSampling::PreSamplerE];
459 theCluster->insertMoment(CaloCluster::ENG_RECO_EME0,
x);
460
461 x = (
double)smp_energy[CaloSampling::TileGap3];
462 theCluster->insertMoment(CaloCluster::ENG_RECO_TILEG3,
x);
463
464 x = (
double)(dm_total+wrong_dm_energy);
465 theCluster->insertMoment(CaloCluster::ENG_RECO_DEAD_TOT,
x);
466
467 for(size_t i_dm=0; i_dm < areas.size(); i_dm++) {
468 Area& dma = areas[i_dm];
471 if(i_dm ==
sDM_EMB0_EMB1 && smp_energy[CaloSampling::PreSamplerB]>0.0)
x+= smp_energy[CaloSampling::PreSamplerB];
472 if(i_dm ==
sDM_EME0_EME12 && smp_energy[CaloSampling::PreSamplerE] > 0.0)
x+= smp_energy[CaloSampling::PreSamplerE];
473 if(i_dm ==
sDM_SCN && smp_energy[CaloSampling::TileGap3] > 0.0)
x+= smp_energy[CaloSampling::TileGap3];
474 theCluster->insertMoment(m_type, CaloClusterMoment(
x));
475 }
476 }
477
478#endif
479
480
481
482
483
484
485 if ( eWeighted > 0 || eWeighted < 0 ) {
486 weightMom *= theCluster->e()/eWeighted;
487 }
489
490 return StatusCode::SUCCESS;
491}
static void calculateKine(xAOD::CaloCluster *clu, const bool useweight=true, const bool updateLayers=true, const bool useGPUCriteria=false)
Helper class to calculate cluster kinematics based on cells.
static bool Interpolate(const CaloLocalHadCoeff *m_data, const unsigned int n_area, std::vector< float > &x, CaloLocalHadCoeff::LocalHadCoeff &pars, const std::vector< int > &dim, double xfit=0.)
std::vector< float > LocalHadCoeff
Correction parameters for one general bin.
CaloClusterCellLink::iterator cell_iterator
Iterator of the underlying CaloClusterCellLink (non-const version)
MomentType
Enums to identify different moments.
@ EM_PROBABILITY
Classification probability to be em-like.
@ DM_WEIGHT
Dead-material weight (E_dm/E_ooc)
@ CENTER_LAMBDA
Shower depth at Cluster Centroid.
@ ENG_CALIB_DEAD_TOT
Attached Calibration Hit energy in dead material.
::StatusCode StatusCode
StatusCode definition for legacy code.