PFSubtractionEnergyRatioCalculator.cxx 3.65 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"

#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"

#include "AthenaKernel/getMessageSvc.h"

PFSubtractionEnergyRatioCalculator::PFSubtractionEnergyRatioCalculator() : AthMessaging(Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"),
                  "PFSubtractionEnergyRatioCalculator"){}

void PFSubtractionEnergyRatioCalculator::calculateSubtractedEnergyRatios(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterSubtractionList, std::map<xAOD::CaloCluster*, double>& clusterEnergyMap, std::vector<float>& clusterSubtractedEnergyRatios){

  ATH_MSG_DEBUG("Setting subtracted energy ratios here");
  
  for (auto thisCluster: clusterSubtractionList) {
    ATH_MSG_DEBUG("Cluster energies are " << thisCluster.first->e() << " and " << clusterEnergyMap[thisCluster.first]);

    //clusterEnergyMap[thisCluster.first can be zero, but this is only a problem if thisCluster.first.e() < 0 - this can happen if a cluster starts with E =0 
    //from a previous shower subtraction step and then we subtract more such that the new E is < 0. Then the ratio would cause an FPE due to the zero.
    //If both thisCluster.first->e() and clusterEnergyMap[thisCluster.first] are zero we never enter this step.          

    if (std::abs(thisCluster.first->e() - clusterEnergyMap[thisCluster.first]) > 0.0001) {
       if ( clusterEnergyMap[thisCluster.first] >= 0) {
        ATH_MSG_DEBUG("Subtracted energy ratio is " << thisCluster.first->e()/clusterEnergyMap[thisCluster.first]);
        clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/clusterEnergyMap[thisCluster.first]);
       }
       //approximate zero with 0.0001 to avoid FPE and still give a meaningful ratio (e.g -100/0.0001)
       else {
         ATH_MSG_DEBUG("Subtracted energy ratio is " << thisCluster.first->e()/0.0001);
         clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/0.0001);
       }
    }
    //else if the cluster enegry did not change then we use NAN to denote that no charged shower subtraction occurred.
    else {
      clusterSubtractedEnergyRatios.push_back(NAN);
      ATH_MSG_DEBUG("Subtracted energy ratio is NAN ");
    }
  }//Loop over clusterSubtractionList

}

void PFSubtractionEnergyRatioCalculator::calculateSubtractedEnergyRatiosForAnnih(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterSubtractionList, std::map<xAOD::CaloCluster*, double>& clusterEnergyMap, std::vector<float>& clusterSubtractedEnergyRatios){

  ATH_MSG_DEBUG("Setting subtracted energy ratios for annih here");

  for (auto thisCluster: clusterSubtractionList) {
    ATH_MSG_DEBUG("Cluster energies are " << thisCluster.first->e() << " and " << clusterEnergyMap[thisCluster.first]);

    //The energy can be zero if we previously "annihiliated" this cluster because we will have already set the energy to zero.
    //We don't want to set the energy ratio to zero if we just change the energy from "0" to "0" - in effect in this case
    //we don't subtract anything and we will set the value of the ratio to NAN which denotes that no subtraction was performed.
    //So we build a list of indices of clusters which don't need to have the energy set to zero - we compare a floating point
    //directly to zero, because in a previous annihilation we explicity set the energy to be exactly zero.
    if (0 != clusterEnergyMap[thisCluster.first]) {
      ATH_MSG_DEBUG("Setting cluster energy ratio to zero");
      clusterSubtractedEnergyRatios.push_back(0);
    }
    else {
      ATH_MSG_DEBUG("Setting cluster energy ratio to NAN");
      clusterSubtractedEnergyRatios.push_back(NAN);
    }

  }

}