281 ATH_MSG_DEBUG(
"Starting CaloCalibClusterMomentsMaker2::execute");
285 bool foundAllContainers (
true);
286 std::vector<const CaloCalibrationHitContainer *> v_cchc;
295 msg(MSG::ERROR) <<
"SG does not contain calibration hit container " << key.key() <<
endmsg;
297 foundAllContainers =
false;
300 v_cchc.push_back(cchc.
cptr());
304 std::vector<const CaloCalibrationHitContainer *> v_dmcchc;
313 ATH_MSG_ERROR(
"SG does not contain DM calibration hit container " << key.key());
315 foundAllContainers =
false;
318 v_dmcchc.push_back(cchc.
cptr());
325 if ( !foundAllContainers )
return StatusCode::SUCCESS;
340 for( ;clusIter!=clusIterEnd;++clusIter,++iClus) {
346 for(; cellIter != cellIterEnd; cellIter++ ){
347 const CaloCell* pCell = (*cellIter);
351 CellInfoSet_t::iterator bookmark = cellInfo.lower_bound( myId );
352 if(bookmark == cellInfo.end() || bookmark->first != myId) {
354 if (bookmark != cellInfo.begin()) --bookmark;
355 cellInfo.emplace_hint (bookmark, myId, std::move(info));
358 bookmark->second.Add( info );
367 unsigned int nHitsTotal = 0;
368 unsigned int nHitsWithoutParticleUID = 0;
374 CellInfoSet_t::iterator pos = cellInfo.find( myId );
375 if(pos != cellInfo.end() ) {
378 for (
const std::pair<int, double>& p : pos->second) {
380 double weight = p.second;
384 ATH_MSG_ERROR(
"Invalid uniqueID detected - this sample cannot be properly analysed.");
387 clusInfoVec[iClus].Add(weight * hit->energyTotal(), nsmp, (
unsigned int)(
HepMC::uniqueID(hit)));
389 clusInfoVec[iClus].Add(weight * hit->energyTotal(), nsmp);
394 nHitsWithoutParticleUID++;
404 ATH_MSG_INFO(
"Calibration hits do not have ParticleUID, ids of particle-caused hits are always 0. Continuing without ParticleID machinery.");
405 useParticleID =
false;
412 if (doCalibFrac && !truthParticleContainerReadHandle.
isValid()){
417 std::vector<double> engCalibOut[3];
422 std::vector<std::vector <int > > clusListL;
423 std::vector<std::vector <int > > clusListM;
424 std::vector<std::vector <int > > clusListT;
432 for(
unsigned int ii=0;ii<3;ii++) {
436 engCalibOut[ii].resize(theClusColl->
size(),0);
444 if ( theCluster->eta() < 0 ) iEtaSign = -1;
450 unsigned int iEtaBin(jeta);
451 if ( jeta < 0 ) iEtaBin = abs(jeta)-1;
452 const std::vector<CalibHitIPhiIEtaRange>&
bins
455 int jp = range.iPhi+jphi;
458 int jEtaMin = iEtaSign<0?-range.iEtaMax-1:range.iEtaMin;
459 int jEtaMax = iEtaSign<0?-range.iEtaMin-1:range.iEtaMax;
460 for(
int je = jEtaMin;je<=jEtaMax;je++ ) {
481 CellInfoSet_t::iterator pos = cellInfo.find( myId );
482 if(pos == cellInfo.end() ) {
494 for (
unsigned int ii=0;ii<3;ii++) {
495 std::vector<std::vector <int > > *pClusList=
nullptr;
497 pClusList = &clusListL;
499 pClusList = &clusListM;
501 pClusList = &clusListT;
504 double hitClusNorm(0.);
505 std::vector<int> hitClusIndex;
506 std::vector<double > hitClusEffEnergy;
517 hitClusNorm += pos->second.engTot;
518 hitClusEffEnergy.push_back(pos->second.engTot);
519 hitClusIndex.push_back(iClus);
522 if ( hitClusNorm > 0 ) {
523 const double inv_hitClusNorm = 1. / hitClusNorm;
524 for(
unsigned int i_cls=0; i_cls<hitClusIndex.size(); i_cls++){
525 int iClus = hitClusIndex[i_cls];
526 double w = hitClusEffEnergy[i_cls] * inv_hitClusNorm;
527 engCalibOut[ii][iClus] += w*hit->energyTotal();
544 std::vector<std::vector <int > > *pClusList=
nullptr;
546 pClusList = &clusListL;
548 pClusList = &clusListM;
550 pClusList = &clusListT;
575 int hitClusVecIndex = 0;
576 std::vector<int> hitClusIndex(theClusColl->
size());
577 std::vector<double > hitClusEffEnergy(theClusColl->
size());
578 double hitClusNorm = 0.0;
589 double engClusTruthUniqueIDCalib = pos->second.engTot;
592 double sum_smp_energy = 0.0;
599 sum_smp_energy += pos->second.engSmp[nsmp];
602 if(sum_smp_energy > 0.0) {
603 double phi_diff=myCDDE->
phi()-theCluster->
phi();
604 if(phi_diff <= -
M_PI){
606 }
else if (phi_diff >
M_PI){
609 double eta_diff = (myCDDE->
eta()-theCluster->
eta());
610 float distance=sqrt(eta_diff * eta_diff + phi_diff * phi_diff);
614 hitClusIndex [hitClusVecIndex] = iClus;
615 hitClusEffEnergy [hitClusVecIndex++]= effEner;
616 hitClusNorm += effEner;
621 hitClusIndex.resize(hitClusVecIndex);
622 hitClusEffEnergy.resize(hitClusVecIndex);
626 if(hitClusNorm >0.0) {
627 const double inv_hitClusNorm = 1. / hitClusNorm;
628 for(
unsigned int i_cls=0; i_cls<hitClusIndex.size(); i_cls++){
629 int iClus = hitClusIndex[i_cls];
630 double dm_weight = hitClusEffEnergy[i_cls] * inv_hitClusNorm;
631 if(dm_weight > 1.0 || dm_weight < 0.0 ){
632 std::cout <<
"CaloCalibClusterMomentsMaker2::execute() ->Error! Strange weight " << dm_weight<< std::endl;
633 std::cout << hitClusEffEnergy[i_cls] <<
" " << hitClusNorm << std::endl;
635 clusInfoVec[iClus].engCalibDeadInArea[nDmArea] += hit->energyTotal()*dm_weight;
649 std::vector<double> engCalibFrac;
652 std::map<unsigned int,int> truthIDToPdgCodeMap;
655 for (
const auto *thisTruthParticle : *truthParticleContainerReadHandle){
657 if (!thisTruthParticle){
662 truthIDToPdgCodeMap[
HepMC::uniqueID(thisTruthParticle)] = thisTruthParticle->pdgId();
666 for( clusIter = theClusColl->
begin(),iClus=0; clusIter!=clusIterEnd;++clusIter,++iClus) {
696 double eng_calib_dead_unclass = eng_calib_dead_tot - eng_calib_dead_emb0 - eng_calib_dead_tile0
697 - eng_calib_dead_tileg3 - eng_calib_dead_eme0 - eng_calib_dead_hec0 - eng_calib_dead_fcal
698 - eng_calib_dead_leakage;
709 if (
auto it = truthIDToPdgCodeMap.find(p.first); it != truthIDToPdgCodeMap.end()) {
712 ATH_MSG_WARNING(
"truthIDToPdgCodeMap cannot find an entry with uniqueID " << p.first);
715 if( std::abs(pdg_id) == 211) {
717 }
else if( pdg_id == 111 || pdg_id == 22 || std::abs(pdg_id)==11) {
723 for(
size_t i=0; i<engCalibFrac.size(); i++) engCalibFrac[i] = engCalibFrac[i]/clusInfo.
engCalibIn.
engTot;
730 moment_name_set::const_iterator vMomentsIter =
m_validMoments.begin();
731 moment_name_set::const_iterator vMomentsIterEnd =
m_validMoments.end();
734 for(; vMomentsIter!=vMomentsIterEnd; ++vMomentsIter,++iMoment) {
736 switch (vMomentsIter->second) {
742 myMoments[iMoment] = engCalibOut[0][iClus];
745 myMoments[iMoment] = engCalibOut[1][iClus];
748 myMoments[iMoment] = engCalibOut[2][iClus];
760 myMoments[iMoment] = eng_calib_dead_tot;
763 myMoments[iMoment] = eng_calib_dead_emb0;
766 myMoments[iMoment] = eng_calib_dead_tile0;
769 myMoments[iMoment] = eng_calib_dead_tileg3;
772 myMoments[iMoment] = eng_calib_dead_eme0;
775 myMoments[iMoment] = eng_calib_dead_hec0;
778 myMoments[iMoment] = eng_calib_dead_fcal;
781 myMoments[iMoment] = eng_calib_dead_leakage;
784 myMoments[iMoment] = eng_calib_dead_unclass;
800 theCluster->
insertMoment(vMomentsIter->second, myMoments[iMoment]);
805 return StatusCode::SUCCESS;